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