MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
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
8module mod_afld
9 use mod_comm_lib, only: mpistop
10 implicit none
11
12 !> source split for energy interact and radforce:
13 logical :: fld_eint_split = .false.
14 logical :: fld_radforce_split = .false.
15 !> Opacity per unit of unit_density
16 double precision, public :: fld_kappa0 = 0.34d0
17 !> mean particle mass
18 double precision, public :: afld_mu = 0.6d0
19 !> Tolerance for bisection method for Energy sourceterms
20 !> This is a percentage of the minimum of gas- and radiation energy
21 double precision, public :: fld_bisect_tol = 1.d-4
22 !> Tolerance for adi method for radiative Energy diffusion
23 double precision, public :: fld_diff_tol = 1.d-4
24 !> Number for splitting the diffusion module
25 double precision, public :: diff_crit
26 !> Use constant Opacity?
27 character(len=8), allocatable :: fld_opacity_law(:)
28 character(len=50) :: fld_opal_table = 'Y09800'
29 !> Diffusion limit lambda = 0.33
30 character(len=16) :: fld_fluxlimiter = 'Pomraning'
31 ! !> diffusion coefficient for multigrid method
32 integer, allocatable :: i_diff_mg(:)
33 !> Which method to solve diffusion part
34 character(len=8) :: afld_diff_scheme = 'mg'
35 !> Which method to find the root for the energy interaction polynomial
36 character(len=8) :: fld_interaction_method = 'Halley'
37 !> Set Diffusion coefficient to unity
38 logical :: fld_diff_testcase = .false.
39 !> Take running average for Diffusion coefficient
40 logical :: diff_coef_filter = .false.
41 integer :: size_d_filter = 1
42 !> Take a running average over the fluxlimiter
43 logical :: flux_lim_filter = .false.
44 integer :: size_l_filter = 1
45 !> Use or dont use lineforce opacities
46 logical :: lineforce_opacities = .false.
47 !> Resume run when multigrid returns error
48 logical :: diffcrash_resume = .true.
49 !> A copy of rhd_Gamma
50 double precision, private, protected :: afld_gamma
51 !> running timestep for diffusion solver, initialised as zero
52 double precision :: dt_diff = 0.d0
53 !> public methods
54 !> these are called in mod_rhd_phys
55 public :: get_afld_rad_force
57 public :: afld_radforce_get_dt
58 public :: afld_init
59 public :: afld_get_radflux
60 public :: afld_get_radpress
61 public :: afld_get_fluxlimiter
62 public :: afld_get_opacity
63 contains
64
65 !> Reading in fld-list parameters from .par file
66 subroutine afld_params_read(files)
69 character(len=*), intent(in) :: files(:)
70 integer :: n
71
72 namelist /fld_list/ fld_kappa0, fld_eint_split, fld_radforce_split, &
77
78 do n = 1, size(files)
79 open(unitpar, file=trim(files(n)), status="old")
80 read(unitpar, fld_list, end=111)
81 111 close(unitpar)
82 end do
83 end subroutine afld_params_read
84
85 !> Initialising FLD-module:
86 !> Read opacities
87 !> Initialise Multigrid
88 !> adimensionalise kappa
89 !> Add extra variables to w-array, flux, kappa, eddington Tensor
90 !> Lambda and R
91 !> ...
92 subroutine afld_init(He_abundance, rhd_radiation_diffusion, afld_gamma)
95 use mod_physics
98 double precision, intent(in) :: he_abundance, afld_gamma
99 logical, intent(in) :: rhd_radiation_diffusion
100 double precision :: sigma_thomson
101 integer :: idir,jdir
102 character(len=1) :: ind_1
103 character(len=1) :: ind_2
104 character(len=2) :: cmp_f
105 character(len=5) :: cmp_e
106
107 {^ifoned
108 call mpistop("Anisotropic in 1D... really?")
109 }
110
111 !> allocate boundary conditions
112 allocate(fld_opacity_law(1:ndim))
113 fld_opacity_law(1:ndim) = 'const'
114
115 !> read par files
117
118 if (rhd_radiation_diffusion) then
119 if (afld_diff_scheme .eq. 'mg') then
120 use_multigrid = .true.
123 mg%n_extra_vars = 1
124 mg%operator_type = mg_ahelmholtz
125 endif
126 endif
127 allocate(i_diff_mg(ndim))
128 do idir = 1,ndim
129 write(ind_1,'(I1)') idir
130 cmp_f = 'D' // ind_1
131 i_diff_mg(idir) = var_set_extravar(cmp_f,cmp_f)
132 enddo
133 !> Need mean molecular weight
134 afld_mu = (1.+4*he_abundance)/(2.+3.*he_abundance)
135 end subroutine afld_init
136
137 !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
138 !> This subroutine handles the radiation force
139 subroutine get_afld_rad_force(qdt,ixI^L,ixO^L,wCT,w,x,&
140 energy,qsourcesplit,active)
141 use mod_constants
144 use mod_geometry
145 integer, intent(in) :: ixi^l, ixo^l
146 double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
147 double precision, intent(in) :: wct(ixi^s,1:nw)
148 double precision, intent(inout) :: w(ixi^s,1:nw)
149 logical, intent(in) :: energy,qsourcesplit
150 logical, intent(inout) :: active
151 double precision :: radiation_forcect(ixo^s,1:ndim)
152 double precision :: kappact(ixo^s,1:ndim)
153 double precision :: rad_fluxct(ixo^s,1:ndim)
154 double precision :: div_v(ixi^s,1:ndim,1:ndim)
155 double precision :: edd(ixo^s,1:ndim,1:ndim)
156 double precision :: nabla_vp(ixo^s)
157 double precision :: vel(ixi^s), grad_v(ixi^s), grad0_v(ixo^s)
158 double precision :: grad_e(ixi^s)
159 integer :: idir, jdir
160 double precision :: afld_r(ixo^s,1:ndim), lambda(ixo^s,1:ndim)
161
162 !> Calculate and add sourceterms
163 if(qsourcesplit .eqv. fld_radforce_split) then
164 active = .true.
165 call afld_get_fluxlimiter(wct,x,ixi^l,ixo^l,lambda,afld_r,nghostcells)
166 do idir = 1,ndim
167 !> Radiation force = kappa*rho/c *Flux = lambda gradE
168 call gradient(wct(ixi^s,iw_r_e),ixi^l,ixo^l,idir,grad_e,nghostcells)
169 radiation_forcect(ixo^s,idir) = -lambda(ixo^s,idir)*grad_e(ixo^s)
170 !> Momentum equation source term
171 w(ixo^s,iw_mom(idir)) = w(ixo^s,iw_mom(idir)) &
172 + qdt*radiation_forcect(ixo^s,idir)
173 !> Energy equation source term (kinetic energy)
174 w(ixo^s,iw_e) = w(ixo^s,iw_e) &
175 + qdt*wct(ixo^s,iw_mom(idir))/wct(ixo^s,iw_rho)*radiation_forcect(ixo^s,idir)
176 enddo
177 !> Photon tiring
178 !> calculate tensor div_v
179 !> !$OMP PARALLEL DO
180 do idir = 1,ndim
181 do jdir = 1,ndim
182 vel(ixi^s) = wct(ixi^s,iw_mom(jdir))/wct(ixi^s,iw_rho)
183 call gradient(vel,ixi^l,ixo^l,idir,grad_v)
184 div_v(ixo^s,idir,jdir) = grad_v(ixo^s)
185 enddo
186 enddo
187 !> !$OMP END PARALLEL DO
188 call afld_get_eddington(wct, x, ixi^l, ixo^l, edd, nghostcells)
189 !> VARIABLE NAMES DIV ARE ACTUALLY GRADIENTS
190 {^ifoned
191 nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1)
192 }
193 {^iftwod
194 !>eq 34 Turner and stone (Only 2D)
195 nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1) &
196 + div_v(ixo^s,1,2)*edd(ixo^s,1,2) &
197 + div_v(ixo^s,2,1)*edd(ixo^s,2,1) &
198 + div_v(ixo^s,2,2)*edd(ixo^s,2,2)
199 }
200 {^ifthreed
201 nabla_vp(ixo^s) = div_v(ixo^s,1,1)*edd(ixo^s,1,1) &
202 + div_v(ixo^s,1,2)*edd(ixo^s,1,2) &
203 + div_v(ixo^s,1,3)*edd(ixo^s,1,3) &
204 + div_v(ixo^s,2,1)*edd(ixo^s,2,1) &
205 + div_v(ixo^s,2,2)*edd(ixo^s,2,2) &
206 + div_v(ixo^s,2,3)*edd(ixo^s,2,3) &
207 + div_v(ixo^s,3,1)*edd(ixo^s,3,1) &
208 + div_v(ixo^s,3,2)*edd(ixo^s,3,2) &
209 + div_v(ixo^s,3,3)*edd(ixo^s,3,3)
210 }
211 w(ixo^s,iw_r_e) = w(ixo^s,iw_r_e) &
212 - qdt * nabla_vp(ixo^s)*wct(ixo^s,iw_r_e)
213 end if
214 end subroutine get_afld_rad_force
215
216 subroutine afld_radforce_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
219 use mod_geometry
220 integer, intent(in) :: ixi^l, ixo^l
221 double precision, intent(in) :: dx^d, x(ixi^s,1:ndim), w(ixi^s,1:nw)
222 double precision, intent(inout) :: dtnew
223 integer :: idir
224 double precision :: radiation_force(ixo^s,1:ndim)
225 double precision :: lambda(ixo^s,1:ndim), grad_e(ixi^s)
226 double precision :: dxinv(1:ndim), afld_r(ixo^s,1:ndim)
227 double precision :: max_grad
228
229 ^d&dxinv(^d)=one/dx^d;
230 call afld_get_fluxlimiter(w,x,ixi^l,ixo^l,lambda,afld_r,nghostcells)
231 do idir = 1, ndim
232 call gradient(w(ixi^s,iw_r_e),ixi^l,ixo^l,idir,grad_e,nghostcells)
233 radiation_force(ixo^s,idir) = -lambda(ixo^s,idir)*grad_e(ixo^s)
234 max_grad = maxval(abs(radiation_force(ixo^s,idir)))
235 max_grad = max(max_grad, epsilon(1.0d0))
236 dtnew = min(dtnew, courantpar/sqrt(max_grad*dxinv(idir)))
237 end do
238 end subroutine afld_radforce_get_dt
239
240 !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
241 !> This subroutine handles the energy exchange between gas and radiation
242 subroutine get_afld_energy_interact(qdt,ixI^L,ixO^L,wCT,w,x,&
243 energy,qsourcesplit,active)
244 use mod_constants
247 use mod_physics, only: phys_get_pthermal !needed to get temp
248 integer, intent(in) :: ixi^l, ixo^l
249 double precision, intent(in) :: qdt, x(ixi^s,1:ndim)
250 double precision, intent(in) :: wct(ixi^s,1:nw)
251 double precision, intent(inout) :: w(ixi^s,1:nw)
252 logical, intent(in) :: energy,qsourcesplit
253 logical, intent(inout) :: active
254
255 !> Calculate and add sourceterms
256 if(qsourcesplit .eqv. fld_eint_split) then
257 active = .true.
258 !> Add energy sourceterms
259 call energy_interaction(w,w,x,ixi^l,ixo^l)
260 end if
261 end subroutine get_afld_energy_interact
262
263 !> Sets the opacity in the w-array
264 !> by calling mod_opal_opacity
265 subroutine afld_get_opacity(w, x, ixI^L, ixO^L, fld_kappa)
268 use mod_physics, only: phys_get_tgas
271 integer, intent(in) :: ixi^l, ixo^l
272 double precision, intent(in) :: w(ixi^s, 1:nw)
273 double precision, intent(in) :: x(ixi^s, 1:ndim)
274 double precision, intent(out) :: fld_kappa(ixo^s,1:ndim)
275 integer :: i,j,ix^d, idir
276 double precision :: temp(ixi^s), pth(ixi^s), a2(ixo^s)
277 double precision :: rho0,temp0,n,sigma_b
278 double precision :: akram, bkram
279 double precision :: vth(ixo^s), gradv(ixi^s), eta(ixo^s), t(ixo^s)
280
281 do idir = 1, ndim
282 select case (fld_opacity_law(idir))
283 case('const')
284 fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity
285 case('thomson')
286 fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity
287 case('kramers')
288 rho0 = half !> Take lower value of rho in domain
289 fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity*((w(ixo^s,iw_rho)/rho0))
290 case('bump')
291 !> Opacity bump
292 rho0 = 0.2d0
293 n = 7.d0
294 sigma_b = 2.d-2
295 fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity*(one + n*dexp(-one/sigma_b*(dlog(w(ixo^s,iw_rho)/rho0))**two))
296 case('non_iso')
297 call phys_get_pthermal(w,x,ixi^l,ixo^l,temp)
298 temp(ixo^s) = temp(ixo^s)/w(ixo^s,iw_rho)
299 rho0 = 0.5d0 !> Take lower value of rho in domain
300 temp0 = one
301 n = -7.d0/two
302 fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity*(w(ixo^s,iw_rho)/rho0)*(temp(ixo^s)/temp0)**n
303 case('fastwind')
304 call phys_get_pthermal(w,x,ixi^l,ixo^l,pth)
305 a2(ixo^s) = pth(ixo^s)/w(ixo^s,iw_rho)*unit_velocity**2.d0
306 akram = 13.1351597305d0
307 bkram = -4.5182188206d0
308 fld_kappa(ixo^s,idir) = fld_kappa0/unit_opacity &
309 * (1.d0+10.d0**akram*w(ixo^s,iw_rho)*unit_density*(a2(ixo^s)/1.d12)**bkram)
310 {do ix^d=ixomin^d,ixomax^d\ }
311 !> Hard limit on kappa
312 fld_kappa(ix^d,idir) = min(fld_kappa(ix^d,idir),2.3d0*fld_kappa0/unit_opacity)
313 {enddo\}
314 case('opal')
315 call phys_get_tgas(w,x,ixi^l,ixo^l,temp)
316 {do ix^d=ixomin^d,ixomax^d\ }
317 rho0 = w(ix^d,iw_rho)*unit_density
318 temp0 = temp(ix^d)*unit_temperature
319 call set_opal_opacity(rho0,temp0,n)
320 fld_kappa(ix^d,idir) = n/unit_opacity
321 {enddo\}
322 case('special')
323 if (.not. associated(usr_special_aniso_opacity)) then
324 call mpistop("special opacity not defined")
325 endif
326 call usr_special_aniso_opacity(ixi^l,ixo^l,w,x,fld_kappa(ixo^s,idir),idir)
327 case default
328 call mpistop("Doesn't know opacity law")
329 end select
330 end do
331 end subroutine afld_get_opacity
332
333 !> Calculate fld flux limiter
334 !> This subroutine calculates flux limiter lambda using the prescription
335 !> stored in fld_fluxlimiter.
336 !> It also calculates the ratio of radiation scaleheight and mean free path
337 subroutine afld_get_fluxlimiter(w,x,ixI^L,ixO^L,fld_lambda,fld_R,nth)
339 use mod_geometry
341 integer, intent(in) :: ixi^l, ixo^l, nth
342 double precision, intent(in) :: w(ixi^s, 1:nw)
343 double precision, intent(in) :: x(ixi^s, 1:ndim)
344 double precision, intent(out) :: fld_r(ixo^s,1:ndim), fld_lambda(ixo^s,1:ndim)
345 integer :: idir, ix^d
346 integer :: filter, idim
347 double precision :: kappa(ixo^s,1:ndim)
348 double precision :: normgrad2(ixo^s)
349 double precision :: grad_r_e(ixi^s)
350 double precision :: rad_e(ixi^s)
351 double precision :: tmp_l(ixi^s), filtered_l(ixi^s)
352
353 select case (fld_fluxlimiter)
354 case('Diffusion')
355 fld_lambda(ixo^s,1:ndim) = 1.d0/3.d0
356 fld_r(ixo^s,1:ndim) = zero
357 case('FreeStream')
358 !> Calculate R everywhere
359 !> |grad E|/(rho kappa E)
360 rad_e(ixi^s) = w(ixi^s, iw_r_e)
361 call afld_get_opacity(w,x,ixi^l,ixo^l,kappa)
362 do idir=1,ndim
363 call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e,nth)
364 normgrad2(ixo^s) = grad_r_e(ixo^s)**2
365 fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
366 !> Calculate the flux limiter, lambda
367 fld_lambda(ixo^s,idir) = one/fld_r(ixo^s,idir)
368 end do
369 case('Pomraning')
370 !> Calculate R everywhere
371 !> |grad E|/(rho kappa E)
372 normgrad2(ixo^s) = zero
373 rad_e(ixi^s) = w(ixi^s,iw_r_e)
374 call afld_get_opacity(w,x,ixi^l,ixo^l,kappa)
375 do idir = 1,ndim
376 call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e,nth)
377 normgrad2(ixo^s) = grad_r_e(ixo^s)**2
378 fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
379 !> Calculate the flux limiter, lambda
380 !> Levermore and Pomraning: lambda = (2 + R)/(6 + 3R + R^2)
381 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)
382 end do
383 case('Minerbo')
384 !> Calculate R everywhere
385 !> |grad E|/(rho kappa E)
386 normgrad2(ixo^s) = zero
387 rad_e(ixi^s) = w(ixi^s,iw_r_e)
388 call afld_get_opacity(w,x,ixi^l,ixo^l,kappa)
389 do idir = 1,ndim
390 call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e,nth)
391 normgrad2(ixo^s) = grad_r_e(ixo^s)**2
392 fld_r(ixo^s,idir) = dsqrt(normgrad2(ixo^s))/(kappa(ixo^s,idir)*w(ixo^s,iw_rho)*w(ixo^s,iw_r_e))
393 !> Calculate the flux limiter, lambda
394 !> Minerbo:
395 {do ix^d=ixomin^d,ixomax^d\ }
396 if(fld_r(ix^d,idir) .lt. 3.d0/2.d0) then
397 fld_lambda(ix^d,idir) = 2.d0/(3.d0 + dsqrt(9.d0 + 12.d0*fld_r(ix^d,idir)**2.d0))
398 else
399 fld_lambda(ix^d,idir) = 1.d0/(1.d0 + fld_r(ix^d,idir) + dsqrt(1.d0 + 2.d0*fld_r(ix^d,idir)))
400 endif
401 {enddo\}
402 end do
403 case('special')
404 if (.not. associated(usr_special_fluxlimiter)) then
405 call mpistop("special fluxlimiter not defined")
406 endif
407 call usr_special_fluxlimiter(ixi^l, ixo^l, w, x, fld_lambda, fld_r)
408 case default
409 call mpistop('Fluxlimiter unknown')
410 end select
411
412 if(flux_lim_filter) then
413 if(size_l_filter .lt. 1) call mpistop("D filter of size < 1 makes no sense")
414 if(size_l_filter .gt. nghostcells) call mpistop("D filter of size > nghostcells makes no sense")
415 do idir=1,ndim
416 tmp_l(ixo^s) = fld_lambda(ixo^s,idir)
417 filtered_l(ixo^s) = zero
418 do filter = 1,size_l_filter
419 {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_l_filter\}
420 do idim = 1,ndim
421 filtered_l(ix^d) = filtered_l(ix^d) &
422 + tmp_l(ix^d+filter*kr(idim,^d)) &
423 + tmp_l(ix^d-filter*kr(idim,^d))
424 enddo
425 {enddo\}
426 enddo
427 {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
428 tmp_l(ix^d) = (tmp_l(ix^d)+filtered_l(ix^d))/(1+2*size_l_filter*ndim)
429 {enddo\}
430 enddo
431 fld_lambda(ixo^s,idir) = tmp_l(ixo^s)
432 endif
433 end subroutine afld_get_fluxlimiter
434
435 !> Calculate Radiation Flux
436 !> Calculate the Flux using the fld closure relation
437 !> F = -c*lambda/(kappa*rho) *grad E
438 subroutine afld_get_radflux(w, x, ixI^L, ixO^L, rad_flux, nth)
440 use mod_geometry
441 integer, intent(in) :: ixi^l, ixo^l, nth
442 double precision, intent(in) :: w(ixi^s, 1:nw)
443 double precision, intent(in) :: x(ixi^s, 1:ndim)
444 double precision, intent(out) :: rad_flux(ixo^s, 1:ndim)
445 integer :: ix^d, idir
446 double precision :: cc
447 double precision :: grad_r_e(ixi^s)
448 double precision :: rad_e(ixi^s)
449 double precision :: kappa(ixo^s,1:ndim), lambda(ixo^s,1:ndim), fld_r(ixo^s,1:ndim)
450
451 cc = const_c/unit_velocity
452 rad_e(ixi^s) = w(ixi^s, iw_r_e)
453 call afld_get_opacity(w, x, ixi^l, ixo^l, kappa)
454 call afld_get_fluxlimiter(w, x, ixi^l, ixo^l, lambda, fld_r, nghostcells)
455 do idir = 1,ndim
456 call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e,nghostcells)
457 rad_flux(ixo^s, idir) = -cc*lambda(ixo^s,idir)/(kappa(ixo^s,idir)*w(ixo^s,iw_rho))*grad_r_e(ixo^s)
458 end do
459 end subroutine afld_get_radflux
460
461 !> Calculate Eddington-tensor
462 !> Stores Eddington-tensor in w-array
463 subroutine afld_get_eddington(w, x, ixI^L, ixO^L, eddington_tensor, nth)
465 use mod_geometry
466 integer, intent(in) :: ixI^L, ixO^L, nth
467 double precision, intent(in) :: w(ixI^S, 1:nw)
468 double precision, intent(in) :: x(ixI^S, 1:ndim)
469 double precision, intent(out) :: eddington_tensor(ixO^S,1:ndim,1:ndim)
470 double precision :: tnsr2(ixO^S,1:ndim,1:ndim)
471 double precision :: normgrad2(ixO^S), f(ixO^S,1:ndim)
472 double precision :: grad_r_e(ixI^S,1:ndim), rad_e(ixI^S)
473 double precision :: lambda(ixO^S,1:ndim), fld_R(ixO^S,1:ndim)
474 integer :: i,j, idir,jdir
475
476 call afld_get_fluxlimiter(w,x,ixi^l,ixo^l,lambda,fld_r, nth)
477 !> Calculate R everywhere
478 !> |grad E|/(rho kappa E)
479 normgrad2(ixo^s) = zero
480 rad_e(ixi^s) = w(ixi^s, iw_r_e)
481 do idir = 1,ndim
482 call gradient(rad_e,ixi^l,ixo^l,idir,grad_r_e(ixi^s,idir),nth)
483 normgrad2(ixo^s) = normgrad2(ixo^s)+grad_r_e(ixo^s,idir)**2
484 !> Calculate radiation pressure
485 !> P = (lambda + lambda^2 R^2)*E
486 f(ixo^s,idir) = lambda(ixo^s,idir)+lambda(ixo^s,idir)**two*fld_r(ixo^s,idir)**two
487 f(ixo^s,idir) = half*(one-f(ixo^s,idir))+half*(3.d0*f(ixo^s,idir)-one)
488 end do
489
490 {^ifoned
491 eddington_tensor(ixo^s,1,1) = f(ixo^s,1)
492 }
493 {^nooned
494 do idir = 1,ndim
495 eddington_tensor(ixo^s,idir,idir) = half*(one-f(ixo^s,idir))
496 enddo
497 do idir = 1,ndim
498 do jdir = 1,ndim
499 if(idir .ne. jdir) eddington_tensor(ixo^s,idir,jdir) = zero
500 tnsr2(ixo^s,idir,jdir) = half*(3.d0*f(ixo^s,idir) - 1)&
501 * grad_r_e(ixo^s,idir)*grad_r_e(ixo^s,jdir)/normgrad2(ixo^s)
502 enddo
503 enddo
504 do idir = 1,ndim
505 do jdir = 1,ndim
506 where((tnsr2(ixo^s,idir,jdir) .eq. tnsr2(ixo^s,idir,jdir)) &
507 .and. (normgrad2(ixo^s) .gt. smalldouble))
508 eddington_tensor(ixo^s,idir,jdir) = eddington_tensor(ixo^s,idir,jdir)+tnsr2(ixo^s,idir,jdir)
509 endwhere
510 enddo
511 enddo
512 }
513 end subroutine afld_get_eddington
514
515 !> Calculate Radiation Pressure
516 !> Returns Radiation Pressure as tensor
517 subroutine afld_get_radpress(w, x, ixI^L, ixO^L, rad_pressure, nth)
519 integer, intent(in) :: ixi^l, ixo^l, nth
520 double precision, intent(in) :: w(ixi^s, 1:nw)
521 double precision, intent(in) :: x(ixi^s, 1:ndim)
522 double precision, intent(out):: rad_pressure(ixo^s,1:ndim,1:ndim)
523 integer :: i,j
524 double precision :: eddington_tensor(ixo^s,1:ndim,1:ndim)
525
526 call afld_get_eddington(w,x,ixi^l,ixo^l,eddington_tensor, nth)
527 do i=1,ndim
528 do j=1,ndim
529 rad_pressure(ixo^s,i,j) = eddington_tensor(ixo^s,i,j)*w(ixo^s,iw_r_e)
530 enddo
531 enddo
532 end subroutine afld_get_radpress
533
534!> Calling all subroutines to perform the multigrid method
535!> Communicates rad_e and diff_coeff to multigrid library
536 subroutine diffuse_e_rad_mg(dtfactor,qdt,qtC,psa,psb)
538 use mod_forest
542 type(state), target :: psa(max_blocks)
543 type(state), target :: psb(max_blocks)
544 double precision, intent(in) :: qdt
545 double precision, intent(in) :: qtC
546 double precision, intent(in) :: dtfactor
547 integer, parameter :: max_its = 50
548 integer :: n
549 integer :: iw_to,iw_from
550 integer :: iigrid, igrid, id
551 integer :: nc, lvl, i
552 type(tree_node), pointer :: pnode
553 double precision :: fac, facD
554 double precision :: res, max_residual, lambda
555
556 !> Set diffusion timestep, add previous timestep if mg did not converge:
557 dt_diff = dt_diff + qdt
558 ! Avoid setting a very restrictive limit to the residual when the time step
559 ! is small (as the operator is ~ 1/(D * qdt))
560 if (qdt < dtmin) then
561 if(mype==0)then
562 print *,'skipping implicit solve: dt too small!'
563 print *,'Currently at time=',global_time,' time step=',qdt,' dtmin=',dtmin
564 endif
565 return
566 endif
567 ! max_residual = 1d-7/qdt
568 max_residual = fld_diff_tol !1d-7/qdt
569 mg%operator_type = mg_ahelmholtz
570 mg%smoother_type = mg_smoother_gsrb
571 call mg_set_methods(mg)
572 if (.not. mg%is_allocated) call mpistop("multigrid tree not allocated yet")
573 lambda = 1.d0/(dtfactor * dt_diff)
574 call ahelmholtz_set_lambda(lambda)
575 call update_diffcoeff(psb)
576 fac = 1.d0
577 facd = 1.d0
578 !This is mg_copy_to_tree from psb state
579 {^ifoned
580 call mg_copy_to_tree(i_diff_mg(1), mg_iveps, factor=facd, state_from=psb)
581 if (time_advance)then
582 call mg_restrict(mg, mg_iveps)
583 call mg_fill_ghost_cells(mg, mg_iveps)
584 endif
585 }
586 {^iftwod
587 call mg_copy_to_tree(i_diff_mg(1), mg_iveps1, factor=facd, state_from=psb)
588 call mg_copy_to_tree(i_diff_mg(2), mg_iveps2, factor=facd, state_from=psb)
589 if (time_advance)then
590 call mg_restrict(mg, mg_iveps1)
591 call mg_fill_ghost_cells(mg, mg_iveps1)
592
593 call mg_restrict(mg, mg_iveps2)
594 call mg_fill_ghost_cells(mg, mg_iveps2)
595 endif
596 }
597 !This is mg_copy_to_tree from psb state
598 call mg_copy_to_tree(iw_r_e, mg_iphi, factor=fac, state_from=psb)
599 !>replace call set_rhs(mg, -1/dt, 0.0_dp)
600 call mg_copy_to_tree(iw_r_e, mg_irhs, factor=-1/(dtfactor*dt_diff), state_from=psb)
601 call phys_set_mg_bounds()
602 call mg_fas_fmg(mg, .true., max_res=res)
603 do n = 1, max_its
604 ! print*, n, res
605 if (res < max_residual) exit
606 call mg_fas_vcycle(mg, max_res=res)
607 end do
608 if (res .le. 0.d0) then
609 if (diffcrash_resume) then
610 if (mg%my_rank == 0) &
611 write(*,*) it, ' resiudal zero ', res
612 return
613 endif
614 if (mg%my_rank == 0) then
615 print*, res
616 error stop "Diffusion residual to zero"
617 endif
618 endif
619 if(n == max_its + 1) then
620 if(diffcrash_resume) then
621 if(mg%my_rank == 0) &
622 write(*,*) it, ' resiudal high ', res
623 return
624 endif
625 if(mg%my_rank == 0) then
626 print *, "Did you specify boundary conditions correctly?"
627 print *, "Or is the variation in diffusion too large?"
628 print*, n, res
629 print *, mg%bc(1, mg_iphi)%bc_value, mg%bc(2, mg_iphi)%bc_value
630 end if
631 error stop "diffusion_solve: no convergence"
632 end if
633 !> Reset dt_diff when diffusion worked out
634 dt_diff = 0.d0
635 ! !This is mg_copy_from_tree_gc for psa state
636 call mg_copy_from_tree_gc(mg_iphi, iw_r_e, state_to=psa)
637 end subroutine diffuse_e_rad_mg
638
639 !> inplace update of psa==>F_im(psa)
640 subroutine evaluate_e_rad_mg(qtC,psa)
643 type(state), target :: psa(max_blocks)
644 double precision, intent(in) :: qtC
645 integer :: iigrid, igrid, level
646 integer :: ixO^L
647
648 call update_diffcoeff(psa)
649 ixo^l=ixg^ll^lsub1;
650 !$OMP PARALLEL DO PRIVATE(igrid)
651 do iigrid=1,igridstail; igrid=igrids(iigrid);
652 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
653 ! call afld_get_diffcoef_central(psa(igrid)%w, psa(igrid)%w, psa(igrid)%x, ixG^LL, ixO^L)
654 call put_diffterm_onegrid(ixg^ll,ixo^l,psa(igrid)%w)
655 end do
656 !$OMP END PARALLEL DO
657 ! enforce boundary conditions for psa
658 call getbc(qtc,0.d0,psa,1,nwflux+nwaux)
659 end subroutine evaluate_e_rad_mg
660
661 !> inplace update of psa==>F_im(psa)
662 subroutine put_diffterm_onegrid(ixI^L,ixO^L,w)
664 integer, intent(in) :: ixI^L, ixO^L
665 double precision, intent(inout) :: w(ixI^S, 1:nw)
666 double precision :: gradE(ixO^S),divF(ixO^S)
667 double precision :: divF_h(ixO^S),divF_j(ixO^S)
668 double precision :: diff_term(ixO^S)
669 integer :: idir, jxO^L, hxO^L
670
671 divf(ixo^s) = 0.d0
672 do idir = 1,ndim
673 hxo^l=ixo^l-kr(idir,^d);
674 jxo^l=ixo^l+kr(idir,^d);
675 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))
676 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))
677 divf(ixo^s) = divf(ixo^s) + 2.d0*(divf_h(ixo^s) + divf_j(ixo^s))/dxlevel(idir)**2
678 enddo
679 w(ixo^s,iw_r_e) = divf(ixo^s)
680 end subroutine put_diffterm_onegrid
681
682 !> Calculates cell-centered diffusion coefficient to be used in multigrid
683 subroutine afld_get_diffcoef_central(w, wCT, x, ixI^L, ixO^L)
685 use mod_geometry
687 integer, intent(in) :: ixI^L, ixO^L
688 double precision, intent(inout) :: w(ixI^S, 1:nw)
689 double precision, intent(in) :: wCT(ixI^S, 1:nw)
690 double precision, intent(in) :: x(ixI^S, 1:ndim)
691 double precision :: kappa(ixO^S,1:ndim), lambda(ixO^S,1:ndim), fld_R(ixO^S,1:ndim)
692 double precision :: max_D(ixI^S), grad_r_e(ixI^S), rad_e(ixI^S)
693 double precision :: cc
694 integer :: idir,i,j,ix^D
695
696 cc = const_c/unit_velocity
697 if(fld_diff_testcase) then
698 w(ixo^s,i_diff_mg(:)) = 1.d0
699 else
700 call afld_get_opacity(wct, x, ixi^l, ixo^l, kappa)
701 call afld_get_fluxlimiter(wct, x, ixi^l, ixo^l, lambda, fld_r, nghostcells-1)
702 do idir = 1,ndim
703 !> calculate diffusion coefficient
704 w(ixo^s,i_diff_mg(idir)) = cc*lambda(ixo^s,idir)/(kappa(ixo^s,idir)*wct(ixo^s,iw_rho))
705 where(w(ixo^s,i_diff_mg(idir)) .lt. 0.d0)
706 w(ixo^s,i_diff_mg(idir)) = smalldouble
707 end where
708 if(diff_coef_filter) then
709 call afld_smooth_diffcoef(w,ixi^l,ixo^l,idir)
710 endif
711 enddo
712 endif
713 if(associated(usr_special_diffcoef)) &
714 call usr_special_diffcoef(w, wct, x, ixi^l, ixo^l)
715 end subroutine afld_get_diffcoef_central
716
717 subroutine update_diffcoeff(psa)
719 type(state), target :: psa(max_blocks)
720 integer :: iigrid, igrid, level
721 integer :: ixO^L
722
723 ixo^l=ixg^ll^lsub1;
724 !$OMP PARALLEL DO PRIVATE(igrid)
725 do iigrid=1,igridstail; igrid=igrids(iigrid);
726 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
727 call afld_get_diffcoef_central(psa(igrid)%w, psa(igrid)%w, psa(igrid)%x, ixg^ll, ixo^l)
728 end do
729 !$OMP END PARALLEL DO
730 end subroutine update_diffcoeff
731
732 !> Use running average on Diffusion coefficient
733 subroutine afld_smooth_diffcoef(w, ixI^L, ixO^L,idir)
735 integer, intent(in) :: ixI^L, ixO^L, idir
736 double precision, intent(inout) :: w(ixI^S, 1:nw)
737 double precision :: tmp_D(ixI^S), filtered_D(ixI^S)
738 integer :: ix^D, filter, idim
739
740 if(size_d_filter .lt. 1) call mpistop("D filter of size < 1 makes no sense")
741 if(size_d_filter .gt. nghostcells) call mpistop("D filter of size > nghostcells makes no sense")
742 tmp_d(ixo^s) = w(ixo^s,i_diff_mg(idir))
743 filtered_d(ixo^s) = zero
744 do filter = 1,size_d_filter
745 {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
746 do idim = 1,ndim
747 filtered_d(ix^d) = filtered_d(ix^d) &
748 + tmp_d(ix^d+filter*kr(idim,^d)) &
749 + tmp_d(ix^d-filter*kr(idim,^d))
750 enddo
751 {enddo\}
752 enddo
753 {do ix^d = ixomin^d+size_d_filter,ixomax^d-size_d_filter\}
754 tmp_d(ix^d) = (tmp_d(ix^d)+filtered_d(ix^d))/(1+2*size_d_filter*ndim)
755 {enddo\}
756 w(ixo^s,i_diff_mg(idir)) = tmp_d(ixo^s)
757 end subroutine afld_smooth_diffcoef
758
759 !> This subroutine calculates the radiation heating, radiation cooling
760 !> and photon tiring using an implicit scheme.
761 !> These sourceterms are applied using the root-finding of a 4th order polynomial
762 !> This routine loops over every cell in the domain
763 !> and computes the coefficients of the polynomials in every cell
764 subroutine energy_interaction(w, wCT, x, ixI^L, ixO^L)
766 use mod_geometry
767 use mod_physics
769 integer, intent(in) :: ixI^L, ixO^L
770 double precision, intent(in) :: x(ixI^S,1:ndim)
771 double precision, intent(in) :: wCT(ixI^S,1:nw)
772 double precision, intent(inout) :: w(ixI^S,1:nw)
773 double precision :: a1(ixO^S), a2(ixO^S)
774 double precision :: c0(ixO^S), c1(ixO^S)
775 double precision :: e_gas(ixO^S), E_rad(ixO^S)
776 double precision :: kappa(ixO^S)
777 double precision :: sigma_b,cc
778 integer :: i,j,ix^D
779
780 cc = const_c/unit_velocity
781 sigma_b = const_rad_a*cc/4.d0*(unit_temperature**4.d0)/(unit_pressure)
782 !> e_gas is the INTERNAL ENERGY without KINETIC ENERGY
783 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)
784 if(allocated(iw_mag)) then
785 e_gas(ixo^s) = e_gas(ixo^s)-half*sum(wct(ixo^s,iw_mag(:))**2,dim=ndim+1)
786 endif
787 {do ix^d=ixomin^d,ixomax^d\ }
788 e_gas(ix^d) = max(e_gas(ix^d),small_pressure/(afld_gamma-1.d0))
789 {enddo\}
790 e_rad(ixo^s) = wct(ixo^s,iw_r_e)
791 if(associated(usr_special_opacity_qdot)) then
792 call usr_special_opacity_qdot(ixi^l,ixo^l,w,x,kappa)
793 else
794 call afld_get_opacity(wct,x,ixi^l,ixo^l,kappa)
795 endif
796 if(fld_interaction_method .eq. 'Instant') then
797 a1(ixo^s) = 4.d0*sigma_b/cc*(afld_gamma-one)**4.d0/wct(ixo^s,iw_rho)**4.d0
798 a2(ixo^s) = e_gas(ixo^s) + e_rad(ixo^s)
799 c0(ixo^s) = a2(ixo^s)/a1(ixo^s)
800 c1(ixo^s) = 1.d0/a1(ixo^s)
801 else
802 !> Calculate coefficients for polynomial
803 a1(ixo^s) = 4.d0*kappa(ixo^s)*sigma_b*(afld_gamma-one)**4.d0/wct(ixo^s,iw_rho)**3.d0*dt
804 a2(ixo^s) = (const_c/unit_velocity)*kappa(ixo^s)*wct(ixo^s,iw_rho)*dt
805 c0(ixo^s) = ((one+a2(ixo^s))*e_gas(ixo^s)+a2(ixo^s)*e_rad(ixo^s))/a1(ixo^s)
806 c1(ixo^s) = (one+a2(ixo^s))/a1(ixo^s)
807 endif
808 !> Loop over every cell for rootfinding method
809 {do ix^d=ixomin^d,ixomax^d\ }
810 select case(fld_interaction_method)
811 case('Bisect')
812 call bisection_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
813 case('Newton')
814 call newton_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
815 case('Halley')
816 call halley_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
817 case('Instant')
818 call halley_method(e_gas(ix^d), e_rad(ix^d), c0(ix^d), c1(ix^d))
819 case default
820 call mpistop('root-method not known')
821 end select
822 {enddo\}
823
824 if(fld_interaction_method .eq. 'Instant') then
825 e_rad(ixo^s) = a2(ixo^s)-e_gas(ixo^s)
826 else
827 !> advance E_rad
828 e_rad(ixo^s) = (a1(ixo^s)*e_gas(ixo^s)**4.d0+e_rad(ixo^s))/(one + a2(ixo^s))
829 endif
830 !> Update gas-energy in w, internal + kinetic
831 w(ixo^s,iw_e) = e_gas(ixo^s)
832 !> Beginning of module substracted wCT Ekin
833 !> So now add wCT Ekin
834 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)
835 if(allocated(iw_mag)) then
836 w(ixo^s,iw_e) = w(ixo^s,iw_e)+half*sum(wct(ixo^s, iw_mag(:))**2,dim=ndim+1)
837 endif
838 !> Update rad-energy in w
839 w(ixo^s,iw_r_e) = e_rad(ixo^s)
840 end subroutine energy_interaction
841
842 !> Find the root of the 4th degree polynomial using the bisection method
843 subroutine bisection_method(e_gas, E_rad, c0, c1)
845 double precision, intent(in) :: c0, c1
846 double precision, intent(in) :: E_rad
847 double precision, intent(inout) :: e_gas
848 double precision :: bisect_a, bisect_b, bisect_c
849 integer :: n, max_its
850
851 n = 0
852 max_its = 10000000
853 bisect_a = zero
854 bisect_b = 1.2d0*max(abs(c0/c1),abs(c0)**(1.d0/4.d0))
855 do while(abs(bisect_b-bisect_a) .ge. fld_bisect_tol*min(e_gas,e_rad))
856 bisect_c = (bisect_a + bisect_b)/two
857 n = n +1
858 if(n .gt. max_its) then
859 goto 2435
860 call mpistop('No convergece in bisection scheme')
861 endif
862 if(polynomial_bisection(bisect_a, c0, c1)*&
863 polynomial_bisection(bisect_b, c0, c1) .lt. zero) then
864 if(polynomial_bisection(bisect_a, c0, c1)*&
865 polynomial_bisection(bisect_c, c0, c1) .lt. zero) then
866 bisect_b = bisect_c
867 elseif(polynomial_bisection(bisect_b, c0, c1)*&
868 polynomial_bisection(bisect_c, c0, c1) .lt. zero) then
869 bisect_a = bisect_c
870 elseif(polynomial_bisection(bisect_a, c0, c1) .eq. zero) then
871 bisect_b = bisect_a
872 bisect_c = bisect_a
873 goto 2435
874 elseif(polynomial_bisection(bisect_b, c0, c1) .eq. zero) then
875 bisect_a = bisect_b
876 bisect_c = bisect_b
877 goto 2435
878 elseif(polynomial_bisection(bisect_c, c0, c1) .eq. zero) then
879 bisect_a = bisect_c
880 bisect_b = bisect_c
881 goto 2435
882 else
883 call mpistop("Problem with fld bisection method")
884 endif
885 elseif(polynomial_bisection(bisect_a, c0, c1) &
886 - polynomial_bisection(bisect_b, c0, c1) .lt. fld_bisect_tol*polynomial_bisection(bisect_a, c0, c1)) then
887 goto 2435
888 else
889 bisect_a = e_gas
890 bisect_b = e_gas
891 print*, "IGNORING GAS-RAD ENERGY EXCHANGE ", c0, c1
892 print*, polynomial_bisection(bisect_a, c0, c1), polynomial_bisection(bisect_b, c0, c1)
893 if(polynomial_bisection(bisect_a, c0, c1) .le. smalldouble) then
894 bisect_b = bisect_a
895 elseif(polynomial_bisection(bisect_a, c0, c1) .le. smalldouble) then
896 bisect_a = bisect_b
897 endif
898 goto 2435
899 endif
900 enddo
901 2435 e_gas = (bisect_a + bisect_b)/two
902 end subroutine bisection_method
903
904 !> Find the root of the 4th degree polynomial using the Newton-Ralphson method
905 subroutine newton_method(e_gas, E_rad, c0, c1)
907 double precision, intent(in) :: c0, c1
908 double precision, intent(in) :: E_rad
909 double precision, intent(inout) :: e_gas
910 double precision :: xval, yval, der, deltax
911 integer :: ii
912
913 yval = bigdouble
914 xval = e_gas
915 der = one
916 deltax = one
917 ii = 0
918 !> Compare error with dx = dx/dy dy
919 do while (abs(deltax) .gt. fld_bisect_tol)
920 yval = polynomial_bisection(xval, c0, c1)
921 der = dpolynomial_bisection(xval, c0, c1)
922 deltax = -yval/der
923 xval = xval + deltax
924 ii = ii + 1
925 if (ii .gt. 1d3) then
926 call bisection_method(e_gas, e_rad, c0, c1)
927 return
928 endif
929 enddo
930 e_gas = xval
931 end subroutine newton_method
932
933 !> Find the root of the 4th degree polynomial using the Halley method
934 subroutine halley_method(e_gas, E_rad, c0, c1)
936 double precision, intent(in) :: c0, c1
937 double precision, intent(in) :: E_rad
938 double precision, intent(inout) :: e_gas
939 double precision :: xval, yval, der, dder, deltax
940 integer :: ii
941
942 yval = bigdouble
943 xval = e_gas
944 der = one
945 dder = one
946 deltax = one
947 ii = 0
948 !> Compare error with dx = dx/dy dy
949 do while (abs(deltax) .gt. fld_bisect_tol)
950 yval = polynomial_bisection(xval, c0, c1)
951 der = dpolynomial_bisection(xval, c0, c1)
952 dder = ddpolynomial_bisection(xval, c0, c1)
953 deltax = -two*yval*der/(two*der**2 - yval*dder)
954 xval = xval + deltax
955 ii = ii + 1
956 if (ii .gt. 1d3) then
957 call newton_method(e_gas, e_rad, c0, c1)
958 return
959 endif
960 enddo
961 e_gas = xval
962 end subroutine halley_method
963
964 !> Evaluate polynomial at argument e_gas
965 function polynomial_bisection(e_gas, c0, c1) result(val)
967 double precision, intent(in) :: e_gas
968 double precision, intent(in) :: c0, c1
969 double precision :: val
970
971 val = e_gas**4.d0 + c1*e_gas - c0
972 end function polynomial_bisection
973
974 !> Evaluate first derivative of polynomial at argument e_gas
975 function dpolynomial_bisection(e_gas, c0, c1) result(der)
977 double precision, intent(in) :: e_gas
978 double precision, intent(in) :: c0, c1
979 double precision :: der
980
981 der = 4.d0*e_gas**3.d0 + c1
982 end function dpolynomial_bisection
983
984 !> Evaluate second derivative of polynomial at argument e_gas
985 function ddpolynomial_bisection(e_gas, c0, c1) result(dder)
987 double precision, intent(in) :: e_gas
988 double precision, intent(in) :: c0, c1
989 double precision :: dder
990
991 dder = 4.d0*3.d0*e_gas**2.d0
992 end function ddpolynomial_bisection
993end 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:46
subroutine afld_smooth_diffcoef(w, ixil, ixol, idir)
Use running average on Diffusion coefficient.
Definition mod_afld.t:734
character(len=8) fld_interaction_method
Which method to find the root for the energy interaction polynomial.
Definition mod_afld.t:36
logical fld_eint_split
source split for energy interact and radforce:
Definition mod_afld.t:13
subroutine afld_get_diffcoef_central(w, wct, x, ixil, ixol)
Calculates cell-centered diffusion coefficient to be used in multigrid.
Definition mod_afld.t:684
logical flux_lim_filter
Take a running average over the fluxlimiter.
Definition mod_afld.t:43
subroutine, public afld_get_radflux(w, x, ixil, ixol, rad_flux, nth)
Calculate Radiation Flux Calculate the Flux using the fld closure relation F = -c*lambda/(kappa*rho) ...
Definition mod_afld.t:439
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:21
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:141
double precision, public fld_kappa0
Opacity per unit of unit_density.
Definition mod_afld.t:16
subroutine, public afld_get_fluxlimiter(w, x, ixil, ixol, fld_lambda, fld_r, nth)
Calculate fld flux limiter This subroutine calculates flux limiter lambda using the prescription stor...
Definition mod_afld.t:338
double precision, public fld_diff_tol
Tolerance for adi method for radiative Energy diffusion.
Definition mod_afld.t:23
character(len=16) fld_fluxlimiter
Diffusion limit lambda = 0.33.
Definition mod_afld.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_afld.t:537
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:93
double precision function ddpolynomial_bisection(e_gas, c0, c1)
Evaluate second derivative of polynomial at argument e_gas.
Definition mod_afld.t:986
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:844
logical fld_diff_testcase
Set Diffusion coefficient to unity.
Definition mod_afld.t:38
logical diffcrash_resume
Resume run when multigrid returns error.
Definition mod_afld.t:48
double precision function dpolynomial_bisection(e_gas, c0, c1)
Evaluate first derivative of polynomial at argument e_gas.
Definition mod_afld.t:976
character(len=8), dimension(:), allocatable fld_opacity_law
Use constant Opacity?
Definition mod_afld.t:27
subroutine, public afld_radforce_get_dt(w, ixil, ixol, dtnew, dxd, x)
Definition mod_afld.t:217
subroutine evaluate_e_rad_mg(qtc, psa)
inplace update of psa==>F_im(psa)
Definition mod_afld.t:641
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:765
subroutine, public afld_get_radpress(w, x, ixil, ixol, rad_pressure, nth)
Calculate Radiation Pressure Returns Radiation Pressure as tensor.
Definition mod_afld.t:518
subroutine update_diffcoeff(psa)
Definition mod_afld.t:718
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:935
double precision function polynomial_bisection(e_gas, c0, c1)
Evaluate polynomial at argument e_gas.
Definition mod_afld.t:966
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:906
character(len=50) fld_opal_table
Definition mod_afld.t:28
character(len=8) afld_diff_scheme
Which method to solve diffusion part.
Definition mod_afld.t:34
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:244
double precision, public diff_crit
Number for splitting the diffusion module.
Definition mod_afld.t:25
logical fld_radforce_split
Definition mod_afld.t:14
integer size_d_filter
Definition mod_afld.t:41
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:266
logical diff_coef_filter
Take running average for Diffusion coefficient.
Definition mod_afld.t:40
double precision, public afld_mu
mean particle mass
Definition mod_afld.t:18
integer, dimension(:), allocatable i_diff_mg
Definition mod_afld.t:32
integer size_l_filter
Definition mod_afld.t:44
subroutine put_diffterm_onegrid(ixil, ixol, w)
inplace update of psa==>F_im(psa)
Definition mod_afld.t:663
subroutine afld_params_read(files)
public methods these are called in mod_rhd_phys
Definition mod_afld.t:67
double precision dt_diff
running timestep for diffusion solver, initialised as zero
Definition mod_afld.t:52
subroutine afld_get_eddington(w, x, ixil, ixol, eddington_tensor, nth)
Calculate Eddington-tensor Stores Eddington-tensor in w-array.
Definition mod_afld.t:464
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
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, nth_in)
update ghost cells of all blocks including physical boundaries
subroutine getbc(time, qdt, psb, nwstart, nwbc)
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 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.
double precision, dimension(:), allocatable, parameter d
double precision dt
global time step
double precision courantpar
The Courant (CFL) number used for the simulation.
double precision unit_velocity
Physical scaling factor for velocity.
logical time_advance
do time evolving
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.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
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.
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:74
procedure(sub_evaluate_implicit), pointer phys_evaluate_implicit
Definition mod_physics.t:83
procedure(sub_get_pthermal), pointer phys_get_pthermal
Definition mod_physics.t:73
procedure(sub_set_mg_bounds), pointer phys_set_mg_bounds
Definition mod_physics.t:53
procedure(sub_implicit_update), pointer phys_implicit_update
Definition mod_physics.t:82
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