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