MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_cak_force.t
Go to the documentation of this file.
1!> Module to include CAK radiation line force in (magneto)hydrodynamic models
2!> Computes both the force from free electrons and the force from an ensemble of
3!> lines (various possibilities for the latter).
4!> There is an option to only simulate the pure radial CAK force (with various
5!> corrections applied) as well as the full vector CAK force. Depending on the
6!> chosen option additional output are the CAK line force component(s) and,
7!> when doing a 1-D radial force, the finite disc factor.
8!>
9!> USAGE:
10!>
11!> 1. Include a cak_list in the .par file and activate (m)hd_cak_force in the
12!> (m)hd_list
13!> 2. Create a mod_usr.t file for the problem with appropriate initial and
14!> boundary conditions
15!> 3. In the mod_usr.t header call the mod_cak_force module to have access to
16!> global variables from mod_cak_force, which may be handy for printing or
17!> the computation of other variables inside mod_usr.t
18!> 4. In usr_init of mod_usr.t call the set_cak_force_norm routine and pass
19!> along the stellar radius and wind temperature---this is needed to
20!> correctly compute the (initial) force normalisation inside mod_cak_force
21!> 5. Ensure that the order of calls in usr_init is similar as for test problem
22!> CAKwind_spherical_1D: first reading usr.par list; then set unit scales;
23!> then call (M)HD_activate; then call set_cak_force_norm. This order avoids
24!> an incorrect force normalisation and code crash
25!>
26!> Developed by Florian Driessen (2022)
29 implicit none
30
31 !> Line-ensemble parameters in the Gayley (1995) formalism
32 real(8), public :: cak_alpha, gayley_qbar, gayley_q0
33
34 !> Ray positions + weights for impact parameter and azimuthal radiation angle
35 real(8), allocatable, private :: ay(:), wy(:), aphi(:), wphi(:)
36
37 !> The adiabatic index
38 real(8), private :: cak_gamma
39
40 !> Variables needed to compute force normalisation fnorm in initialisation
41 real(8), private :: lum, dlum, drstar, dke, dclight
42
43 !> To enforce a floor temperature when doing adiabatic (M)HD
44 real(8), private :: tfloor
45
46 !> Switch to choose between the 1-D CAK line force options
47 integer :: cak_1d_opt
48
49 ! Avoid magic numbers in code for 1-D CAK line force option
50 integer, parameter, private :: radstream=0, fdisc=1, fdisc_cutoff=2
51
52 !> Amount of rays in radiation polar and radiation azimuthal direction
53 integer :: nthetaray, nphiray
54
55 !> Extra slots to store quantities in w-array
56 integer :: gcak1_, gcak2_, gcak3_, fdf_
57
58 !> To treat source term in split or unsplit (default) fashion
59 logical :: cak_split=.false.
60
61 !> To activate the original CAK 1-D line force computation
62 logical :: cak_1d_force=.false.
63
64 !> To activate the vector CAK line force computation
65 logical :: cak_vector_force=.false.
66
67 !> To activate the pure radial vector CAK line force computation
68 logical :: fix_vector_force_1d=.false.
69
70 !> Public method
71 public :: set_cak_force_norm
72
73contains
74
75 !> Read this module's parameters from a file
76 subroutine cak_params_read(files)
78
79 character(len=*), intent(in) :: files(:)
80
81 ! Local variable
82 integer :: n
83
84 namelist /cak_list/ cak_alpha, gayley_qbar, gayley_q0, cak_1d_opt, &
87
88 do n = 1,size(files)
89 open(unitpar, file=trim(files(n)), status="old")
90 read(unitpar, cak_list, end=111)
91 111 close(unitpar)
92 enddo
93
94 end subroutine cak_params_read
95
96 !> Initialize the module
97 subroutine cak_init(phys_gamma)
99 use mod_comm_lib, only: mpistop
100
101 real(8), intent(in) :: phys_gamma
102
103 cak_gamma = phys_gamma
104
105 ! Set some defaults when user does not
106 cak_alpha = 0.65d0
107 gayley_qbar = 2000.0d0
108 gayley_q0 = 2000.0d0
109 cak_1d_opt = fdisc
110 nthetaray = 6
111 nphiray = 6
112
114
115 if (cak_1d_force) then
116 gcak1_ = var_set_extravar("gcak1", "gcak1")
117 fdf_ = var_set_extravar("fdfac", "fdfac")
118 endif
119
120 if (cak_vector_force) then
121 gcak1_ = var_set_extravar("gcak1", "gcak1")
122 gcak2_ = var_set_extravar("gcak2", "gcak2")
123 gcak3_ = var_set_extravar("gcak3", "gcak3")
125 endif
126
127 ! Some sanity checks
128 if ((cak_alpha <= 0.0d0) .or. (cak_alpha > 1.0d0)) then
129 call mpistop('CAK error: choose alpha in [0,1[')
130 endif
131
132 if ((gayley_qbar < 0.0d0) .or. (gayley_q0 < 0.0d0)) then
133 call mpistop('CAK error: chosen Qbar or Q0 is < 0')
134 endif
135
136 if (cak_1d_force .and. cak_vector_force) then
137 call mpistop('CAK error: choose either 1-D or vector force')
138 endif
139
140 end subroutine cak_init
141
142 !> Compute some (unitless) variables for CAK force normalisation
143 subroutine set_cak_force_norm(rstar,twind)
145 use mod_constants
146
147 real(8), intent(in) :: rstar, twind
148
149 lum = 4.0d0*dpi * rstar**2.0d0 * const_sigma * twind**4.0d0
151 dclight = const_c/unit_velocity
152 dlum = lum/(unit_density * unit_length**5.0d0 / unit_time**3.0d0)
153 drstar = rstar/unit_length
154 tfloor = twind/unit_temperature
155
156 end subroutine set_cak_force_norm
157
158 !> w[iw]=w[iw]+qdt*S[wCT,qtC,x] where S is the source based on wCT within ixO
159 subroutine cak_add_source(qdt,ixI^L,ixO^L,wCT,w,x,energy,qsourcesplit,active)
161 use mod_comm_lib, only: mpistop
162
163 integer, intent(in) :: ixI^L, ixO^L
164 real(8), intent(in) :: qdt, x(ixI^S,1:ndim), wCT(ixI^S,1:nw)
165 real(8), intent(inout) :: w(ixI^S,1:nw)
166 logical, intent(in) :: energy, qsourcesplit
167 logical, intent(inout) :: active
168
169 ! Local variables
170 real(8) :: gl(ixO^S,1:3), ge(ixO^S), ptherm(ixI^S), pmin(ixI^S)
171 integer :: idir
172
173 ! By default add source in unsplit fashion together with the fluxes
174 if (qsourcesplit .eqv. cak_split) then
175
176 ! Thomson force
177 call get_gelectron(ixi^l,ixo^l,wct,x,ge)
178
179 ! CAK line force
180 gl(ixo^s,1:3) = 0.0d0
181
182 if (cak_1d_force) then
183 call get_cak_force_radial(ixi^l,ixo^l,wct,w,x,gl)
184 elseif (cak_vector_force) then
185 call get_cak_force_vector(ixi^l,ixo^l,wct,w,x,gl)
186 else
187 call mpistop("No valid force option")
188 endif
189
190 ! Update conservative vars: w = w + qdt*gsource
191 do idir = 1,ndir
192 if (idir == 1) gl(ixo^s,idir) = gl(ixo^s,idir) + ge(ixo^s)
193
194 w(ixo^s,iw_mom(idir)) = w(ixo^s,iw_mom(idir)) &
195 + qdt * gl(ixo^s,idir) * wct(ixo^s,iw_rho)
196
197 if (energy) then
198 w(ixo^s,iw_e) = w(ixo^s,iw_e) + qdt * gl(ixo^s,idir) * wct(ixo^s,iw_mom(idir))
199 endif
200 enddo
201
202 ! Impose fixed floor temperature to mimic stellar heating
203 if (energy) then
204 call phys_get_pthermal(w,x,ixi^l,ixo^l,ptherm)
205 pmin(ixo^s) = w(ixo^s,iw_rho) * tfloor
206
207 where (ptherm(ixo^s) < pmin(ixo^s))
208 w(ixo^s,iw_e) = w(ixo^s,iw_e) + (pmin(ixo^s) - ptherm(ixo^s))/(cak_gamma - 1.0d0)
209 endwhere
210 endif
211 endif
212
213 end subroutine cak_add_source
214
215 !> 1-D CAK line force in the Gayley line-ensemble distribution parametrisation
216 subroutine get_cak_force_radial(ixI^L,ixO^L,wCT,w,x,gcak)
218 use mod_comm_lib, only: mpistop
219
220 integer, intent(in) :: ixI^L, ixO^L
221 real(8), intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
222 real(8), intent(inout) :: w(ixI^S,1:nw)
223 real(8), intent(inout) :: gcak(ixO^S,1:3)
224
225 ! Local variables
226 real(8) :: vr(ixI^S), dvrdr(ixO^S)
227 real(8) :: beta_fd(ixO^S), fdfac(ixO^S), taus(ixO^S), ge(ixO^S)
228
229 vr(ixi^s) = wct(ixi^s,iw_mom(1)) / wct(ixi^s,iw_rho)
230 call get_velocity_gradient(ixi^l,ixo^l,vr,x,1,dvrdr)
231
232 if (physics_type == 'hd') then
233 ! Monotonic flow to avoid multiple resonances and radiative coupling
234 dvrdr(ixo^s) = abs(dvrdr(ixo^s))
235 else
236 ! Allow material to fallback to the star in a magnetosphere model
237 dvrdr(ixo^s) = max(dvrdr(ixo^s), smalldouble)
238 endif
239
240 ! Thomson force
241 call get_gelectron(ixi^l,ixo^l,wct,x,ge)
242
243 ! Sobolev optical depth for line ensemble (tau = Qbar * t_r) and the force
244 select case (cak_1d_opt)
245 case(radstream, fdisc)
246 taus(ixo^s) = gayley_qbar * dke * dclight * wct(ixo^s,iw_rho)/dvrdr(ixo^s)
247 gcak(ixo^s,1) = gayley_qbar/(1.0d0 - cak_alpha) &
248 * ge(ixo^s)/taus(ixo^s)**cak_alpha
249
250 case(fdisc_cutoff)
251 taus(ixo^s) = gayley_q0 * dke * dclight * wct(ixo^s,iw_rho)/dvrdr(ixo^s)
252 gcak(ixo^s,1) = gayley_qbar * ge(ixo^s) &
253 * ( (1.0d0 + taus(ixo^s))**(1.0d0 - cak_alpha) - 1.0d0 ) &
254 / ( (1.0d0 - cak_alpha) * taus(ixo^s) )
255 case default
256 call mpistop("Error in force computation.")
257 end select
258
259 ! Finite disk factor parameterisation (Owocki & Puls 1996)
260 beta_fd(ixo^s) = ( 1.0d0 - vr(ixo^s)/(x(ixo^s,1) * dvrdr(ixo^s)) ) &
261 * (drstar/x(ixo^s,1))**2.0d0
262
263 select case (cak_1d_opt)
264 case(radstream)
265 fdfac(ixo^s) = 1.0d0
266 case(fdisc, fdisc_cutoff)
267 where (beta_fd(ixo^s) >= 1.0d0)
268 fdfac(ixo^s) = 1.0d0/(1.0d0 + cak_alpha)
269 elsewhere (beta_fd(ixo^s) < -1.0d10)
270 fdfac(ixo^s) = abs(beta_fd(ixo^s))**cak_alpha / (1.0d0 + cak_alpha)
271 elsewhere (abs(beta_fd(ixo^s)) > 1.0d-3)
272 fdfac(ixo^s) = (1.0d0 - (1.0d0 - beta_fd(ixo^s))**(1.0d0 + cak_alpha)) &
273 / (beta_fd(ixo^s)*(1.0d0 + cak_alpha))
274 elsewhere
275 fdfac(ixo^s) = 1.0d0 - 0.5d0*cak_alpha*beta_fd(ixo^s) &
276 * (1.0d0 + 1.0d0/3.0d0 * (1.0d0 - cak_alpha)*beta_fd(ixo^s))
277 endwhere
278 end select
279
280 ! Correct radial line force for finite disc (if applicable)
281 gcak(ixo^s,1) = gcak(ixo^s,1) * fdfac(ixo^s)
282 gcak(ixo^s,2) = 0.0d0
283 gcak(ixo^s,3) = 0.0d0
284
285 ! Fill the nwextra slots for output
286 w(ixo^s,gcak1_) = gcak(ixo^s,1)
287 w(ixo^s,fdf_) = fdfac(ixo^s)
288
289 end subroutine get_cak_force_radial
290
291 !> Vector CAK line force in the Gayley line-ensemble distribution parametrisation
292 subroutine get_cak_force_vector(ixI^L,ixO^L,wCT,w,x,gcak)
295
296 ! Subroutine arguments
297 integer, intent(in) :: ixI^L, ixO^L
298 real(8), intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
299 real(8), intent(inout) :: w(ixI^S,1:nw)
300 real(8), intent(inout) :: gcak(ixO^S,1:3)
301
302 ! Local variables
303 real(8) :: a1, a2, a3, wyray, y, wpray, phiray, wtot, mustar, dvndn
304 real(8) :: costp, costp2, sintp, cospp, sinpp, cott0
305 real(8) :: vr(ixI^S), vt(ixI^S), vp(ixI^S)
306 real(8) :: vrr(ixI^S), vtr(ixI^S), vpr(ixI^S)
307 real(8) :: dvrdr(ixO^S), dvtdr(ixO^S), dvpdr(ixO^S)
308 real(8) :: dvrdt(ixO^S), dvtdt(ixO^S), dvpdt(ixO^S)
309 real(8) :: dvrdp(ixO^S), dvtdp(ixO^S), dvpdp(ixO^S)
310 integer :: ix^D, itray, ipray
311
312 ! Initialisation to have full velocity strain tensor expression at all times
313 vt(ixo^s) = 0.0d0; vtr(ixo^s) = 0.0d0
314 vp(ixo^s) = 0.0d0; vpr(ixo^s) = 0.0d0
315 cott0 = 0.0d0
316 dvrdr(ixo^s) = 0.0d0; dvtdr(ixo^s) = 0.0d0; dvpdr(ixo^s) = 0.0d0
317 dvrdt(ixo^s) = 0.0d0; dvtdt(ixo^s) = 0.0d0; dvpdt(ixo^s) = 0.0d0
318 dvrdp(ixo^s) = 0.0d0; dvtdp(ixo^s) = 0.0d0; dvpdp(ixo^s) = 0.0d0
319
320 ! Populate velocity field(s) depending on dimensions and directions
321 vr(ixi^s) = wct(ixi^s,iw_mom(1)) / wct(ixi^s,iw_rho)
322 vrr(ixi^s) = vr(ixi^s) / x(ixi^s,1)
323
324 {^nooned
325 vt(ixi^s) = wct(ixi^s,iw_mom(2)) / wct(ixi^s,iw_rho)
326 vtr(ixi^s) = vt(ixi^s) / x(ixi^s,1)
327
328 if (ndir > 2) then
329 vp(ixi^s) = wct(ixi^s,iw_mom(3)) / wct(ixi^s,iw_rho)
330 vpr(ixi^s) = vp(ixi^s) / x(ixi^s,1)
331 endif
332 }
333
334 ! Derivatives of velocity field in each coordinate direction (r=1,t=2,p=3)
335 call get_velocity_gradient(ixi^l,ixo^l,vr,x,1,dvrdr)
336
337 {^nooned
338 call get_velocity_gradient(ixi^l,ixo^l,vr,x,2,dvrdt)
339 call get_velocity_gradient(ixi^l,ixo^l,vt,x,1,dvtdr)
340 call get_velocity_gradient(ixi^l,ixo^l,vt,x,2,dvtdt)
341
342 if (ndir > 2) then
343 call get_velocity_gradient(ixi^l,ixo^l,vp,x,1,dvpdr)
344 call get_velocity_gradient(ixi^l,ixo^l,vp,x,2,dvpdt)
345 endif
346 }
347 {^ifthreed
348 call get_velocity_gradient(ixi^l,ixo^l,vr,x,3,dvrdp)
349 call get_velocity_gradient(ixi^l,ixo^l,vt,x,3,dvtdp)
350 call get_velocity_gradient(ixi^l,ixo^l,vp,x,3,dvpdp)
351 }
352
353 ! Get total acceleration from all rays at a certain grid point
354 {do ix^db=ixomin^db,ixomax^db\}
355 ! Loop over the rays; first theta then phi radiation angle
356 ! Get weights from current ray and their position
357 do itray = 1,nthetaray
358 wyray = wy(itray)
359 y = ay(itray)
360
361 do ipray = 1,nphiray
362 wpray = wphi(ipray)
363 phiray = aphi(ipray)
364
365 ! Redistribute the phi rays by a small offset
366 ! if (mod(itp,3) == 1) then
367 ! phip = phip + dphi/3.0d0
368 ! elseif (mod(itp,3) == 2) then
369 ! phip = phip - dphi/3.0d0
370 ! endif
371
372 ! === Geometrical factors ===
373 ! Make y quadrature linear in mu, not mu**2; better for gtheta,gphi
374 ! y -> mu quadrature is preserved; y=0 <=> mu=1; y=1 <=> mu=mustar
375 mustar = sqrt(max(1.0d0 - (drstar/x(ix^d,1))**2.0d0, 0.0d0))
376 costp = 1.0d0 - y*(1.0d0 - mustar)
377 costp2 = costp*costp
378 sintp = sqrt(max(1.0d0 - costp2, 0.0d0))
379 sinpp = sin(phiray)
380 cospp = cos(phiray)
381 {^nooned cott0 = cos(x(ix^d,2))/sin(x(ix^d,2))}
382
383 ! More weight close to star, less farther away
384 wtot = wyray * wpray * (1.0d0 - mustar)
385
386 ! Convenients a la Cranmer & Owocki (1995)
387 a1 = costp
388 a2 = sintp * cospp
389 a3 = sintp * sinpp
390
391 ! Get total velocity gradient for one ray with given (theta', phi')
392 dvndn = a1*a1 * dvrdr(ix^d) + a2*a2 * (dvtdt(ix^d) + vrr(ix^d)) &
393 + a3*a3 * (dvpdp(ix^d) + cott0 * vtr(ix^d) + vrr(ix^d)) &
394 + a1*a2 * (dvtdr(ix^d) + dvrdt(ix^d) - vtr(ix^d)) &
395 + a1*a3 * (dvpdr(ix^d) + dvrdp(ix^d) - vpr(ix^d)) &
396 + a2*a3 * (dvpdt(ix^d) + dvtdp(ix^d) - cott0 * vpr(ix^d))
397
398 ! No multiple resonances in CAK
399 dvndn = abs(dvndn)
400
401 ! Convert gradient back from wind coordinates (r',theta',phi') to
402 ! stellar coordinates (r,theta,phi)
403 gcak(ix^d,1) = gcak(ix^d,1) + (dvndn/wct(ix^d,iw_rho))**cak_alpha * a1 * wtot
404 gcak(ix^d,2) = gcak(ix^d,2) + (dvndn/wct(ix^d,iw_rho))**cak_alpha * a2 * wtot
405 gcak(ix^d,3) = gcak(ix^d,3) + (dvndn/wct(ix^d,iw_rho))**cak_alpha * a3 * wtot
406 enddo
407 enddo
408 {enddo\}
409
410 ! Normalisation for line force
411 ! NOTE: extra 1/pi factor comes from integration in radiation Phi angle
412 gcak(ixo^s,:) = (dke*gayley_qbar)**(1.0d0 - cak_alpha)/(1.0d0 - cak_alpha) &
413 * dlum/(4.0d0*dpi*drstar**2.0d0 * dclight**(1.0d0+cak_alpha)) &
414 * gcak(ixo^s,:)/dpi
415
416 if (fix_vector_force_1d) then
417 gcak(ixo^s,2) = 0.0d0
418 gcak(ixo^s,3) = 0.0d0
419 endif
420
421 ! Fill the nwextra slots for output
422 w(ixo^s,gcak1_) = gcak(ixo^s,1)
423 w(ixo^s,gcak2_) = gcak(ixo^s,2)
424 w(ixo^s,gcak3_) = gcak(ixo^s,3)
425
426 end subroutine get_cak_force_vector
427
428 !> Compute continuum radiation force from Thomson scattering
429 subroutine get_gelectron(ixI^L,ixO^L,w,x,ge)
431
432 integer, intent(in) :: ixI^L, ixO^L
433 real(8), intent(in) :: w(ixI^S,1:nw), x(ixI^S,1:ndim)
434 real(8), intent(out):: ge(ixO^S)
435
436 ge(ixo^s) = dke * dlum/(4.0d0*dpi * dclight * x(ixo^s,1)**2.0d0)
437
438 end subroutine get_gelectron
439
440 !> Check time step for total radiation contribution
441 subroutine cak_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
443
444 integer, intent(in) :: ixI^L, ixO^L
445 real(8), intent(in) :: dx^D, x(ixI^S,1:ndim)
446 real(8), intent(in) :: w(ixI^S,1:nw)
447 real(8), intent(inout) :: dtnew
448
449 ! Local variables
450 real(8) :: ge(ixO^S), max_gr, dt_cak
451
452 call get_gelectron(ixi^l,ixo^l,w,x,ge)
453
454 dtnew = bigdouble
455
456 ! Get dt from line force that is saved in the w-array in nwextra slot
457 max_gr = max( maxval(abs(ge(ixo^s) + w(ixo^s,gcak1_))), epsilon(1.0d0) )
458 dt_cak = minval( sqrt(block%dx(ixo^s,1)/max_gr) )
459 dtnew = min(dtnew, courantpar*dt_cak)
460
461 {^nooned
462 if (cak_vector_force) then
463 max_gr = max( maxval(abs(w(ixo^s,gcak2_))), epsilon(1.0d0) )
464 dt_cak = minval( sqrt(block%dx(ixo^s,1) * block%dx(ixo^s,2)/max_gr) )
465 dtnew = min(dtnew, courantpar*dt_cak)
466
467 {^ifthreed
468 max_gr = max( maxval(abs(w(ixo^s,gcak3_))), epsilon(1.0d0) )
469 dt_cak = minval( sqrt(block%dx(ixo^s,1) * sin(block%dx(ixo^s,3))/max_gr) )
470 dtnew = min(dtnew, courantpar*dt_cak)
471 }
472 endif
473 }
474
475 end subroutine cak_get_dt
476
477 !> Compute velocity gradient in direction 'idir' on a non-uniform grid
478 subroutine get_velocity_gradient(ixI^L,ixO^L,vfield,x,idir,grad_vn)
480
481 integer, intent(in) :: ixI^L, ixO^L, idir
482 real(8), intent(in) :: vfield(ixI^S), x(ixI^S,1:ndim)
483 real(8), intent(out) :: grad_vn(ixO^S)
484
485 ! Local variables
486 real(8) :: forw(ixO^S), backw(ixO^S), cent(ixO^S)
487 integer :: jrx^L, hrx^L{^NOONED,jtx^L, htx^L}{^IFTHREED,jpx^L, hpx^L}
488
489 ! Index +1 (j) and index -1 (h) in radial direction; kr(dir,dim)=1, dir=dim
490 jrx^l=ixo^l+kr(1,^d);
491 hrx^l=ixo^l-kr(1,^d);
492
493 {^nooned
494 ! Index +1 (j) and index -1 (h) in polar direction
495 jtx^l=ixo^l+kr(2,^d);
496 htx^l=ixo^l-kr(2,^d);
497 }
498
499 {^ifthreed
500 ! Index +1 (j) and index -1 (h) in azimuthal direction
501 jpx^l=ixo^l+kr(3,^d);
502 hpx^l=ixo^l-kr(3,^d);
503 }
504
505 ! grad(v.n) on non-uniform grid according to Sundqvist & Veronis (1970)
506 select case (idir)
507 case(1) ! Radial forward, backward, and central derivatives
508 forw(ixo^s) = (x(ixo^s,1) - x(hrx^s,1)) * vfield(jrx^s) &
509 / ((x(jrx^s,1) - x(ixo^s,1)) * (x(jrx^s,1) - x(hrx^s,1)))
510
511 backw(ixo^s) = -(x(jrx^s,1) - x(ixo^s,1)) * vfield(hrx^s) &
512 / ((x(ixo^s,1) - x(hrx^s,1)) * (x(jrx^s,1) - x(hrx^s,1)))
513
514 cent(ixo^s) = (x(jrx^s,1) + x(hrx^s,1) - 2.0d0*x(ixo^s,1)) * vfield(ixo^s) &
515 / ((x(ixo^s,1) - x(hrx^s,1)) * (x(jrx^s,1) - x(ixo^s,1)))
516 {^nooned
517 case(2) ! Polar forward, backward, and central derivatives
518 forw(ixo^s) = (x(ixo^s,2) - x(htx^s,2)) * vfield(jtx^s) &
519 / (x(ixo^s,1) * (x(jtx^s,2) - x(ixo^s,2)) * (x(jtx^s,2) - x(htx^s,2)))
520
521 backw(ixo^s) = -(x(jtx^s,2) - x(ixo^s,2)) * vfield(htx^s) &
522 / ( x(ixo^s,1) * (x(ixo^s,2) - x(htx^s,2)) * (x(jtx^s,2) - x(htx^s,2)))
523
524 cent(ixo^s) = (x(jtx^s,2) + x(htx^s,2) - 2.0d0*x(ixo^s,2)) * vfield(ixo^s) &
525 / ( x(ixo^s,1) * (x(ixo^s,2) - x(htx^s,2)) * (x(jtx^s,2) - x(ixo^s,2)))
526 }
527 {^ifthreed
528 case(3) ! Azimuthal forward, backward, and central derivatives
529 forw(ixo^s) = (x(ixo^s,3) - x(hpx^s,3)) * vfield(jpx^s) &
530 / ( x(ixo^s,1)*sin(x(ixo^s,2)) * (x(jpx^s,3) - x(ixo^s,3)) * (x(jpx^s,3) - x(hpx^s,3)))
531
532 backw(ixo^s) = -(x(jpx^s,3) - x(ixo^s,3)) * vfield(hpx^s) &
533 / ( x(ixo^s,1)*sin(x(ixo^s,2)) * (x(ixo^s,3) - x(hpx^s,3)) * (x(jpx^s,3) - x(hpx^s,3)))
534
535 cent(ixo^s) = (x(jpx^s,3) + x(hpx^s,3) - 2.0d0*x(ixo^s,3)) * vfield(ixo^s) &
536 / ( x(ixo^s,1)*sin(x(ixo^s,2)) * (x(ixo^s,3) - x(hpx^s,3)) * (x(jpx^s,3) - x(ixo^s,3)))
537 }
538 end select
539
540 ! Total gradient for given velocity field
541 grad_vn(ixo^s) = backw(ixo^s) + cent(ixo^s) + forw(ixo^s)
542
543 end subroutine get_velocity_gradient
544
545 !> Initialise (theta',phi') radiation angles coming from stellar disc
546 subroutine rays_init(ntheta_point,nphi_point)
548
549 ! Subroutine arguments
550 integer, intent(in) :: ntheta_point, nphi_point
551
552 ! Local variables
553 real(8) :: ymin, ymax, phipmin, phipmax, adum
554 integer :: ii
555
556 ! Minimum and maximum range of theta and phi rays
557 ! NOTE: theta points are cast into y-space
558 ymin = 0.0d0
559 ymax = 1.0d0
560 phipmin = -dpi !0.0d0
561 phipmax = dpi !2.0d0*dpi
562 ! dphi = (phipmax - phipmin) / nphi_point
563
564 if (mype == 0) then
565 allocate(ay(ntheta_point))
566 allocate(wy(ntheta_point))
567 allocate(aphi(nphi_point))
568 allocate(wphi(nphi_point))
569
570 ! theta and phi ray positions and weights: Gauss-Legendre
571 call gauss_legendre_quadrature(ymin,ymax,ntheta_point,ay,wy)
572 call gauss_legendre_quadrature(phipmin,phipmax,nphi_point,aphi,wphi)
573
574 ! theta rays and weights: uniform
575 ! dth = 1.0d0 / nthetap
576 ! adum = ymin + 0.5d0*dth
577 ! do ip = 1,nthetap
578 ! ay(ip) = adum
579 ! wy(ip) = 1.0d0/nthetap
580 ! adum = adum + dth
581 ! !print*,'phipoints'
582 ! !print*,ip,aphi(ip),wphi(ip),dphi
583 ! enddo
584
585 ! phi ray position and weights: uniform
586 ! adum = phipmin + 0.5d0*dphi
587 ! do ii = 1,nphi_point
588 ! aphi(ii) = adum
589 ! wphi(ii) = 1.0d0/nphi_point
590 ! adum = adum + dphi
591 ! enddo
592
593 print*, '==========================='
594 print*, ' Radiation ray setup '
595 print*, '==========================='
596 print*, 'Theta ray points + weights '
597 do ii = 1,ntheta_point
598 print*,ii,ay(ii),wy(ii)
599 enddo
600 print*
601 print*, 'Phi ray points + weights '
602 do ii = 1,nphi_point
603 print*,ii,aphi(ii),wphi(ii)
604 enddo
605 print*
606 endif
607
608 call mpi_barrier(icomm,ierrmpi)
609
610 !===========================
611 ! Broadcast what mype=0 read
612 !===========================
613 if (npe > 1) then
614 call mpi_bcast(ntheta_point,1,mpi_integer,0,icomm,ierrmpi)
615 call mpi_bcast(nphi_point,1,mpi_integer,0,icomm,ierrmpi)
616
617 if (mype /= 0) then
618 allocate(ay(ntheta_point))
619 allocate(wy(ntheta_point))
620 allocate(aphi(nphi_point))
621 allocate(wphi(nphi_point))
622 endif
623
624 call mpi_bcast(ay,ntheta_point,mpi_double_precision,0,icomm,ierrmpi)
625 call mpi_bcast(wy,ntheta_point,mpi_double_precision,0,icomm,ierrmpi)
626 call mpi_bcast(aphi,nphi_point,mpi_double_precision,0,icomm,ierrmpi)
627 call mpi_bcast(wphi,nphi_point,mpi_double_precision,0,icomm,ierrmpi)
628 endif
629
630 end subroutine rays_init
631
632 !> Fast Gauss-Legendre N-point quadrature algorithm by G. Rybicki
633 subroutine gauss_legendre_quadrature(xlow,xhi,n,x,w)
634 ! Given the lower and upper limits of integration xlow and xhi, and given n,
635 ! this routine returns arrays x and w of length n, containing the abscissas
636 ! and weights of the Gauss-Legendre N-point quadrature
638
639 ! Subroutine arguments
640 real(8), intent(in) :: xlow, xhi
641 integer, intent(in) :: n
642 real(8), intent(out) :: x(n), w(n)
643
644 ! Local variables
645 real(8) :: p1, p2, p3, pp, xl, xm, z, z1
646 real(8), parameter :: error=3.0d-14
647 integer :: i, j, m
648
649 m = (n + 1)/2
650 xm = 0.5d0*(xhi + xlow)
651 xl = 0.5d0*(xhi - xlow)
652
653 do i = 1,m
654 z = cos( dpi * (i - 0.25d0)/(n + 0.5d0) )
655 z1 = 2.0d0 * z
656
657 do while (abs(z1 - z) > error)
658 p1 = 1.0d0
659 p2 = 0.0d0
660
661 do j = 1,n
662 p3 = p2
663 p2 = p1
664 p1 = ( (2.0d0*j - 1.0d0)*z*p2 - (j - 1.0d0)*p3 )/j
665 enddo
666
667 pp = n*(z*p1 - p2) / (z*z - 1.0d0)
668 z1 = z
669 z = z1 - p1/pp
670 enddo
671
672 x(i) = xm - xl*z
673 x(n+1-i) = xm + xl*z
674 w(i) = 2.0d0*xl / ((1.0d0 - z*z) * pp*pp)
675 w(n+1-i) = w(i)
676 enddo
677
678 end subroutine gauss_legendre_quadrature
679
680end module mod_cak_force
Module to include CAK radiation line force in (magneto)hydrodynamic models Computes both the force fr...
real(8), public gayley_qbar
subroutine cak_get_dt(w, ixil, ixol, dtnew, dxd, x)
Check time step for total radiation contribution.
real(8), public gayley_q0
logical cak_split
To treat source term in split or unsplit (default) fashion.
subroutine cak_init(phys_gamma)
Initialize the module.
subroutine cak_params_read(files)
Public method.
subroutine gauss_legendre_quadrature(xlow, xhi, n, x, w)
Fast Gauss-Legendre N-point quadrature algorithm by G. Rybicki.
real(8), public cak_alpha
Line-ensemble parameters in the Gayley (1995) formalism.
subroutine cak_add_source(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
subroutine, public set_cak_force_norm(rstar, twind)
Compute some (unitless) variables for CAK force normalisation.
logical fix_vector_force_1d
To activate the pure radial vector CAK line force computation.
integer gcak1_
Extra slots to store quantities in w-array.
logical cak_vector_force
To activate the vector CAK line force computation.
subroutine get_velocity_gradient(ixil, ixol, vfield, x, idir, grad_vn)
Compute velocity gradient in direction 'idir' on a non-uniform grid.
integer cak_1d_opt
Switch to choose between the 1-D CAK line force options.
subroutine get_cak_force_radial(ixil, ixol, wct, w, x, gcak)
1-D CAK line force in the Gayley line-ensemble distribution parametrisation
subroutine get_gelectron(ixil, ixol, w, x, ge)
Compute continuum radiation force from Thomson scattering.
subroutine get_cak_force_vector(ixil, ixol, wct, w, x, gcak)
Vector CAK line force in the Gayley line-ensemble distribution parametrisation.
logical cak_1d_force
To activate the original CAK 1-D line force computation.
integer nthetaray
Amount of rays in radiation polar and radiation azimuthal direction.
subroutine rays_init(ntheta_point, nphi_point)
Initialise (theta',phi') radiation angles coming from stellar disc.
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for physical and numeric constants.
double precision, parameter dpi
Pi.
double precision, parameter const_kappae
double precision, parameter const_c
double precision, parameter const_sigma
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, dimension(3, 3) kr
Kronecker delta tensor.
double precision unit_length
Physical scaling factor for length.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
double precision courantpar
The Courant (CFL) number used for the simulation.
integer ierrmpi
A global MPI error return code.
integer npe
The number of MPI tasks.
double precision unit_velocity
Physical scaling factor for velocity.
double precision unit_temperature
Physical scaling factor for temperature.
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition mod_physics.t:4
procedure(sub_get_pthermal), pointer phys_get_pthermal
Definition mod_physics.t:73
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition mod_physics.t:50
Module with all the methods that users can customize in AMRVAC.