MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_afld.t
Go to the documentation of this file.
1 !> Module for including anisotropic flux limited diffusion (AFLD)-approximation in Radiation-hydrodynamics simulations using mod_rhd
2 !> Based on Turner and stone 2001. See
3 !> [1]Moens, N., Sundqvist, J. O., El Mellah, I., Poniatowski, L., Teunissen, J., and Keppens, R.,
4 !> “Radiation-hydrodynamics with MPI-AMRVAC . Flux-limited diffusion”,
5 !> <i>Astronomy and Astrophysics</i>, vol. 657, 2022. doi:10.1051/0004-6361/202141023.
6 !> For more information.
7 
8 module mod_afld
9 
10  use mod_comm_lib, only: mpistop
11 
12  implicit none
13 
14  !> source split for energy interact and radforce:
15  logical :: fld_eint_split = .false.
16  logical :: fld_radforce_split = .false.
17 
18  !> Opacity per unit of unit_density
19  double precision, public :: fld_kappa0 = 0.34d0
20 
21  !> mean particle mass
22  double precision, public :: afld_mu = 0.6d0
23 
24  !> Tolerance for bisection method for Energy sourceterms
25  !> This is a percentage of the minimum of gas- and radiation energy
26  double precision, public :: fld_bisect_tol = 1.d-4
27 
28  !> Tolerance for adi method for radiative Energy diffusion
29  double precision, public :: fld_diff_tol = 1.d-4
30 
31  !> Number for splitting the diffusion module
32  double precision, public :: diff_crit
33 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DELETE WHEN DONE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
34  !> Index for testvariable
35  !!
36  !!
37  !!
38  integer, public :: i_test
39 
40  !> Use constant Opacity?
41  character(len=8), allocatable :: fld_opacity_law(:)
42  character(len=50) :: fld_opal_table = 'Y09800'
43 
44  !> Diffusion limit lambda = 0.33
45  character(len=16) :: fld_fluxlimiter = 'Pomraning'
46 
47  ! !> diffusion coefficient for multigrid method
48  ! integer :: i_diff_mg
49  !> Index for Diffusion coeficients
50  integer, allocatable, public :: i_diff_mg(:)
51 
52  !> Which method to solve diffusion part
53  character(len=8) :: afld_diff_scheme = 'mg'
54 
55  !> Which method to find the root for the energy interaction polynomial
56  character(len=8) :: fld_interaction_method = 'Halley'
57 
58  !> Set Diffusion coefficient to unity
59  logical :: fld_diff_testcase = .false.
60 
61  !> Take running average for Diffusion coefficient
62  logical :: diff_coef_filter = .false.
63  integer :: size_d_filter = 1
64 
65  !> Take a running average over the fluxlimiter
66  logical :: flux_lim_filter = .false.
67  integer :: size_l_filter = 1
68 
69  !> Use or dont use lineforce opacities
70  logical :: lineforce_opacities = .false.
71 
72  !> Resume run when multigrid returns error
73  logical :: diffcrash_resume = .true.
74 
75  !> A copy of rhd_Gamma
76  double precision, private, protected :: afld_gamma
77 
78  !> running timestep for diffusion solver, initialised as zero
79  double precision :: dt_diff = 0.d0
80 
81  !> public methods
82  !> these are called in mod_rhd_phys
83  public :: get_afld_rad_force
84  public :: get_afld_energy_interact
85  public :: afld_radforce_get_dt
86  public :: afld_init
87  public :: afld_get_radflux
88  public :: afld_get_radpress
89  public :: afld_get_fluxlimiter
90  public :: afld_get_opacity
91 
92  contains
93 
94 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
95 !!!!!!!!!!!!!!!!!!! GENERAL
96 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
97 
98  !> Reading in fld-list parameters from .par file
99  subroutine afld_params_read(files)
100  use mod_global_parameters, only: unitpar
101  use mod_constants
102  character(len=*), intent(in) :: files(:)
103  integer :: n
104 
105  namelist /fld_list/ fld_kappa0, fld_eint_split, fld_radforce_split, &
110 
111  do n = 1, size(files)
112  open(unitpar, file=trim(files(n)), status="old")
113  read(unitpar, fld_list, end=111)
114  111 close(unitpar)
115  end do
116  end subroutine afld_params_read
117 
118  !> Initialising FLD-module:
119  !> Read opacities
120  !> Initialise Multigrid
121  !> adimensionalise kappa
122  !> Add extra variables to w-array, flux, kappa, eddington Tensor
123  !> Lambda and R
124  !> ...
125  subroutine afld_init(He_abundance, rhd_radiation_diffusion, afld_gamma)
127  use mod_variables
128  use mod_physics
131 
132  double precision, intent(in) :: he_abundance, afld_gamma
133  logical, intent(in) :: rhd_radiation_diffusion
134  double precision :: sigma_thomson
135  integer :: idir,jdir
136 
137  character(len=1) :: ind_1
138  character(len=1) :: ind_2
139  character(len=2) :: cmp_f
140  character(len=5) :: cmp_e
141 
142  {^ifoned
143  call mpistop("Anisotropic in 1D... really?")
144  }
145 
146  !> allocate boundary conditions
147  allocate(fld_opacity_law(1:ndim))
148  fld_opacity_law(1:ndim) = 'const'
149 
150  !> read par files
152 
153  ! !> Set lineforce opacities as variable
154  ! if (lineforce_opacities) then
155  ! allocate(i_opf(ndim))
156  ! do idir = 1,ndim
157  ! write(ind_1,'(I1)') idir
158  ! cmp_f = 'k' // ind_1
159  ! i_opf(idir) = var_set_extravar(cmp_f,cmp_f)
160  ! enddo
161  ! endif
162 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! DELETE WHEN DONE!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
163  !> Introduce test variable globally
164  !!
165  !!
166  !!
167  ! i_test = var_set_extravar('test','test')
168 
169  if (rhd_radiation_diffusion) then
170  if (afld_diff_scheme .eq. 'mg') then
171 
172  use_multigrid = .true.
173 
176 
177  mg%n_extra_vars = 1
178  mg%operator_type = mg_ahelmholtz
179 
180  endif
181  endif
182 
183  allocate(i_diff_mg(ndim))
184  do idir = 1,ndim
185  write(ind_1,'(I1)') idir
186  cmp_f = 'D' // ind_1
187  i_diff_mg(idir) = var_set_extravar(cmp_f,cmp_f)
188  enddo
189 
190  !> Need mean molecular weight
191  afld_mu = (1.+4*he_abundance)/(2.+3.*he_abundance)
192 
193  !> set rhd_gamma
194  ! afld_gamma = phys_gamma
195 
196  ! !> Read in opacity table if necesary
197  ! if (any(fld_opacity_law) .eq. 'opal') call init_opal_table(fld_opal_table)
198  !
199  ! sigma_thomson = 6.6524585d-25
200  ! fld_kappa0 = sigma_thomson/const_mp * (1.+2.*He_abundance)/(1.+4.*He_abundance)
201 
202  end subroutine afld_init
203 
204  !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
205  !> This subroutine handles the radiation force
206  subroutine get_afld_rad_force(qdt,ixI^L,ixO^L,wCT,w,x,&
207  energy,qsourcesplit,active)
208  use mod_constants
210  use mod_usr_methods
211  use mod_geometry
212 
213  integer, intent(in) :: ixi^l, ixo^l
214  double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
215  double precision, intent(in) :: wct(ixi^s,1:nw)
216  double precision, intent(inout) :: w(ixi^s,1:nw)
217  logical, intent(in) :: energy,qsourcesplit
218  logical, intent(inout) :: active
219 
220  double precision :: radiation_forcect(ixo^s,1:ndim)
221  double precision :: kappact(ixo^s,1:ndim)
222  double precision :: rad_fluxct(ixo^s,1:ndim)
223 
224  double precision :: div_v(ixi^s,1:ndim,1:ndim)
225  double precision :: edd(ixo^s,1:ndim,1:ndim)
226  double precision :: nabla_vp(ixo^s)
227  double precision :: vel(ixi^s), grad_v(ixi^s), grad0_v(ixo^s)
228  double precision :: grad_e(ixo^s)
229 
230  integer :: idir, jdir
231 
232  double precision :: fld_r(ixo^s,1:ndim), lambda(ixo^s,1:ndim)
233 
234  !> Calculate and add sourceterms
235  if(qsourcesplit .eqv. fld_radforce_split) then
236  active = .true.
237 
238  call afld_get_opacity(wct, x, ixi^l, ixo^l, kappact)
239  call afld_get_radflux(wct, x, ixi^l, ixo^l, rad_fluxct)
240  call afld_get_fluxlimiter(wct, x, ixi^l, ixo^l, lambda, fld_r)
241 
242  do idir = 1,ndim
243  !> Radiation force = kappa*rho/c *Flux = lambda gradE
244  radiation_forcect(ixo^s,idir) = kappact(ixo^s,idir)*rad_fluxct(ixo^s,idir)/(const_c/unit_velocity)
245 
246  ! call gradientO(wCT(ixI^S,iw_r_e),x,ixI^L,ixO^L,idir,grad_E,nghostcells)
247  ! radiation_forceCT(ixO^S,idir) = lambda(ixO^S)*grad_E(ixO^S)
248 
249  !> Momentum equation source term
250  w(ixo^s,iw_mom(idir)) = w(ixo^s,iw_mom(idir)) &
251  + qdt * wct(ixo^s,iw_rho)*radiation_forcect(ixo^s,idir)
252  ! w(ixO^S,iw_mom(idir)) = w(ixO^S,iw_mom(idir)) &
253  ! + qdt *radiation_forceCT(ixO^S,idir)
254 
255  ! if (energy .and. .not. block%e_is_internal) then
256  !> Energy equation source term (kinetic energy)
257  w(ixo^s,iw_e) = w(ixo^s,iw_e) &
258  + qdt * wct(ixo^s,iw_mom(idir))*radiation_forcect(ixo^s,idir)
259  ! w(ixO^S,iw_e) = w(ixO^S,iw_e) &
260  ! + qdt * wCT(ixO^S,iw_mom(idir))/wCT(ixO^S,iw_rho)*radiation_forceCT(ixO^S,idir)
261  ! endif
262  enddo
263 
264  !> Photon tiring
265  !> calculate tensor div_v
266  !> !$OMP PARALLEL DO
267  do idir = 1,ndim
268  do jdir = 1,ndim
269  vel(ixi^s) = wct(ixi^s,iw_mom(jdir))/wct(ixi^s,iw_rho)
270 
271  call gradient(vel,ixi^l,ixo^l,idir,grad_v)
272  div_v(ixo^s,idir,jdir) = grad_v(ixo^s)
273 
274  ! call gradientO(vel,x,ixI^L,ixO^L,idir,grad0_v,nghostcells)
275  ! div_v(ixO^S,idir,jdir) = grad0_v(ixO^S)
276  enddo
277  enddo
278  !> !$OMP END PARALLEL DO
279 
280  call afld_get_eddington(wct, x, ixi^l, ixo^l, edd)
281 
282  !> VARIABLE NAMES DIV ARE ACTUALLY GRADIENTS
283  {^ifoned
284  nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1) &
285  }
286 
287  {^iftwod
288  !>eq 34 Turner and stone (Only 2D)
289  nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1) &
290  + div_v(ixo^s,1,2)*edd(ixo^s,1,2) &
291  + div_v(ixo^s,2,1)*edd(ixo^s,2,1) &
292  + div_v(ixo^s,2,2)*edd(ixo^s,2,2)
293  }
294 
295  {^ifthreed
296  nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1) &
297  + div_v(ixo^s,1,2)*edd(ixo^s,1,2) &
298  + div_v(ixo^s,1,3)*edd(ixo^s,1,3) &
299  + div_v(ixo^s,2,1)*edd(ixo^s,2,1) &
300  + div_v(ixo^s,2,2)*edd(ixo^s,2,2) &
301  + div_v(ixo^s,2,3)*edd(ixo^s,2,3) &
302  + div_v(ixo^s,3,1)*edd(ixo^s,3,1) &
303  + div_v(ixo^s,3,2)*edd(ixo^s,3,2) &
304  + div_v(ixo^s,3,3)*edd(ixo^s,3,3)
305  }
306 
307  w(ixo^s,iw_r_e) = w(ixo^s,iw_r_e) &
308  - qdt * nabla_vp(ixo^s)*wct(ixo^s,iw_r_e)
309 
310  end if
311 
312  end subroutine get_afld_rad_force
313 
314 
315  subroutine afld_radforce_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
317  use mod_usr_methods
318 
319  integer, intent(in) :: ixi^l, ixo^l
320  double precision, intent(in) :: dx^d, x(ixi^s,1:ndim), w(ixi^s,1:nw)
321  double precision, intent(inout) :: dtnew
322 
323  double precision :: radiation_force(ixo^s,1:ndim)
324  double precision :: rad_flux(ixo^s,1:ndim)
325  double precision :: kappa(ixo^s,1:ndim)
326  double precision :: dxinv(1:ndim), max_grad
327  integer :: idim
328 
329  ^d&dxinv(^d)=one/dx^d;
330 
331  call afld_get_opacity(w, x, ixi^l, ixo^l, kappa)
332  call afld_get_radflux(w, x, ixi^l, ixo^l, rad_flux)
333 
334  do idim = 1, ndim
335  radiation_force(ixo^s,idim) = w(ixo^s,iw_rho)*kappa(ixo^s,idim)*rad_flux(ixo^s,idim)/(const_c/unit_velocity)
336  max_grad = maxval(abs(radiation_force(ixo^s,idim)))
337  max_grad = max(max_grad, epsilon(1.0d0))
338  dtnew = min(dtnew, courantpar / sqrt(max_grad * dxinv(idim)))
339  end do
340 
341  end subroutine afld_radforce_get_dt
342 
343  !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
344  !> This subroutine handles the energy exchange between gas and radiation
345  subroutine get_afld_energy_interact(qdt,ixI^L,ixO^L,wCT,w,x,&
346  energy,qsourcesplit,active)
347  use mod_constants
349  use mod_usr_methods
350 
351  use mod_physics, only: phys_get_pthermal !needed to get temp
352 
353  integer, intent(in) :: ixi^l, ixo^l
354  double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
355  double precision, intent(in) :: wct(ixi^s,1:nw)
356  double precision, intent(inout) :: w(ixi^s,1:nw)
357  logical, intent(in) :: energy,qsourcesplit
358  logical, intent(inout) :: active
359 
360  !> Calculate and add sourceterms
361  if(qsourcesplit .eqv. fld_eint_split) then
362  active = .true.
363  !> Add energy sourceterms
364  call energy_interaction(w, w, x, ixi^l, ixo^l)
365  end if
366  end subroutine get_afld_energy_interact
367 
368  !> Sets the opacity in the w-array
369  !> by calling mod_opal_opacity
370  subroutine afld_get_opacity(w, x, ixI^L, ixO^L, fld_kappa)
372  use mod_physics, only: phys_get_pthermal
373  use mod_physics, only: phys_get_tgas
374  use mod_usr_methods
375  use mod_opal_opacity
376 
377  integer, intent(in) :: ixi^l, ixo^l
378  double precision, intent(in) :: w(ixi^s, 1:nw)
379  double precision, intent(in) :: x(ixi^s, 1:ndim)
380  double precision, intent(out) :: fld_kappa(ixo^s,1:ndim)
381  double precision :: temp(ixi^s), pth(ixi^s), a2(ixo^s)
382  double precision :: rho0,temp0,n,sigma_b
383  double precision :: akram, bkram
384  double precision :: vth(ixo^s), gradv(ixi^s), eta(ixo^s), t(ixo^s)
385 
386  integer :: i,j,ix^d, idim
387 
388  do idim = 1, ndim
389 
390  select case (fld_opacity_law(idim))
391  case('const')
392  fld_kappa(ixo^s,idim) = fld_kappa0/unit_opacity
393  case('thomson')
394  fld_kappa(ixo^s,idim) = fld_kappa0/unit_opacity
395  case('kramers')
396  rho0 = half !> Take lower value of rho in domain
397  fld_kappa(ixo^s,idim) = fld_kappa0/unit_opacity*((w(ixo^s,iw_rho)/rho0))
398  case('bump')
399  !> Opacity bump
400  rho0 = 0.2d0 !0.5d-1
401  n = 7.d0
402  sigma_b = 2.d-2
403  !fld_kappa(ixO^S) = fld_kappa0/unit_opacity*(one + n*dexp(-((rho0 - w(ixO^S,iw_rho))**two)/rho0))
404  fld_kappa(ixo^s,idim) = fld_kappa0/unit_opacity*(one + n*dexp(-one/sigma_b*(dlog(w(ixo^s,iw_rho)/rho0))**two))
405  case('non_iso')
406  call phys_get_pthermal(w,x,ixi^l,ixo^l,temp)
407  temp(ixo^s)=temp(ixo^s)/w(ixo^s,iw_rho)
408 
409  rho0 = 0.5d0 !> Take lower value of rho in domain
410  temp0 = one
411  n = -7.d0/two
412  fld_kappa(ixo^s,idim) = fld_kappa0/unit_opacity*(w(ixo^s,iw_rho)/rho0)*(temp(ixo^s)/temp0)**n
413  case('fastwind')
414  call phys_get_pthermal(w,x,ixi^l,ixo^l,pth)
415  a2(ixo^s) = pth(ixo^s)/w(ixo^s,iw_rho)*unit_velocity**2.d0
416 
417  akram = 13.1351597305
418  bkram = -4.5182188206
419 
420  fld_kappa(ixo^s,idim) = fld_kappa0/unit_opacity &
421  * (1.d0+10.d0**akram*w(ixo^s,iw_rho)*unit_density*(a2(ixo^s)/1.d12)**bkram)
422 
423  {do ix^d=ixomin^d,ixomax^d\ }
424  !> Hard limit on kappa
425  fld_kappa(ix^d,idim) = min(fld_kappa(ix^d,idim),2.3d0*fld_kappa0/unit_opacity)
426 
427  !> Limit kappa through T
428  ! fld_kappa(ix^D) = fld_kappa0/unit_opacity &
429  ! * (1.d0+10.d0**akram*w(ix^D,iw_rho)*unit_density &
430  ! * (max(a2(ix^D),const_kB*5.9d4/(afld_mu*const_mp))/1.d12)**bkram)
431  {enddo\ }
432 
433  case('opal')
434  call phys_get_tgas(w,x,ixi^l,ixo^l,temp)
435  {do ix^d=ixomin^d,ixomax^d\ }
436  rho0 = w(ix^d,iw_rho)*unit_density
437  temp0 = temp(ix^d)*unit_temperature
438  call set_opal_opacity(rho0,temp0,n)
439  fld_kappa(ix^d,idim) = n/unit_opacity
440  {enddo\ }
441 
442  case('special')
443  if (.not. associated(usr_special_aniso_opacity)) then
444  call mpistop("special opacity not defined")
445  endif
446  call usr_special_aniso_opacity(ixi^l, ixo^l, w, x, fld_kappa(ixo^s,idim),idim)
447 
448  case default
449  call mpistop("Doesn't know opacity law")
450  end select
451 
452  end do
453 
454  end subroutine afld_get_opacity
455 
456  !> Calculate fld flux limiter
457  !> This subroutine calculates flux limiter lambda using the prescription
458  !> stored in fld_fluxlimiter.
459  !> It also calculates the ratio of radiation scaleheight and mean free path
460  subroutine afld_get_fluxlimiter(w, x, ixI^L, ixO^L, fld_lambda, fld_R)
462  use mod_geometry
463  use mod_usr_methods
464 
465  integer, intent(in) :: ixi^l, ixo^l
466  double precision, intent(in) :: w(ixi^s, 1:nw)
467  double precision, intent(in) :: x(ixi^s, 1:ndim)
468  double precision, intent(out) :: fld_r(ixo^s,1:ndim), fld_lambda(ixo^s,1:ndim)
469  double precision :: kappa(ixo^s,1:ndim)
470  double precision :: normgrad2(ixo^s)
471  double precision :: grad_r_e(ixi^s), grad_r_eo(ixo^s)
472  double precision :: rad_e(ixi^s)
473  integer :: idir, ix^d
474 
475  double precision :: tmp_l(ixi^s), filtered_l(ixi^s)
476  integer :: filter, idim
477 
478  select case (fld_fluxlimiter)
479  case('Diffusion')
480  fld_lambda(ixo^s,1:ndim) = 1.d0/3.d0
481  fld_r(ixo^s,1:ndim) = zero
482 
483  case('FreeStream')
484  !> Calculate R everywhere
485  !> |grad E|/(rho kappa E)
486  rad_e(ixi^s) = w(ixi^s, iw_r_e)
487 
488  call afld_get_opacity(w, x, ixi^l, ixo^l, kappa)
489 
490  do idir = 1,ndim
491  call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo,nghostcells)
492  normgrad2(ixo^s) = grad_r_eo(ixo^s)**2
493 
494  ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e)
495  ! normgrad2(ixO^S) = normgrad2(ixO^S) + grad_r_e(ixO^S)**2
496  enddo
497 
498  do idir = 1,ndim
499  fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
500 
501  !> Calculate the flux limiter, lambda
502  fld_lambda(ixo^s,idir) = one/fld_r(ixo^s,idir)
503  end do
504 
505 
506 
507  case('Pomraning')
508  !> Calculate R everywhere
509  !> |grad E|/(rho kappa E)
510  normgrad2(ixo^s) = 0.d0 !smalldouble !*w(ixO^S,iw_r_e)
511 
512  rad_e(ixi^s) = w(ixi^s, iw_r_e)
513 
514  call afld_get_opacity(w, x, ixi^l, ixo^l, kappa)
515 
516  do idir = 1,ndim
517  call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo,nghostcells)
518  normgrad2(ixo^s) = grad_r_eo(ixo^s)**2
519  enddo
520 
521  do idir = 1,ndim
522  ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e)
523  ! normgrad2(ixO^S) = normgrad2(ixO^S) + grad_r_e(ixO^S)**2
524  fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
525 
526  !> Calculate the flux limiter, lambda
527  !> Levermore and Pomraning: lambda = (2 + R)/(6 + 3R + R^2)
528  fld_lambda(ixo^s,idir) = (2.d0+fld_r(ixo^s,idir))/(6.d0+3*fld_r(ixo^s,idir)+fld_r(ixo^s,idir)**2.d0)
529  end do
530 
531 
532  case('Minerbo')
533  !> Calculate R everywhere
534  !> |grad E|/(rho kappa E)
535  normgrad2(ixo^s) = 0.d0 !smalldouble
536 
537  rad_e(ixi^s) = w(ixi^s, iw_r_e)
538 
539  call afld_get_opacity(w, x, ixi^l, ixo^l, kappa)
540 
541  do idir = 1,ndim
542  call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo,nghostcells)
543  normgrad2(ixo^s) = grad_r_eo(ixo^s)**2
544  enddo
545 
546  do idir = 1,ndim
547  ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e)
548  ! normgrad2(ixO^S) = normgrad2(ixO^S) + grad_r_e(ixO^S)**2
549 
550  fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
551 
552  !> Calculate the flux limiter, lambda
553  !> Minerbo:
554  {do ix^d=ixomin^d,ixomax^d\ }
555  if (fld_r(ix^d,idir) .lt. 3.d0/2.d0) then
556  fld_lambda(ix^d,idir) = 2.d0/(3.d0 + dsqrt(9.d0 + 12.d0*fld_r(ix^d,idir)**2.d0))
557  else
558  fld_lambda(ix^d,idir) = 1.d0/(1.d0 + fld_r(ix^d,idir) + dsqrt(1.d0 + 2.d0*fld_r(ix^d,idir)))
559  endif
560  {enddo\}
561  end do
562 
563 
564  case('special')
565  if (.not. associated(usr_special_fluxlimiter)) then
566  call mpistop("special fluxlimiter not defined")
567  endif
568  call usr_special_fluxlimiter(ixi^l, ixo^l, w, x, fld_lambda, fld_r)
569  case default
570  call mpistop('Fluxlimiter unknown')
571  end select
572 
573 
574  if (flux_lim_filter) then
575  if (size_l_filter .lt. 1) call mpistop("D filter of size < 1 makes no sense")
576  if (size_l_filter .gt. nghostcells) call mpistop("D filter of size > nghostcells makes no sense")
577 
578  do idir = 1,ndim
579 
580  tmp_l(ixo^s) = fld_lambda(ixo^s, idir)
581  filtered_l(ixo^s) = zero
582 
583  do filter = 1,size_l_filter
584  {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_l_filter\}
585  do idim = 1,ndim
586  filtered_l(ix^d) = filtered_l(ix^d) &
587  + tmp_l(ix^d+filter*kr(idim,^d)) &
588  + tmp_l(ix^d-filter*kr(idim,^d))
589  enddo
590  {enddo\}
591  enddo
592 
593  {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
594  tmp_l(ix^d) = (tmp_l(ix^d)+filtered_l(ix^d))/(1+2*size_l_filter*ndim)
595  {enddo\}
596 
597  enddo
598 
599  fld_lambda(ixo^s,idir) = tmp_l(ixo^s)
600  endif
601 
602  end subroutine afld_get_fluxlimiter
603 
604  !> Calculate Radiation Flux
605  !> stores radiation flux in w-array
606  subroutine afld_get_radflux(w, x, ixI^L, ixO^L, rad_flux)
608  use mod_geometry
609 
610  integer, intent(in) :: ixi^l, ixo^l
611  double precision, intent(in) :: w(ixi^s, 1:nw)
612  double precision, intent(in) :: x(ixi^s, 1:ndim)
613  double precision, intent(out) :: rad_flux(ixo^s, 1:ndim)
614 
615  double precision :: grad_r_e(ixi^s), grad_r_eo(ixo^s)
616  double precision :: rad_e(ixi^s)
617  double precision :: kappa(ixo^s,1:ndim), lambda(ixo^s,1:ndim), fld_r(ixo^s,1:ndim)
618  integer :: ix^d, idir
619 
620  rad_e(ixi^s) = w(ixi^s, iw_r_e)
621 
622  call afld_get_opacity(w, x, ixi^l, ixo^l, kappa)
623  call afld_get_fluxlimiter(w, x, ixi^l, ixo^l, lambda, fld_r)
624 
625  !> Calculate the Flux using the fld closure relation
626  !> F = -c*lambda/(kappa*rho) *grad E
627  do idir = 1,ndim
628  call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo,nghostcells)
629  rad_flux(ixo^s, idir) = -(const_c/unit_velocity)*lambda(ixo^s,idir)/(kappa(ixo^s,idir)*w(ixo^s,iw_rho))*grad_r_eo(ixo^s)
630 
631  ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e)
632  ! rad_flux(ixO^S, idir) = -(const_c/unit_velocity)*lambda(ixO^S)/(kappa(ixO^S)*w(ixO^S,iw_rho))*grad_r_e(ixO^S)
633  end do
634 
635  end subroutine afld_get_radflux
636 
637  !> Calculate Eddington-tensor
638  !> Stores Eddington-tensor in w-array
639  subroutine afld_get_eddington(w, x, ixI^L, ixO^L, eddington_tensor)
641  use mod_geometry
642 
643  integer, intent(in) :: ixI^L, ixO^L
644  double precision, intent(in) :: w(ixI^S, 1:nw)
645  double precision, intent(in) :: x(ixI^S, 1:ndim)
646  double precision, intent(out) :: eddington_tensor(ixO^S,1:ndim,1:ndim)
647  double precision :: tnsr2(ixO^S,1:ndim,1:ndim)
648  double precision :: normgrad2(ixO^S), f(ixO^S,1:ndim)
649  double precision :: grad_r_e(ixI^S, 1:ndim), rad_e(ixI^S)
650  double precision :: grad_r_eO(ixO^S, 1:ndim)
651  double precision :: lambda(ixO^S,1:ndim), fld_R(ixO^S,1:ndim)
652  integer :: i,j, idir,jdir
653 
654  call afld_get_fluxlimiter(w, x, ixi^l, ixo^l, lambda, fld_r)
655 
656  !> Calculate R everywhere
657  !> |grad E|/(rho kappa E)
658  normgrad2(ixo^s) = zero
659 
660  rad_e(ixi^s) = w(ixi^s, iw_r_e)
661  do idir = 1,ndim
662  call gradiento(rad_e,x,ixi^l,ixo^l,idir,grad_r_eo(ixo^s, idir),nghostcells)
663  normgrad2(ixo^s) = normgrad2(ixo^s) + grad_r_eo(ixo^s,idir)**2
664 
665  ! call gradient(rad_e,ixI^L,ixO^L,idir,grad_r_e(ixI^S,idir))
666  ! normgrad2(ixO^S) = normgrad2(ixO^S) + grad_r_e(ixO^S,idir)**two
667 
668  !> Calculate radiation pressure
669  !> P = (lambda + lambda^2 R^2)*E
670  f(ixo^s,idir) = lambda(ixo^s,idir) + lambda(ixo^s,idir)**two * fld_r(ixo^s,idir)**two
671  f(ixo^s,idir) = one/two*(one-f(ixo^s,idir)) + one/two*(3.d0*f(ixo^s,idir) - one)
672  end do
673 
674 
675 
676  {^ifoned
677  eddington_tensor(ixo^s,1,1) = f(ixo^s,1)
678  }
679 
680  {^nooned
681  do idir = 1,ndim
682  eddington_tensor(ixo^s,idir,idir) = half*(one-f(ixo^s,idir))
683  enddo
684 
685  do idir = 1,ndim
686  do jdir = 1,ndim
687  if (idir .ne. jdir) eddington_tensor(ixo^s,idir,jdir) = zero
688  tnsr2(ixo^s,idir,jdir) = half*(3.d0*f(ixo^s,idir) - 1)&
689  *grad_r_eo(ixo^s,idir)*grad_r_eo(ixo^s,jdir)/normgrad2(ixo^s)
690  enddo
691  enddo
692 
693  ! do idir = 1,ndim
694  ! do jdir = 1,ndim
695  ! if (idir .ne. jdir) eddington_tensor(ixO^S,idir,jdir) = zero
696  ! tnsr2(ixO^S,idir,jdir) = half*(3.d0*f(ixO^S) - 1)&
697  ! *grad_r_e(ixO^S,idir)*grad_r_e(ixO^S,jdir)/normgrad2(ixO^S)
698  ! enddo
699  ! enddo
700 
701  do idir = 1,ndim
702  do jdir = 1,ndim
703  where ((tnsr2(ixo^s,idir,jdir) .eq. tnsr2(ixo^s,idir,jdir)) &
704  .and. (normgrad2(ixo^s) .gt. smalldouble))
705  eddington_tensor(ixo^s,idir,jdir) = eddington_tensor(ixo^s,idir,jdir) + tnsr2(ixo^s,idir,jdir)
706  endwhere
707  enddo
708  enddo
709  }
710 
711  end subroutine afld_get_eddington
712 
713  !> Calculate Radiation Pressure
714  !> Returns Radiation Pressure as tensor
715  subroutine afld_get_radpress(w, x, ixI^L, ixO^L, rad_pressure)
717 
718  integer, intent(in) :: ixi^l, ixo^l
719  double precision, intent(in) :: w(ixi^s, 1:nw)
720  double precision, intent(in) :: x(ixi^s, 1:ndim)
721  double precision :: eddington_tensor(ixo^s,1:ndim,1:ndim)
722  double precision, intent(out):: rad_pressure(ixo^s,1:ndim,1:ndim)
723 
724  integer i,j
725 
726  call afld_get_eddington(w, x, ixi^l, ixo^l, eddington_tensor)
727 
728  do i=1,ndim
729  do j=1,ndim
730  rad_pressure(ixo^s,i,j) = eddington_tensor(ixo^s,i,j)* w(ixo^s,iw_r_e)
731  enddo
732  enddo
733  end subroutine afld_get_radpress
734 
735  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
736  !!!!!!!!!!!!!!!!!!! Multigrid diffusion
737  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
738 
739 
740 ! !> Calling all subroutines to perform the multigrid method
741 ! !> Communicates rad_e and diff_coeff to multigrid library
742 ! subroutine Diffuse_E_rad_mg(dtfactor,qdt,qtC,psa,psb)
743 ! use mod_global_parameters
744 ! use mod_forest
745 ! use mod_ghostcells_update
746 ! use mod_multigrid_coupling
747 ! use mod_physics, only: phys_set_mg_bounds, phys_req_diagonal
748 !
749 ! type(state), target :: psa(max_blocks)
750 ! type(state), target :: psb(max_blocks)
751 ! double precision, intent(in) :: qdt
752 ! double precision, intent(in) :: qtC
753 ! double precision, intent(in) :: dtfactor
754 !
755 ! integer :: n
756 ! double precision :: res, max_residual, lambda
757 ! integer, parameter :: max_its = 50
758 !
759 ! integer :: iw_to,iw_from
760 ! integer :: iigrid, igrid, id
761 ! integer :: nc, lvl, i
762 ! type(tree_node), pointer :: pnode
763 ! real(dp) :: fac, facD
764 !
765 !
766 ! !> Set diffusion timestep, add previous timestep if mg did not converge:
767 ! dt_diff = dt_diff + qdt
768 !
769 ! ! Avoid setting a very restrictive limit to the residual when the time step
770 ! ! is small (as the operator is ~ 1/(D * qdt))
771 ! if (qdt < dtmin) then
772 ! if(mype==0)then
773 ! print *,'skipping implicit solve: dt too small!'
774 ! print *,'Currently at time=',global_time,' time step=',qdt,' dtmin=',dtmin
775 ! endif
776 ! return
777 ! endif
778 ! ! max_residual = 1d-7/qdt
779 ! max_residual = fld_diff_tol !1d-7/qdt
780 !
781 ! mg%operator_type = mg_ahelmholtz
782 ! mg%smoother_type = mg_smoother_gsrb
783 ! call mg_set_methods(mg)
784 !
785 ! if (.not. mg%is_allocated) call mpistop("multigrid tree not allocated yet")
786 !
787 ! ! lambda = 1.d0/(dtfactor * qdt)
788 ! lambda = 1.d0/(dtfactor * dt_diff)
789 ! call ahelmholtz_set_lambda(lambda)
790 !
791 ! call update_diffcoeff(psb)
792 !
793 ! fac = 1.d0
794 ! facD = 1.d0
795 !
796 ! !This is mg_copy_to_tree from psb state
797 ! {^IFONED
798 ! ! call mg_copy_to_tree(i_diff_mg(1), mg_iveps, factor=facD, state_from=psb)
799 ! !This is mg_copy_to_tree from psb state
800 ! !!! replaces:: call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
801 ! iw_from=i_diff_mg(1)
802 ! iw_to=mg_iveps
803 ! fac=facD
804 ! do iigrid=1,igridstail; igrid=igrids(iigrid);
805 ! pnode => igrid_to_node(igrid, mype)%node
806 ! id = pnode%id
807 ! lvl = mg%boxes(id)%lvl
808 ! nc = mg%box_size_lvl(lvl)
809 ! ! Include one layer of ghost cells on grid leaves
810 !
811 ! mg%boxes(id)%cc(0:nc+1, iw_to) = fac * &
812 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, iw_from)
813 !
814 ! end do
815 !
816 ! if (time_advance)then
817 ! call mg_restrict(mg, mg_iveps)
818 ! call mg_fill_ghost_cells(mg, mg_iveps)
819 ! endif
820 ! }
821 !
822 ! {^IFTWOD
823 ! ! call mg_copy_to_tree(i_diff_mg(1), mg_iveps1, factor=facD, state_from=psb)
824 ! !This is mg_copy_to_tree from psb state
825 ! !!! replaces:: call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
826 ! iw_from=i_diff_mg(1)
827 ! iw_to=mg_iveps1
828 ! fac=facD
829 ! do iigrid=1,igridstail; igrid=igrids(iigrid);
830 ! pnode => igrid_to_node(igrid, mype)%node
831 ! id = pnode%id
832 ! lvl = mg%boxes(id)%lvl
833 ! nc = mg%box_size_lvl(lvl)
834 ! ! Include one layer of ghost cells on grid leaves
835 !
836 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_to) = fac * &
837 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, iw_from)
838 !
839 ! end do
840 !
841 ! ! call mg_copy_to_tree(i_diff_mg(2), mg_iveps2, factor=facD, state_from=psb)
842 ! !This is mg_copy_to_tree from psb state
843 ! !!! replaces:: call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
844 ! iw_from=i_diff_mg(2)
845 ! iw_to=mg_iveps2
846 ! fac=facD
847 ! do iigrid=1,igridstail; igrid=igrids(iigrid);
848 ! pnode => igrid_to_node(igrid, mype)%node
849 ! id = pnode%id
850 ! lvl = mg%boxes(id)%lvl
851 ! nc = mg%box_size_lvl(lvl)
852 ! ! Include one layer of ghost cells on grid leaves
853 !
854 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_to) = fac * &
855 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, iw_from)
856 !
857 ! end do
858 !
859 ! if (time_advance)then
860 ! call mg_restrict(mg, mg_iveps1)
861 ! call mg_fill_ghost_cells(mg, mg_iveps1)
862 !
863 ! call mg_restrict(mg, mg_iveps2)
864 ! call mg_fill_ghost_cells(mg, mg_iveps2)
865 ! endif
866 ! }
867 !
868 ! {^IFTHREED
869 ! ! call mg_copy_to_tree(i_diff_mg(1), mg_iveps1, factor=facD, state_from=psb)
870 ! !This is mg_copy_to_tree from psb state
871 ! !!! replaces:: call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
872 ! iw_from=i_diff_mg(1)
873 ! iw_to=mg_iveps1
874 ! fac=facD
875 ! do iigrid=1,igridstail; igrid=igrids(iigrid);
876 ! pnode => igrid_to_node(igrid, mype)%node
877 ! id = pnode%id
878 ! lvl = mg%boxes(id)%lvl
879 ! nc = mg%box_size_lvl(lvl)
880 ! ! Include one layer of ghost cells on grid leaves
881 !
882 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_to) = fac * &
883 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
884 ! ixMlo3-1:ixMhi3+1, iw_from)
885 !
886 ! end do
887 !
888 ! ! call mg_copy_to_tree(i_diff_mg(2), mg_iveps2, factor=facD, state_from=psb)
889 ! !This is mg_copy_to_tree from psb state
890 ! !!! replaces:: call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
891 ! iw_from=i_diff_mg(2)
892 ! iw_to=mg_iveps2
893 ! fac=facD
894 ! do iigrid=1,igridstail; igrid=igrids(iigrid);
895 ! pnode => igrid_to_node(igrid, mype)%node
896 ! id = pnode%id
897 ! lvl = mg%boxes(id)%lvl
898 ! nc = mg%box_size_lvl(lvl)
899 ! ! Include one layer of ghost cells on grid leaves
900 !
901 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_to) = fac * &
902 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
903 ! ixMlo3-1:ixMhi3+1, iw_from)
904 !
905 ! end do
906 !
907 ! ! call mg_copy_to_tree(i_diff_mg(3), mg_iveps2, factor=facD, state_from=psb)
908 ! !This is mg_copy_to_tree from psb state
909 ! !!! replaces:: call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
910 ! iw_from=i_diff_mg(3)
911 ! iw_to=mg_iveps3
912 ! fac=facD
913 ! do iigrid=1,igridstail; igrid=igrids(iigrid);
914 ! pnode => igrid_to_node(igrid, mype)%node
915 ! id = pnode%id
916 ! lvl = mg%boxes(id)%lvl
917 ! nc = mg%box_size_lvl(lvl)
918 ! ! Include one layer of ghost cells on grid leaves
919 !
920 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_to) = fac * &
921 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
922 ! ixMlo3-1:ixMhi3+1, iw_from)
923 !
924 ! end do
925 !
926 ! if (time_advance)then
927 ! call mg_restrict(mg, mg_iveps1)
928 ! call mg_fill_ghost_cells(mg, mg_iveps1)
929 !
930 ! call mg_restrict(mg, mg_iveps2)
931 ! call mg_fill_ghost_cells(mg, mg_iveps2)
932 !
933 ! call mg_restrict(mg, mg_iveps3)
934 ! call mg_fill_ghost_cells(mg, mg_iveps3)
935 ! endif
936 ! }
937 !
938 !
939 ! !This is mg_copy_to_tree from psb state
940 ! ! call mg_copy_to_tree(iw_r_e, mg_iphi, factor=fac, state_from=psb)
941 ! !This is mg_copy_to_tree from psb state
942 ! !!! replaces:: call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
943 ! iw_from=iw_r_e
944 ! iw_to=mg_iphi
945 ! fac=-lambda
946 ! do iigrid=1,igridstail; igrid=igrids(iigrid);
947 ! pnode => igrid_to_node(igrid, mype)%node
948 ! id = pnode%id
949 ! lvl = mg%boxes(id)%lvl
950 ! nc = mg%box_size_lvl(lvl)
951 ! ! Include one layer of ghost cells on grid leaves
952 ! {^IFONED
953 ! mg%boxes(id)%cc(0:nc+1, iw_to) = fac * &
954 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, iw_from)
955 ! }
956 ! {^IFTWOD
957 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_to) = fac * &
958 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, iw_from)
959 ! }
960 ! {^IFTHREED
961 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_to) = fac * &
962 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
963 ! ixMlo3-1:ixMhi3+1, iw_from)
964 ! }
965 ! end do
966 !
967 !
968 ! !>replace call set_rhs(mg, -1/dt, 0.0_dp)
969 ! ! call mg_copy_to_tree(iw_r_e, mg_irhs, factor=-1/(dtfactor*dt_diff), state_from=psb)
970 ! !This is mg_copy_to_tree from psb state
971 ! !!! replaces:: call mg_copy_to_tree(su_, mg_irhs, factor=-lambda)
972 ! iw_from=iw_r_e
973 ! iw_to=mg_irhs
974 ! fac=-1/(dtfactor*dt_diff)
975 ! do iigrid=1,igridstail; igrid=igrids(iigrid);
976 ! pnode => igrid_to_node(igrid, mype)%node
977 ! id = pnode%id
978 ! lvl = mg%boxes(id)%lvl
979 ! nc = mg%box_size_lvl(lvl)
980 ! ! Include one layer of ghost cells on grid leaves
981 ! {^IFONED
982 ! mg%boxes(id)%cc(0:nc+1, iw_to) = fac * &
983 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, iw_from)
984 ! }
985 ! {^IFTWOD
986 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_to) = fac * &
987 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, iw_from)
988 ! }
989 ! {^IFTHREED
990 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_to) = fac * &
991 ! psb(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
992 ! ixMlo3-1:ixMhi3+1, iw_from)
993 ! }
994 ! end do
995 !
996 !
997 ! call phys_set_mg_bounds()
998 !
999 ! call mg_fas_fmg(mg, .true., max_res=res)
1000 ! do n = 1, max_its
1001 ! ! print*, n, res
1002 ! if (res < max_residual) exit
1003 ! call mg_fas_vcycle(mg, max_res=res)
1004 ! end do
1005 !
1006 ! if (res .le. 0.d0) then
1007 ! if (diffcrash_resume) then
1008 ! if (mg%my_rank == 0) &
1009 ! write(*,*) it, ' resiudal zero ', res
1010 ! goto 0888
1011 ! endif
1012 ! if (mg%my_rank == 0) then
1013 ! print*, res
1014 ! error stop "Diffusion residual to zero"
1015 ! endif
1016 ! endif
1017 !
1018 ! if (n == max_its + 1) then
1019 ! if (diffcrash_resume) then
1020 ! if (mg%my_rank == 0) &
1021 ! write(*,*) it, ' resiudal high ', res
1022 ! if (res .lt. 1.d3*max_residual) goto 0887
1023 ! goto 0888
1024 ! endif
1025 ! if (mg%my_rank == 0) then
1026 ! print *, "Did you specify boundary conditions correctly?"
1027 ! print *, "Or is the variation in diffusion too large?"
1028 ! print*, n, res
1029 ! print *, mg%bc(1, mg_iphi)%bc_value, mg%bc(2, mg_iphi)%bc_value
1030 ! end if
1031 ! error stop "diffusion_solve: no convergence"
1032 ! end if
1033 !
1034 ! !> Reset dt_diff when diffusion worked out
1035 ! 0887 dt_diff = 0.d0
1036 !
1037 ! ! !This is mg_copy_from_tree_gc for psa state
1038 ! ! call mg_copy_from_tree_gc(mg_iphi, iw_r_e, state_to=psa)
1039 ! !This is mg_copy_from_tree_gc for psa state
1040 ! !!! replaces:: call mg_copy_from_tree_gc(mg_iphi, su_)
1041 ! iw_from=mg_iphi
1042 ! iw_to=iw_r_e
1043 ! do iigrid=1,igridstail; igrid=igrids(iigrid);
1044 ! pnode => igrid_to_node(igrid, mype)%node
1045 ! id = pnode%id
1046 ! lvl = mg%boxes(id)%lvl
1047 ! nc = mg%box_size_lvl(lvl)
1048 ! ! Include one layer of ghost cells on grid leaves
1049 ! {^IFONED
1050 ! psa(igrid)%w(ixMlo1-1:ixMhi1+1, iw_to) = &
1051 ! mg%boxes(id)%cc(0:nc+1, iw_from)
1052 ! }
1053 ! {^IFTWOD
1054 ! psa(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, iw_to) = &
1055 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, iw_from)
1056 ! }
1057 ! {^IFTHREED
1058 ! psa(igrid)%w(ixMlo1-1:ixMhi1+1, ixMlo2-1:ixMhi2+1, &
1059 ! ixMlo3-1:ixMhi3+1, iw_to) = &
1060 ! mg%boxes(id)%cc(0:nc+1, 0:nc+1, 0:nc+1, iw_from)
1061 ! }
1062 ! end do
1063 !
1064 ! ! enforce boundary conditions for psa
1065 ! 0888 call getbc(qtC,0.d0,psa,1,nwflux+nwaux,phys_req_diagonal)
1066 !
1067 ! end subroutine Diffuse_E_rad_mg
1068 
1069 !> Calling all subroutines to perform the multigrid method
1070 !> Communicates rad_e and diff_coeff to multigrid library
1071 subroutine diffuse_e_rad_mg(dtfactor,qdt,qtC,psa,psb)
1073  use mod_forest
1077 
1078  type(state), target :: psa(max_blocks)
1079  type(state), target :: psb(max_blocks)
1080  double precision, intent(in) :: qdt
1081  double precision, intent(in) :: qtC
1082  double precision, intent(in) :: dtfactor
1083 
1084  integer :: n
1085  double precision :: res, max_residual, lambda
1086  integer, parameter :: max_its = 50
1087 
1088  integer :: iw_to,iw_from
1089  integer :: iigrid, igrid, id
1090  integer :: nc, lvl, i
1091  type(tree_node), pointer :: pnode
1092  real(dp) :: fac, facD
1093 
1094 
1095  !> Set diffusion timestep, add previous timestep if mg did not converge:
1096  dt_diff = dt_diff + qdt
1097 
1098  ! Avoid setting a very restrictive limit to the residual when the time step
1099  ! is small (as the operator is ~ 1/(D * qdt))
1100  if (qdt < dtmin) then
1101  if(mype==0)then
1102  print *,'skipping implicit solve: dt too small!'
1103  print *,'Currently at time=',global_time,' time step=',qdt,' dtmin=',dtmin
1104  endif
1105  return
1106  endif
1107  ! max_residual = 1d-7/qdt
1108  max_residual = fld_diff_tol !1d-7/qdt
1109 
1110  mg%operator_type = mg_ahelmholtz
1111  mg%smoother_type = mg_smoother_gsrb
1112  call mg_set_methods(mg)
1113 
1114  if (.not. mg%is_allocated) call mpistop("multigrid tree not allocated yet")
1115 
1116 ! lambda = 1.d0/(dtfactor * qdt)
1117  lambda = 1.d0/(dtfactor * dt_diff)
1118  call ahelmholtz_set_lambda(lambda)
1119 
1120  call update_diffcoeff(psb)
1121 
1122  fac = 1.d0
1123  facd = 1.d0
1124 
1125  !This is mg_copy_to_tree from psb state
1126  {^ifoned
1127  call mg_copy_to_tree(i_diff_mg(1), mg_iveps, factor=facd, state_from=psb)
1128 
1129  if (time_advance)then
1130  call mg_restrict(mg, mg_iveps)
1131  call mg_fill_ghost_cells(mg, mg_iveps)
1132  endif
1133  }
1134 
1135  {^iftwod
1136  call mg_copy_to_tree(i_diff_mg(1), mg_iveps1, factor=facd, state_from=psb)
1137  call mg_copy_to_tree(i_diff_mg(2), mg_iveps2, factor=facd, state_from=psb)
1138 
1139  if (time_advance)then
1140  call mg_restrict(mg, mg_iveps1)
1141  call mg_fill_ghost_cells(mg, mg_iveps1)
1142 
1143  call mg_restrict(mg, mg_iveps2)
1144  call mg_fill_ghost_cells(mg, mg_iveps2)
1145  endif
1146  }
1147 
1148 
1149  !This is mg_copy_to_tree from psb state
1150  call mg_copy_to_tree(iw_r_e, mg_iphi, factor=fac, state_from=psb)
1151 
1152  !>replace call set_rhs(mg, -1/dt, 0.0_dp)
1153 ! call mg_copy_to_tree(iw_r_e, mg_irhs, factor=-1/(dtfactor*qdt), state_from=psb)
1154  call mg_copy_to_tree(iw_r_e, mg_irhs, factor=-1/(dtfactor*dt_diff), state_from=psb)
1155 
1156  call phys_set_mg_bounds()
1157 
1158  call mg_fas_fmg(mg, .true., max_res=res)
1159  do n = 1, max_its
1160  ! print*, n, res
1161  if (res < max_residual) exit
1162  call mg_fas_vcycle(mg, max_res=res)
1163  end do
1164 
1165  if (res .le. 0.d0) then
1166  if (diffcrash_resume) then
1167  if (mg%my_rank == 0) &
1168  write(*,*) it, ' resiudal zero ', res
1169  return
1170  endif
1171  if (mg%my_rank == 0) then
1172  print*, res
1173  error stop "Diffusion residual to zero"
1174  endif
1175  endif
1176 
1177  if (n == max_its + 1) then
1178  if (diffcrash_resume) then
1179  if (mg%my_rank == 0) &
1180  write(*,*) it, ' resiudal high ', res
1181  return
1182  endif
1183  if (mg%my_rank == 0) then
1184  print *, "Did you specify boundary conditions correctly?"
1185  print *, "Or is the variation in diffusion too large?"
1186  print*, n, res
1187  print *, mg%bc(1, mg_iphi)%bc_value, mg%bc(2, mg_iphi)%bc_value
1188  end if
1189  error stop "diffusion_solve: no convergence"
1190  end if
1191 
1192  !> Reset dt_diff when diffusion worked out
1193  dt_diff = 0.d0
1194 
1195  ! !This is mg_copy_from_tree_gc for psa state
1196  call mg_copy_from_tree_gc(mg_iphi, iw_r_e, state_to=psa)
1197 !
1198 ! ! enforce boundary conditions for psa
1199 ! 0888 call getbc(qtC,0.d0,psa,1,nwflux+nwaux,phys_req_diagonal)
1200 
1201 end subroutine diffuse_e_rad_mg
1202 
1203 
1204  !> inplace update of psa==>F_im(psa)
1205  subroutine evaluate_e_rad_mg(qtC,psa)
1208  use mod_physics, only: phys_req_diagonal
1209 
1210  type(state), target :: psa(max_blocks)
1211  double precision, intent(in) :: qtC
1212 
1213  integer :: iigrid, igrid, level
1214  integer :: ixO^L
1215 
1216  call update_diffcoeff(psa)
1217 
1218  ixo^l=ixg^ll^lsub1;
1219  !$OMP PARALLEL DO PRIVATE(igrid)
1220  do iigrid=1,igridstail; igrid=igrids(iigrid);
1221  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
1222  ! call afld_get_diffcoef_central(psa(igrid)%w, psa(igrid)%w, psa(igrid)%x, ixG^LL, ixO^L)
1223  call put_diffterm_onegrid(ixg^ll,ixo^l,psa(igrid)%w)
1224  end do
1225  !$OMP END PARALLEL DO
1226 
1227  ! enforce boundary conditions for psa
1228  call getbc(qtc,0.d0,psa,1,nwflux+nwaux,phys_req_diagonal)
1229 
1230  end subroutine evaluate_e_rad_mg
1231 
1232  !> inplace update of psa==>F_im(psa)
1233  subroutine put_diffterm_onegrid(ixI^L,ixO^L,w)
1235  integer, intent(in) :: ixI^L, ixO^L
1236  double precision, intent(inout) :: w(ixI^S, 1:nw)
1237 
1238  double precision :: gradE(ixO^S),divF(ixO^S)
1239  double precision :: divF_h(ixO^S),divF_j(ixO^S)
1240  double precision :: diff_term(ixO^S)
1241 
1242  integer :: idir, jxO^L, hxO^L
1243 
1244  ! call mpistop("phys_evaluate_implicit not implemented for FLD")
1245 
1246  divf(ixo^s) = 0.d0
1247  do idir = 1,ndim
1248  hxo^l=ixo^l-kr(idir,^d);
1249  jxo^l=ixo^l+kr(idir,^d);
1250 
1251  divf_h(ixo^s) = w(ixo^s,i_diff_mg(idir))*w(hxo^s,i_diff_mg(idir))/(w(ixo^s,i_diff_mg(idir)) + w(hxo^s,i_diff_mg(idir)))*(w(hxo^s,iw_r_e) - w(ixo^s,iw_r_e))
1252  divf_j(ixo^s) = w(ixo^s,i_diff_mg(idir))*w(jxo^s,i_diff_mg(idir))/(w(ixo^s,i_diff_mg(idir)) + w(jxo^s,i_diff_mg(idir)))*(w(jxo^s,iw_r_e) - w(ixo^s,iw_r_e))
1253  divf(ixo^s) = divf(ixo^s) + 2.d0*(divf_h(ixo^s) + divf_j(ixo^s))/dxlevel(idir)**2
1254  enddo
1255 
1256  w(ixo^s,iw_r_e) = divf(ixo^s)
1257 
1258  end subroutine put_diffterm_onegrid
1259 
1260 
1261  !> Calculates cell-centered diffusion coefficient to be used in multigrid
1262  subroutine afld_get_diffcoef_central(w, wCT, x, ixI^L, ixO^L)
1264  use mod_geometry
1265  use mod_usr_methods
1266 
1267  integer, intent(in) :: ixI^L, ixO^L
1268  double precision, intent(inout) :: w(ixI^S, 1:nw)
1269  double precision, intent(in) :: wCT(ixI^S, 1:nw)
1270  double precision, intent(in) :: x(ixI^S, 1:ndim)
1271 
1272  double precision :: kappa(ixO^S,1:ndim), lambda(ixO^S,1:ndim), fld_R(ixO^S,1:ndim)
1273 
1274  double precision :: max_D(ixI^S), grad_r_e(ixI^S), rad_e(ixI^S)
1275  integer :: idir,i,j, ix^D, idim
1276 
1277  if (fld_diff_testcase) then
1278 
1279  w(ixo^s,i_diff_mg(:)) = 1.d0
1280 
1281  else
1282 
1283  call afld_get_opacity(wct, x, ixi^l, ixo^l, kappa)
1284  call afld_get_fluxlimiter(wct, x, ixi^l, ixo^l, lambda, fld_r)
1285 
1286  do idim = 1,ndim
1287  !> calculate diffusion coefficient
1288  w(ixo^s,i_diff_mg(idim)) = (const_c/unit_velocity)*lambda(ixo^s,idim)/(kappa(ixo^s,idim)*wct(ixo^s,iw_rho))
1289 
1290  where (w(ixo^s,i_diff_mg(idim)) .lt. 0.d0)
1291  w(ixo^s,i_diff_mg(idim)) = smalldouble
1292  end where
1293 
1294  if (diff_coef_filter) then
1295  !call mpistop('Hold your bloody horses, not implemented yet ')
1296  call afld_smooth_diffcoef(w, ixi^l, ixo^l,idim)
1297  endif
1298  enddo
1299  endif
1300 
1301  if (associated(usr_special_diffcoef)) &
1302  call usr_special_diffcoef(w, wct, x, ixi^l, ixo^l)
1303 
1304  end subroutine afld_get_diffcoef_central
1305 
1306  subroutine update_diffcoeff(psa)
1308 
1309  type(state), target :: psa(max_blocks)
1310 
1311  ! double precision :: wCT(ixG^LL,1:nw)
1312  integer :: iigrid, igrid, level
1313  integer :: ixO^L
1314 
1315  ixo^l=ixg^ll^lsub1;
1316  !$OMP PARALLEL DO PRIVATE(igrid)
1317  do iigrid=1,igridstail; igrid=igrids(iigrid);
1318  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
1319 
1320  ! wCT = psa(igrid)%w
1321  call afld_get_diffcoef_central(psa(igrid)%w, psa(igrid)%w, psa(igrid)%x, ixg^ll, ixo^l)
1322  end do
1323  !$OMP END PARALLEL DO
1324 
1325  end subroutine update_diffcoeff
1326 
1327  !> Use running average on Diffusion coefficient
1328  subroutine afld_smooth_diffcoef(w, ixI^L, ixO^L,idim)
1330 
1331  integer, intent(in) :: ixI^L, ixO^L, idim
1332  double precision, intent(inout) :: w(ixI^S, 1:nw)
1333 
1334  double precision :: tmp_D(ixI^S), filtered_D(ixI^S)
1335  integer :: ix^D, filter, idir
1336 
1337  if (size_d_filter .lt. 1) call mpistop("D filter of size < 1 makes no sense")
1338  if (size_d_filter .gt. nghostcells) call mpistop("D filter of size > nghostcells makes no sense")
1339 
1340  tmp_d(ixo^s) = w(ixo^s,i_diff_mg(idim))
1341  filtered_d(ixo^s) = zero
1342 
1343  do filter = 1,size_d_filter
1344  {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
1345  do idir = 1,ndim
1346  filtered_d(ix^d) = filtered_d(ix^d) &
1347  + tmp_d(ix^d+filter*kr(idir,^d)) &
1348  + tmp_d(ix^d-filter*kr(idir,^d))
1349  enddo
1350  {enddo\}
1351  enddo
1352 
1353  {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
1354  tmp_d(ix^d) = (tmp_d(ix^d)+filtered_d(ix^d))/(1+2*size_d_filter*ndim)
1355  {enddo\}
1356 
1357  w(ixo^s,i_diff_mg(idim)) = tmp_d(ixo^s)
1358  end subroutine afld_smooth_diffcoef
1359 
1360 
1361 
1362  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1363  !!!!!!!!!!!!!!!!!!! Gas-Rad Energy interaction
1364  !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
1365 
1366  !> This subroutine calculates the radiation heating, radiation cooling
1367  !> and photon tiring using an implicit scheme.
1368  !> These sourceterms are applied using the root-finding of a 4th order polynomial
1369  !> This routine loops over every cell in the domain
1370  !> and computes the coefficients of the polynomials in every cell
1371  subroutine energy_interaction(w, wCT, x, ixI^L, ixO^L)
1373  use mod_geometry
1374  use mod_physics
1375  use mod_usr_methods
1376 
1377  integer, intent(in) :: ixI^L, ixO^L
1378  double precision, intent(in) :: x(ixI^S,1:ndim)
1379  double precision, intent(in) :: wCT(ixI^S,1:nw)
1380  double precision, intent(inout) :: w(ixI^S,1:nw)
1381 
1382  double precision :: a1(ixO^S), a2(ixO^S)
1383  double precision :: c0(ixO^S), c1(ixO^S)
1384  double precision :: e_gas(ixO^S), E_rad(ixO^S)
1385  double precision :: kappa(ixO^S)
1386  double precision :: sigma_b
1387 
1388  integer :: i,j,ix^D
1389 
1390  ! if (fld_interaction_method .eq. 'Instant') then
1391  ! call Instant_qdot(w, w, ixI^L, ixO^L)
1392  ! return
1393  ! endif
1394 
1395  !> e_gas is the INTERNAL ENERGY without KINETIC ENERGY
1396  ! if (.not. block%e_is_internal) then
1397  e_gas(ixo^s) = wct(ixo^s,iw_e) - half*sum(wct(ixo^s, iw_mom(:))**2, dim=ndim+1)/wct(ixo^s, iw_rho)
1398  ! else
1399  ! e_gas(ixO^S) = wCT(ixO^S,iw_e)
1400  ! endif
1401 
1402  {do ix^d=ixomin^d,ixomax^d\ }
1403  e_gas(ix^d) = max(e_gas(ix^d),small_pressure/(afld_gamma-1.d0))
1404  {enddo\}
1405 
1406  e_rad(ixo^s) = wct(ixo^s,iw_r_e)
1407 
1408  if (associated(usr_special_opacity_qdot)) then
1409  call usr_special_opacity_qdot(ixi^l,ixo^l,w,x,kappa)
1410  else
1411  call mpistop("apoint qdot opacity")
1412  call afld_get_opacity(wct, x, ixi^l, ixo^l, kappa)
1413  endif
1414 
1415  sigma_b = const_rad_a*const_c/4.d0*(unit_temperature**4.d0)/(unit_velocity*unit_pressure)
1416 
1417  if (fld_interaction_method .eq. 'Instant') then
1418  a1(ixo^s) = const_rad_a*(afld_mu*const_mp/const_kb*(afld_gamma-1))**4/(wct(ixo^s,iw_rho)*unit_density)**4 &
1419  /unit_pressure**3
1420  a2(ixo^s) = e_gas(ixo^s) + e_rad(ixo^s)
1421 
1422  c0(ixo^s) = a2(ixo^s)/a1(ixo^s)
1423  c1(ixo^s) = 1.d0/a1(ixo^s)
1424  else
1425  !> Calculate coefficients for polynomial
1426  a1(ixo^s) = 4*kappa(ixo^s)*sigma_b*(afld_gamma-one)**4/wct(ixo^s,iw_rho)**3*dt
1427  a2(ixo^s) = (const_c/unit_velocity)*kappa(ixo^s)*wct(ixo^s,iw_rho)*dt
1428 
1429  c0(ixo^s) = ((one + a2(ixo^s))*e_gas(ixo^s) + a2(ixo^s)*e_rad(ixo^s))/a1(ixo^s)
1430  c1(ixo^s) = (one + a2(ixo^s))/a1(ixo^s)
1431  endif
1432 
1433  ! w(ixO^S,i_test) = a2(ixO^S)
1434 
1435  !> Loop over every cell for rootfinding method
1436  {do ix^d=ixomin^d,ixomax^d\ }
1437  select case(fld_interaction_method)
1438  case('Bisect')
1439  call bisection_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
1440  case('Newton')
1441  call newton_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
1442  case('Halley')
1443  call halley_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
1444  case('Instant')
1445  call halley_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
1446  case default
1447  call mpistop('root-method not known')
1448  end select
1449  {enddo\}
1450 
1451  if (fld_interaction_method .eq. 'Instant') then
1452  e_rad(ixo^s) = a2(ixo^s) - e_gas(ixo^s)
1453  else
1454  !> advance E_rad
1455  e_rad(ixo^s) = (a1(ixo^s)*e_gas(ixo^s)**4.d0 + e_rad(ixo^s))/(one + a2(ixo^s))
1456  endif
1457 
1458  !> new w = w + dt f(wCT)
1459  !> e_gas,E_rad = wCT + dt f(wCT)
1460  !> dt f(wCT) = e_gas,E_rad - wCT
1461  !> new w = w + e_gas,E_rad - wCT
1462 
1463  !> WAIT A SECOND?! DOES THIS EVEN WORK WITH THESE KIND OF IMPLICIT METHODS?
1464  !> NOT QUITE SURE TO BE HONEST
1465  !> IS IT POSSIBLE TO SHUT DOWN SOURCE SPLITTING FOR JUST THIS TERM?
1466  !> FIX BY PERFORMING Energy_interaction on (w,w,...)
1467 
1468  ! call mpistop('This still has to be fixed somehow')
1469 
1470  !> Update gas-energy in w, internal + kinetic
1471  ! w(ixO^S,iw_e) = w(ixO^S,iw_e) + e_gas(ixO^S) - wCT(ixO^S,iw_e)
1472  ! {do ix^D=ixOmin^D,ixOmax^D\ }
1473  ! e_gas(ix^D) = max(e_gas(ix^D),small_pressure/(afld_gamma-1.d0))
1474  ! {enddo\}
1475 
1476  w(ixo^s,iw_e) = e_gas(ixo^s)
1477 
1478  !> Beginning of module substracted wCT Ekin
1479  !> So now add wCT Ekin
1480  ! if (.not. block%e_is_internal) then
1481  w(ixo^s,iw_e) = w(ixo^s,iw_e) + half*sum(wct(ixo^s, iw_mom(:))**2, dim=ndim+1)/wct(ixo^s, iw_rho)
1482  ! else
1483  ! w(ixO^S,iw_e) = w(ixO^S,iw_e)
1484  ! endif
1485 
1486  !> Update rad-energy in w
1487  ! w(ixO^S,iw_r_e) = w(ixO^S,iw_r_e) + E_rad(ixO^S) - wCT(ixO^S,iw_r_e)
1488  w(ixo^s,iw_r_e) = e_rad(ixo^s)
1489 
1490  end subroutine energy_interaction
1491 
1492 
1493  !> Find the root of the 4th degree polynomial using the bisection method
1494  subroutine bisection_method(e_gas, E_rad, c0, c1)
1496 
1497  double precision, intent(in) :: c0, c1
1498  double precision, intent(in) :: E_rad
1499  double precision, intent(inout) :: e_gas
1500 
1501  double precision :: bisect_a, bisect_b, bisect_c
1502  integer :: n, max_its
1503 
1504  n = 0
1505  max_its = 10000000
1506 
1507  bisect_a = zero
1508  bisect_b = 1.2d0*max(abs(c0/c1),abs(c0)**(1.d0/4.d0))
1509 
1510  ! do while (abs(Polynomial_Bisection(bisect_b, c0, c1)-Polynomial_Bisection(bisect_a, c0, c1))&
1511  ! .ge. fld_bisect_tol*min(e_gas,E_rad))
1512  do while (abs(bisect_b-bisect_a) .ge. fld_bisect_tol*min(e_gas,e_rad))
1513  bisect_c = (bisect_a + bisect_b)/two
1514 
1515  n = n +1
1516  if (n .gt. max_its) then
1517  goto 2435
1518  call mpistop('No convergece in bisection scheme')
1519  endif
1520 
1521  if (polynomial_bisection(bisect_a, c0, c1)*&
1522  polynomial_bisection(bisect_b, c0, c1) .lt. zero) then
1523 
1524  if (polynomial_bisection(bisect_a, c0, c1)*&
1525  polynomial_bisection(bisect_c, c0, c1) .lt. zero) then
1526  bisect_b = bisect_c
1527  elseif (polynomial_bisection(bisect_b, c0, c1)*&
1528  polynomial_bisection(bisect_c, c0, c1) .lt. zero) then
1529  bisect_a = bisect_c
1530  elseif (polynomial_bisection(bisect_a, c0, c1) .eq. zero) then
1531  bisect_b = bisect_a
1532  bisect_c = bisect_a
1533  goto 2435
1534  elseif (polynomial_bisection(bisect_b, c0, c1) .eq. zero) then
1535  bisect_a = bisect_b
1536  bisect_c = bisect_b
1537  goto 2435
1538  elseif (polynomial_bisection(bisect_c, c0, c1) .eq. zero) then
1539  bisect_a = bisect_c
1540  bisect_b = bisect_c
1541  goto 2435
1542  else
1543  call mpistop("Problem with fld bisection method")
1544  endif
1545  elseif (polynomial_bisection(bisect_a, c0, c1) &
1546  - polynomial_bisection(bisect_b, c0, c1) .lt. fld_bisect_tol*polynomial_bisection(bisect_a, c0, c1)) then
1547  goto 2435
1548  else
1549  bisect_a = e_gas
1550  bisect_b = e_gas
1551  print*, "IGNORING GAS-RAD ENERGY EXCHANGE ", c0, c1
1552 
1553  print*, polynomial_bisection(bisect_a, c0, c1), polynomial_bisection(bisect_b, c0, c1)
1554 
1555  if (polynomial_bisection(bisect_a, c0, c1) .le. smalldouble) then
1556  bisect_b = bisect_a
1557  elseif (polynomial_bisection(bisect_a, c0, c1) .le. smalldouble) then
1558  bisect_a = bisect_b
1559  endif
1560 
1561  goto 2435
1562 
1563  endif
1564  enddo
1565 
1566  2435 e_gas = (bisect_a + bisect_b)/two
1567  end subroutine bisection_method
1568 
1569  !> Find the root of the 4th degree polynomial using the Newton-Ralphson method
1570  subroutine newton_method(e_gas, E_rad, c0, c1)
1572 
1573  double precision, intent(in) :: c0, c1
1574  double precision, intent(in) :: E_rad
1575  double precision, intent(inout) :: e_gas
1576 
1577  double precision :: xval, yval, der, deltax
1578 
1579  integer :: ii
1580 
1581  yval = bigdouble
1582  xval = e_gas
1583  der = one
1584  deltax = one
1585 
1586  ii = 0
1587  !> Compare error with dx = dx/dy dy
1588  do while (abs(deltax) .gt. fld_bisect_tol)
1589  yval = polynomial_bisection(xval, c0, c1)
1590  der = dpolynomial_bisection(xval, c0, c1)
1591  deltax = -yval/der
1592  xval = xval + deltax
1593  ii = ii + 1
1594  if (ii .gt. 1d3) then
1595  call bisection_method(e_gas, e_rad, c0, c1)
1596  return
1597  endif
1598  enddo
1599 
1600  e_gas = xval
1601  end subroutine newton_method
1602 
1603  !> Find the root of the 4th degree polynomial using the Halley method
1604  subroutine halley_method(e_gas, E_rad, c0, c1)
1606 
1607  double precision, intent(in) :: c0, c1
1608  double precision, intent(in) :: E_rad
1609  double precision, intent(inout) :: e_gas
1610 
1611  double precision :: xval, yval, der, dder, deltax
1612 
1613  integer :: ii
1614 
1615  yval = bigdouble
1616  xval = e_gas
1617  der = one
1618  dder = one
1619  deltax = one
1620 
1621  ii = 0
1622  !> Compare error with dx = dx/dy dy
1623  do while (abs(deltax) .gt. fld_bisect_tol)
1624  yval = polynomial_bisection(xval, c0, c1)
1625  der = dpolynomial_bisection(xval, c0, c1)
1626  dder = ddpolynomial_bisection(xval, c0, c1)
1627  deltax = -two*yval*der/(two*der**2 - yval*dder)
1628  xval = xval + deltax
1629  ii = ii + 1
1630  if (ii .gt. 1d3) then
1631  ! call mpistop('Halley did not convergggge')
1632  call newton_method(e_gas, e_rad, c0, c1)
1633  return
1634  endif
1635  enddo
1636 
1637  e_gas = xval
1638  end subroutine halley_method
1639 
1640  !> Evaluate polynomial at argument e_gas
1641  function polynomial_bisection(e_gas, c0, c1) result(val)
1643 
1644  double precision, intent(in) :: e_gas
1645  double precision, intent(in) :: c0, c1
1646  double precision :: val
1647 
1648  val = e_gas**4.d0 + c1*e_gas - c0
1649  end function polynomial_bisection
1650 
1651  !> Evaluate first derivative of polynomial at argument e_gas
1652  function dpolynomial_bisection(e_gas, c0, c1) result(der)
1654 
1655  double precision, intent(in) :: e_gas
1656  double precision, intent(in) :: c0, c1
1657  double precision :: der
1658 
1659  der = 4.d0*e_gas**3.d0 + c1
1660  end function dpolynomial_bisection
1661 
1662  !> Evaluate second derivative of polynomial at argument e_gas
1663  function ddpolynomial_bisection(e_gas, c0, c1) result(dder)
1665 
1666  double precision, intent(in) :: e_gas
1667  double precision, intent(in) :: c0, c1
1668  double precision :: dder
1669 
1670  dder = 4.d0*3.d0*e_gas**2.d0
1671  end function ddpolynomial_bisection
1672 
1673  !> Calculate gradient of a scalar q within ixL in direction idir
1674  !> difference with gradient is gradq(ixO^S), NOT gradq(ixI^S)
1675  subroutine gradiento(q,x,ixI^L,ixO^L,idir,gradq,n)
1677 
1678  integer, intent(in) :: ixI^L, ixO^L, idir
1679  integer, intent(in) :: n
1680 
1681  double precision, intent(in) :: q(ixI^S), x(ixI^S,1:ndim)
1682  double precision, intent(out) :: gradq(ixO^S)
1683  integer :: jxO^L, hxO^L
1684 
1685  ! hxO^L=ixO^L-n*kr(idir,^D);
1686  ! jxO^L=ixO^L+n*kr(idir,^D);
1687  !
1688  ! if (n .gt. nghostcells) call mpistop("gradientO stencil too wide")
1689  !
1690  ! gradq(ixO^S)=(q(jxO^S)-q(hxO^S))/(2*n*dxlevel(idir))
1691  ! gradq(ixO^S)=(q(jxO^S)-q(hxO^S))/(x(jxO^S,idir)-x(hxO^S,idir))
1692 
1693  !> Using higher order derivatives with wider stencil according to:
1694  !> https://en.wikipedia.org/wiki/Finite_difference_coefficient
1695 
1696  if (n .gt. nghostcells) then
1697  call mpistop("gradientO stencil too wide")
1698  elseif (n .eq. 1) then
1699  hxo^l=ixo^l-kr(idir,^d);
1700  jxo^l=ixo^l+kr(idir,^d);
1701  gradq(ixo^s)=(q(jxo^s)-q(hxo^s))/(x(jxo^s,idir)-x(hxo^s,idir))
1702  elseif (n .eq. 2) then
1703  gradq(ixo^s) = 0.d0
1704  !> coef 2/3
1705  hxo^l=ixo^l-kr(idir,^d);
1706  jxo^l=ixo^l+kr(idir,^d);
1707  gradq(ixo^s) = gradq(ixo^s) + 2.d0/3.d0*(q(jxo^s)-q(hxo^s))
1708  !> coef -1/12
1709  hxo^l=ixo^l-2*kr(idir,^d);
1710  jxo^l=ixo^l+2*kr(idir,^d);
1711  gradq(ixo^s) = gradq(ixo^s) - 1.d0/12.d0*(q(jxo^s)-q(hxo^s))
1712  !> divide by dx
1713  gradq(ixo^s) = gradq(ixo^s)/dxlevel(idir)
1714  elseif (n .eq. 3) then
1715  gradq(ixo^s) = 0.d0
1716  !> coef 3/4
1717  hxo^l=ixo^l-kr(idir,^d);
1718  jxo^l=ixo^l+kr(idir,^d);
1719  gradq(ixo^s) = gradq(ixo^s) + 3.d0/4.d0*(q(jxo^s)-q(hxo^s))
1720  !> coef -3/20
1721  hxo^l=ixo^l-2*kr(idir,^d);
1722  jxo^l=ixo^l+2*kr(idir,^d);
1723  gradq(ixo^s) = gradq(ixo^s) - 3.d0/20.d0*(q(jxo^s)-q(hxo^s))
1724  !> coef 1/60
1725  hxo^l=ixo^l-3*kr(idir,^d);
1726  jxo^l=ixo^l+3*kr(idir,^d);
1727  gradq(ixo^s) = gradq(ixo^s) + 1.d0/60.d0*(q(jxo^s)-q(hxo^s))
1728  !> divide by dx
1729  gradq(ixo^s) = gradq(ixo^s)/dxlevel(idir)
1730  else
1731  call mpistop("gradientO stencil unknown")
1732  endif
1733 
1734  end subroutine gradiento
1735 
1736 end module mod_afld
Module for including anisotropic flux limited diffusion (AFLD)-approximation in Radiation-hydrodynami...
Definition: mod_afld.t:8
logical lineforce_opacities
Use or dont use lineforce opacities.
Definition: mod_afld.t:70
subroutine afld_smooth_diffcoef(w, ixIL, ixOL, idim)
Use running average on Diffusion coefficient.
Definition: mod_afld.t:1329
character(len=8) fld_interaction_method
Which method to find the root for the energy interaction polynomial.
Definition: mod_afld.t:56
subroutine, public afld_radforce_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Definition: mod_afld.t:316
subroutine, public afld_get_fluxlimiter(w, x, ixIL, ixOL, fld_lambda, fld_R)
Calculate fld flux limiter This subroutine calculates flux limiter lambda using the prescription stor...
Definition: mod_afld.t:461
subroutine, public get_afld_energy_interact(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO This subroutine handles th...
Definition: mod_afld.t:347
subroutine put_diffterm_onegrid(ixIL, ixOL, w)
inplace update of psa==>F_im(psa)
Definition: mod_afld.t:1234
logical fld_eint_split
source split for energy interact and radforce:
Definition: mod_afld.t:15
subroutine bisection_method(e_gas, E_rad, c0, c1)
Find the root of the 4th degree polynomial using the bisection method.
Definition: mod_afld.t:1495
subroutine, public afld_init(He_abundance, rhd_radiation_diffusion, afld_gamma)
Initialising FLD-module: Read opacities Initialise Multigrid adimensionalise kappa Add extra variable...
Definition: mod_afld.t:126
subroutine evaluate_e_rad_mg(qtC, psa)
inplace update of psa==>F_im(psa)
Definition: mod_afld.t:1206
logical flux_lim_filter
Take a running average over the fluxlimiter.
Definition: mod_afld.t:66
integer, dimension(:), allocatable, public i_diff_mg
Index for Diffusion coeficients.
Definition: mod_afld.t:50
double precision, public fld_bisect_tol
Tolerance for bisection method for Energy sourceterms This is a percentage of the minimum of gas- and...
Definition: mod_afld.t:26
double precision, public fld_kappa0
Opacity per unit of unit_density.
Definition: mod_afld.t:19
double precision, public fld_diff_tol
Tolerance for adi method for radiative Energy diffusion.
Definition: mod_afld.t:29
character(len=16) fld_fluxlimiter
Diffusion limit lambda = 0.33.
Definition: mod_afld.t:45
subroutine afld_get_eddington(w, x, ixIL, ixOL, eddington_tensor)
Calculate Eddington-tensor Stores Eddington-tensor in w-array.
Definition: mod_afld.t:640
double precision function ddpolynomial_bisection(e_gas, c0, c1)
Evaluate second derivative of polynomial at argument e_gas.
Definition: mod_afld.t:1664
subroutine energy_interaction(w, wCT, x, ixIL, ixOL)
This subroutine calculates the radiation heating, radiation cooling and photon tiring using an implic...
Definition: mod_afld.t:1372
subroutine afld_get_diffcoef_central(w, wCT, x, ixIL, ixOL)
Calculates cell-centered diffusion coefficient to be used in multigrid.
Definition: mod_afld.t:1263
logical fld_diff_testcase
Set Diffusion coefficient to unity.
Definition: mod_afld.t:59
logical diffcrash_resume
Resume run when multigrid returns error.
Definition: mod_afld.t:73
double precision function dpolynomial_bisection(e_gas, c0, c1)
Evaluate first derivative of polynomial at argument e_gas.
Definition: mod_afld.t:1653
character(len=8), dimension(:), allocatable fld_opacity_law
Use constant Opacity?
Definition: mod_afld.t:41
subroutine update_diffcoeff(psa)
Definition: mod_afld.t:1307
subroutine halley_method(e_gas, E_rad, c0, c1)
Find the root of the 4th degree polynomial using the Halley method.
Definition: mod_afld.t:1605
double precision function polynomial_bisection(e_gas, c0, c1)
Evaluate polynomial at argument e_gas.
Definition: mod_afld.t:1642
character(len=50) fld_opal_table
Definition: mod_afld.t:42
character(len=8) afld_diff_scheme
Which method to solve diffusion part.
Definition: mod_afld.t:53
subroutine gradiento(q, x, ixIL, ixOL, idir, gradq, n)
Calculate gradient of a scalar q within ixL in direction idir difference with gradient is gradq(ixO^S...
Definition: mod_afld.t:1676
subroutine, public get_afld_rad_force(qdt, ixIL, ixOL, wCT, w, x, energy, qsourcesplit, active)
w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO This subroutine handles th...
Definition: mod_afld.t:208
subroutine diffuse_e_rad_mg(dtfactor, qdt, qtC, psa, psb)
Calling all subroutines to perform the multigrid method Communicates rad_e and diff_coeff to multigri...
Definition: mod_afld.t:1072
subroutine, public afld_get_radflux(w, x, ixIL, ixOL, rad_flux)
Calculate Radiation Flux stores radiation flux in w-array.
Definition: mod_afld.t:607
subroutine, public afld_get_opacity(w, x, ixIL, ixOL, fld_kappa)
Sets the opacity in the w-array by calling mod_opal_opacity.
Definition: mod_afld.t:371
double precision, public diff_crit
Number for splitting the diffusion module.
Definition: mod_afld.t:32
logical fld_radforce_split
Definition: mod_afld.t:16
integer size_d_filter
Definition: mod_afld.t:63
subroutine, public afld_get_radpress(w, x, ixIL, ixOL, rad_pressure)
Calculate Radiation Pressure Returns Radiation Pressure as tensor.
Definition: mod_afld.t:716
logical diff_coef_filter
Take running average for Diffusion coefficient.
Definition: mod_afld.t:62
double precision, public afld_mu
mean particle mass
Definition: mod_afld.t:22
integer size_l_filter
Definition: mod_afld.t:67
subroutine afld_params_read(files)
public methods these are called in mod_rhd_phys
Definition: mod_afld.t:100
integer, public i_test
Index for testvariable.
Definition: mod_afld.t:38
double precision dt_diff
running timestep for diffusion solver, initialised as zero
Definition: mod_afld.t:79
subroutine newton_method(e_gas, E_rad, c0, c1)
Find the root of the 4th degree polynomial using the Newton-Ralphson method.
Definition: mod_afld.t:1571
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
Module for physical and numeric constants.
Definition: mod_constants.t:2
double precision, parameter const_c
Definition: mod_constants.t:48
Module with basic grid data structures.
Definition: mod_forest.t:2
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:321
update ghost cells of all blocks including physical boundaries
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
This module contains definitions of global parameters and variables and some generic functions/subrou...
double precision small_pressure
double precision unit_density
Physical scaling factor for density.
double precision unit_opacity
Physical scaling factor for Opacity.
integer, parameter unitpar
file handle for IO
double precision global_time
The global simulation time.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer it
Number of time steps taken.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer mype
The rank of the current MPI task.
integer, dimension(:), allocatable, parameter d
double precision courantpar
The Courant (CFL) number used for the simulation.
double precision unit_velocity
Physical scaling factor for velocity.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
logical use_multigrid
Use multigrid (only available in 2D and 3D)
double precision dtmin
Stop the simulation when the time step becomes smaller than this value.
double precision, dimension(ndim) dxlevel
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
subroutine mg_copy_to_tree(iw_from, iw_to, restrict, restrict_gc, factor, state_from)
Copy a variable to the multigrid tree, including a layer of ghost cells.
subroutine mg_copy_from_tree_gc(iw_from, iw_to, state_to)
Copy from multigrid tree with one layer of ghost cells. Corner ghost cells are not used/set.
This module reads in Rosseland-mean opacities from OPAL tables. Table opacity values are given in bas...
subroutine, public init_opal_table(tablename)
This subroutine is called when the FLD radiation module is initialised. The OPAL tables for different...
subroutine, public set_opal_opacity(rho, temp, kappa)
This subroutine calculates the opacity for a given temperature-density structure. Opacities are read ...
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_get_tgas), pointer phys_get_tgas
Definition: mod_physics.t:78
procedure(sub_evaluate_implicit), pointer phys_evaluate_implicit
Definition: mod_physics.t:87
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl.
Definition: mod_physics.t:24
procedure(sub_get_pthermal), pointer phys_get_pthermal
Definition: mod_physics.t:77
procedure(sub_set_mg_bounds), pointer phys_set_mg_bounds
Definition: mod_physics.t:57
procedure(sub_implicit_update), pointer phys_implicit_update
Definition: mod_physics.t:86
Module with all the methods that users can customize in AMRVAC.
procedure(special_diffcoef), pointer usr_special_diffcoef
procedure(special_fluxlimiter), pointer usr_special_fluxlimiter
procedure(special_opacity_qdot), pointer usr_special_opacity_qdot
procedure(special_aniso_opacity), pointer usr_special_aniso_opacity
integer function var_set_extravar(name_cons, name_prim, ix)
Set extra variable in w, which is not advected and has no boundary conditions. This has to be done af...
The data structure that contains information about a tree node/grid block.
Definition: mod_forest.t:11