DISPATCH
non_ideal_mod.f90
1 !===============================================================================
2 !> non-ideal MHD module
3 !===============================================================================
4 MODULE non_ideal_mod
5  USE io_mod
6  USE trace_mod
7  USE io_unit_mod
8  USE scalar_mod
9  USE stagger_mod
10  implicit none
11  private
12  type, public:: non_ideal_t
13  real:: eta_ohm, gamma_ad, rho_ions
14  logical:: mhd_ohm, mhd_ad, is_used
15  contains
16  procedure:: init
17  procedure:: update
18  end type
19  type(non_ideal_t), public:: non_ideal
20 CONTAINS
21 
22 !===============================================================================
23 !> Read parameters
24 !===============================================================================
25 SUBROUTINE init (self)
26  class(non_ideal_t):: self
27  integer:: iostat
28  real:: gamma_ad=75., rho_ions=1.0, eta_ohm=0.0
29  logical:: mhd_ad=.false., mhd_ohm=.false.
30  namelist /non_ideal_params/mhd_ad, mhd_ohm, gamma_ad, eta_ohm, rho_ions
31  !---------------------------------------------------------------------------
32  rewind(io_unit%input)
33  read (io_unit%input, non_ideal_params,iostat=iostat)
34  if (io%master) write (*, non_ideal_params)
35  self%mhd_Ohm = mhd_ohm
36  self%eta_Ohm = eta_ohm
37  self%mhd_AD = mhd_ad
38  self%gamma_AD = gamma_ad
39  self%rho_ions = rho_ions
40  self%is_used = mhd_ohm .or. mhd_ad
41 END SUBROUTINE init
42 
43 !===============================================================================
44 !> Update heating Q, electric field E, and Courant velocity u_max
45 !===============================================================================
46 SUBROUTINE update (self, gn, ds, d, Jx, Jy, Jz, Bx, By, Bz, Q, E, u_max)
47  class(non_ideal_t):: self
48  integer:: gn(3)
49  real(8):: ds
50  real:: u_max, gamma_ad
51  logical::mhd_ad, mhd_ohm
52  real, dimension(:,:,:,:):: e
53  real, dimension(:,:,:):: d, jx, jy, jz, bx, by, bz, q
54  real, dimension(:,:,:), allocatable:: ud_x, ud_y, ud_z
55  real, dimension(:,:,:), allocatable:: ambi
56  integer, save:: itimer=0
57  !---------------------------------------------------------------------------
58  ! Below is the actual ambipolar drift parameter as given in equation 20 in
59  ! Padoan, Zweibel and Nordlund (2000).
60  ! To only write ambi, the definition in Padoan et al. is already multiplied
61  ! by density**(-3/2) as it appears in the calculation of the drift velocity.
62  ! The term d**1.5 is a cosnequence of the single fluid approximation and the
63  ! relation n_i = K(n_n/10^5cm^-3)^k + K' (n_n/10^3 cm^-3)^-2
64  ! K = 3*10^-3cm^-3, k=0.5 and K'=4.64*10^-4cm^-3.
65  ! based on assuming that dust grains do NOT affect ionization balance significantly
66  ! => second term falls off quickly => n_i~n_n^0.5
67  ! In order to compute this parameter correctly, we need to know TEMPERATURE,
68  ! RESOLUTION (not directly applicable here because we use AMR), BOX LENGTH,
69  ! AVERAGE DENSITY, COSMIC-RAY IONISATION (K_i parameter:
70  ! K_i = n_i*n_H**(-0.5) = x_i*n_H**(-0.5) ) and AVERAGE VOLUMETRIC FLOW RATE av_sigv:
71  !---------------------------------------------------------------------------
72  ! ambi = 0.28 * (0.1*T)**0.5 * (N/128d0) * (L/10d0)**(-1) * (av_n/200d0)**(-0.5)
73  ! * (K_i/1e-5)**(-1) * (av_sigv/2e-9)**(-1) * d**(-1.5)
74  !
75  ! Alternative definition in Duffin and Pudritz (2008) assumes constant parameter
76  ! for gamma_AD
77  ! gamma_AD = av_sigv / (m_i + m_n) = 3.28e13 g^-1cm^3s^-1
78  ! ambi = 1.4 / (gamma_AD*d**1.5)
79  ! rho_i = 1.
80  !---------------------------------------------------------------------------
81  if (.not.self%is_used) return
82  call trace%begin ('non_ideal_t%update', itimer=itimer)
83  if (self%mhd_AD) then
84  call allocate_scalars_a (gn, ambi, ud_x, ud_y, ud_z)
85  ambi = 1./(self%gamma_AD*d*self%rho_ions)
86  ud_x = xdn(ambi) * (zup(jy)*xdn(zup(bz)) - yup(jz)*xdn(yup(by)))
87  ud_y = ydn(ambi) * (xup(jz)*ydn(xup(bx)) - zup(jx)*ydn(zup(bz)))
88  ud_z = zdn(ambi) * (yup(jx)*zdn(yup(by)) - xup(jy)*zdn(xup(bx)))
89  e(:,:,:,1) = e(:,:,:,1) - (ud_y*bz - ud_z*by)
90  e(:,:,:,2) = e(:,:,:,2) - (ud_z*bx - ud_x*bz)
91  e(:,:,:,3) = e(:,:,:,3) - (ud_x*by - ud_y*bx)
92  q = q + self%rho_ions*d*self%gamma_AD*(xup(ud_x)**2+yup(ud_y)**2+zup(ud_z)**2)
93  call deallocate_scalars_a (ambi, ud_x, ud_y, ud_z)
94  u_max = max(u_max,maxval((xup(bx)**2+yup(by)**2+zup(bz)**2)/(self%gamma_AD*self%rho_ions*d*ds)))
95  end if
96  if (self%mhd_Ohm) then
97  e(:,:,:,1) = e(:,:,:,1) + self%eta_Ohm*jx
98  e(:,:,:,2) = e(:,:,:,2) + self%eta_Ohm*jy
99  e(:,:,:,3) = e(:,:,:,3) + self%eta_Ohm*jz
100  q = q + yup(zup(self%eta_Ohm*jx**2)) + zup(xup(self%eta_Ohm*jy**2)) + xup(yup(self%eta_Ohm*jz**2))
101  u_max = max(u_max,self%eta_Ohm/ds)
102  end if
103  call trace%end (itimer)
104 END SUBROUTINE update
105 
106 END MODULE non_ideal_mod
non-ideal MHD module
Definition: io_mod.f90:4
6th order stagger operators, with self-test procedure
Definition: stagger_2nd.f90:4