MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
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  !> Switch to choose between the 1-D CAK line force options
35  integer :: cak_1d_opt
36 
37  ! Avoid magic numbers in code for 1-D CAK line force option
38  integer, parameter, private :: radstream=0, fdisc=1, fdisc_cutoff=2
39 
40  !> To treat source term in split or unsplit (default) fashion
41  logical :: cak_split=.false.
42 
43  !> To activate the original CAK 1-D line force computation
44  logical :: cak_1d_force=.false.
45 
46  !> To activate the vector CAK line force computation
47  logical :: cak_vector_force=.false.
48 
49  !> To activate the pure radial vector CAK line force computation
50  logical :: fix_vector_force_1d=.false.
51 
52  !> Amount of rays in radiation polar and radiation azimuthal direction
53  integer :: nthetaray, nphiray
54 
55  !> Ray positions + weights for impact parameter and azimuthal radiation angle
56  real(8), allocatable, private :: ay(:), wy(:), aphi(:), wphi(:)
57 
58  !> The adiabatic index
59  real(8), private :: cak_gamma
60 
61  !> Variables needed to compute force normalisation fnorm in initialisation
62  real(8), private :: lum, dlum, drstar, dke, dclight
63 
64  !> To enforce a floor temperature when doing adiabatic (M)HD
65  real(8), private :: tfloor
66 
67  !> Extra slots to store quantities in w-array
68  integer :: gcak1_, gcak2_, gcak3_, fdf_
69 
70  !> Public method
71  public :: set_cak_force_norm
72 
73 contains
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  integer :: idir
171  real(8) :: gl(ixO^S,1:3), ge(ixO^S), ptherm(ixI^S), pmin(ixI^S)
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)
294  use mod_usr_methods
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  integer :: ix^D, itray, ipray
304  real(8) :: a1, a2, a3, wyray, y, wpray, phiray, wtot, mustar, dvndn
305  real(8) :: costp, costp2, sintp, cospp, sinpp, cott0
306  real(8) :: vr(ixI^S), vt(ixI^S), vp(ixI^S)
307  real(8) :: vrr(ixI^S), vtr(ixI^S), vpr(ixI^S)
308  real(8) :: dvrdr(ixO^S), dvtdr(ixO^S), dvpdr(ixO^S)
309  real(8) :: dvrdt(ixO^S), dvtdt(ixO^S), dvpdt(ixO^S)
310  real(8) :: dvrdp(ixO^S), dvtdp(ixO^S), dvpdp(ixO^S)
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  integer :: i, j, m
646  real(8) :: p1, p2, p3, pp, xl, xm, z, z1
647  real(8), parameter :: error=3.0d-14
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 
680 end module mod_cak_force
Module to include CAK radiation line force in (magneto)hydrodynamic models Computes both the force fr...
Definition: mod_cak_force.t:27
subroutine get_gelectron(ixIL, ixOL, w, x, ge)
Compute continuum radiation force from Thomson scattering.
real(8), public gayley_qbar
Definition: mod_cak_force.t:32
subroutine get_velocity_gradient(ixIL, ixOL, vfield, x, idir, grad_vn)
Compute velocity gradient in direction 'idir' on a non-uniform grid.
real(8), public gayley_q0
Definition: mod_cak_force.t:32
subroutine cak_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
Check time step for total radiation contribution.
logical cak_split
To treat source term in split or unsplit (default) fashion.
Definition: mod_cak_force.t:41
subroutine cak_init(phys_gamma)
Initialize the module.
Definition: mod_cak_force.t:98
subroutine cak_params_read(files)
Public method.
Definition: mod_cak_force.t:77
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.
Definition: mod_cak_force.t:32
subroutine, public set_cak_force_norm(rstar, twind)
Compute some (unitless) variables for CAK force normalisation.
subroutine get_cak_force_radial(ixIL, ixOL, wCT, w, x, gcak)
1-D CAK line force in the Gayley line-ensemble distribution parametrisation
logical fix_vector_force_1d
To activate the pure radial vector CAK line force computation.
Definition: mod_cak_force.t:50
integer gcak1_
Extra slots to store quantities in w-array.
Definition: mod_cak_force.t:68
logical cak_vector_force
To activate the vector CAK line force computation.
Definition: mod_cak_force.t:47
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
integer gcak2_
Definition: mod_cak_force.t:68
integer cak_1d_opt
Switch to choose between the 1-D CAK line force options.
Definition: mod_cak_force.t:35
integer gcak3_
Definition: mod_cak_force.t:68
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.
Definition: mod_cak_force.t:44
integer nthetaray
Amount of rays in radiation polar and radiation azimuthal direction.
Definition: mod_cak_force.t:53
subroutine rays_init(ntheta_point, nphi_point)
Initialise (theta',phi') radiation angles coming from stellar disc.
integer nphiray
Definition: mod_cak_force.t:53
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
Module for physical and numeric constants.
Definition: mod_constants.t:2
double precision, parameter dpi
Pi.
Definition: mod_constants.t:25
double precision, parameter const_kappae
Definition: mod_constants.t:62
double precision, parameter const_c
Definition: mod_constants.t:48
double precision, parameter const_sigma
Definition: mod_constants.t:59
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.
integer, 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:77
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition: mod_physics.t:16
Module with all the methods that users can customize in AMRVAC.