Streamer fluid modeling - An overview of ARCoS  1.0
 All Classes Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
cdr_advect_diffu_vec.f90
Go to the documentation of this file.
1 
37 
38 ! subroutine cdr_advect_vec
39 ! subroutine mat_efield_r
40 ! subroutine mat_efield_z
41 ! subroutine mat_cdr_f_ad_r
42 ! subroutine mat_cdr_f_ad_z
43 ! subroutine make_vecs_er
44 
45  subroutine cdr_advect_vec(mass,charge,dens_array,d_dens_array, &
46  er_array,ez_array,diffusion_coeff,dr,dz, &
47  sprite_module,r0,r1,z0,z1)
48  implicit none
49 ! ..
50 ! .. Parameters ..
51  double precision,parameter :: one = 1.0d0
52 ! ..
53 ! .. Dummy Variables ..
54  logical,intent(in) :: sprite_module
55  integer,intent(in) :: r0,r1,z0,z1
56  double precision,intent(in) :: charge,dr,dz,mass,diffusion_coeff
57 ! ..
58 ! .. Dummy Arrays ..
59  double precision,dimension(r0-2:r1+1,z0-2:z1+1),intent(in) :: &
60  dens_array,er_array,ez_array
61  double precision,dimension(r0-2:r1+1,z0-2:z1+1),intent(out) :: &
62  d_dens_array
63 ! ..
64 ! .. Local Variables ..
65  integer :: ir,ishift
66  double precision :: cm0,cm,gdz_inv,d_iso0,d_iso
67 ! ..
68 ! .. Local arrays ..
69  double precision,dimension(r0:r1-1) :: er_min
70  double precision,dimension(r0+1:r1) :: er
71  double precision,dimension(z1-z0) :: d_iso_r
72  double precision,dimension(z1-z0+1) :: d_iso_z
73  double precision,dimension(r1-r0+1,z1-z0 ) :: mat_r
74  double precision,dimension(r1-r0 ,z1-z0+1) :: mat_z
75  double precision,dimension(r1-r0+1,z1-z0 ) :: efield_r
76  double precision,dimension(r0:r1-1,z0-1:z1-1) :: efield_z
77 ! ..
78 ! .. External subroutines ..
79  external :: mat_efield_r,mat_efield_z,mat_cdr_f_ad_r,mat_cdr_f_ad_z
80 ! ..
81 
82  if (mass <= 0) return
83 
84  ! <ISOTROPIC DIFFUSSION>
85  d_iso0 = diffusion_coeff / mass
86  ! </ISOTROPIC DIFFUSION>
87 
88  gdz_inv = one / dz
89 
90  cm0 = charge / mass
91 
92  if (.not.sprite_module) then
93  cm = cm0
94  d_iso = d_iso0
95  end if
96 
97 !-----------------------------------------------
98 ! Increase the values r0, r1, z0 and z1 with 3:
99 ! because of 2 buffer cells at the left side and to compensate the C 0
100 !-----------------------------------------------
101  ishift = 3
102 
103  ! Compute the electric field in radial direction
104  call mat_efield_r(er_array(r0-1:r1-1,z0:z1-1), &
105  efield_r,d_iso_r,cm0,d_iso0, &
106  r0,r1,dz,z0,z1,ishift,dr,sprite_module)
107  ! Calculate the radial advection from left and right
108  call mat_cdr_f_ad_r(dens_array(r0-2:r1+1,z0:z1-1), &
109  efield_r(1:r1-r0+1,1:z1-z0),d_iso_r, &
110  mat_r,r0,r1,z0,z1)
111  ! Add the radial advection
112  call make_vecs_er(er_min,er,r0,r1,dr)
113  do ir = r0,r1-1
114  d_dens_array(ir,z0:z1-1) = d_dens_array(ir,z0:z1-1) + &
115  er_min(ir) * mat_r(ir-r0+1,1:z1-z0) - &
116  er(ir+1) * mat_r(ir-r0+2,1:z1-z0)
117  end do
118 
119  ! Compute the electric field in axial direction
120  call mat_efield_z(ez_array(r0 :r1-1,z0-1:z1-1), &
121  efield_z(r0:r1-1,z0-1:z1-1),d_iso_z,cm0,d_iso0, &
122  r0,r1,dz,z0,z1,ishift,dz,sprite_module)
123  ! Calculate the axial advection from the north and south
124  call mat_cdr_f_ad_z(dens_array(r0 :r1-1,z0-2:z1+1), &
125  efield_z(r0:r1-1,z0-1:z1-1),d_iso_z, &
126  mat_z,r0,r1,z0,z1)
127  ! Add the axial advection
128  d_dens_array(r0:r1-1,z0:z1-1) = d_dens_array(r0:r1-1,z0:z1-1) + &
129  gdz_inv * (mat_z(1:,1:z1-z0) - &
130  mat_z(1:,2:z1-z0+1))
131 
132  end subroutine cdr_advect_vec
133 
134  subroutine mat_efield_r(er_array,efield,d_iso,cm0,d_iso0, &
135  r0,r1,dz,z0,z1,z_shift,dx,sprite_module)
136 ! ..
137  implicit none
138 ! ..
139 ! .. Parameters ..
140  double precision,parameter :: half = 0.5d0,one = 1.0d0
141 ! ..
142 ! .. Dummy Variables ..
143  logical,intent(in) :: sprite_module
144  integer,intent(in) :: r0,r1,z0,z1,z_shift
145  double precision,intent(in) :: cm0,d_iso0,dz,dx
146 ! ..
147 ! .. Dummy Arrays ..
148  double precision,dimension(r1-r0+1,z1-z0),intent(in) :: er_array
149  double precision,dimension(r1-r0+1,z1-z0),intent(out) :: efield
150  double precision,dimension(z1-z0),intent(out) :: d_iso
151 ! ..
152 ! .. Local Variables ..
153  integer :: iz
154 ! ..
155 ! .. Local arrays ..
156  double precision,dimension(z1-z0) :: back_dens_inv,cm
157 !-----------------------------------------------
158 
159  ! If the densities are varying...
160  if (sprite_module) then
161  do iz = z0,z1
162  ! spr_density_at(iz) = (iz-z_shift + half) * dz
163  ! back_dens_inv = ONE / spr_density_at
164  back_dens_inv(iz-z0+1) = one / ((iz-z_shift + half) * dz)
165  end do
166  cm = cm0 * back_dens_inv
167  d_iso = d_iso0 * back_dens_inv
168  else
169  back_dens_inv = one
170  cm = cm0
171  d_iso = d_iso0
172  end if
173 
174  d_iso = d_iso / dx
175 
176  ! efield = cm * er_array
177  do iz = 1,z1-z0
178  efield(:,iz) = cm(iz) * er_array(:,iz)
179  end do
180 
181  end subroutine mat_efield_r
182 
183  subroutine mat_efield_z (ez_array,efield,d_iso,cm0,d_iso0, &
184  r0,r1,dz,z0,z1,z_shift,dx,sprite_module)
185 ! ..
186  implicit none
187 ! ..
188 ! .. Parameters ..
189  double precision,parameter :: half = 0.5d0,one = 1.0d0
190 ! ..
191 ! .. Dummy Variables ..
192  logical,intent(in) :: sprite_module
193  integer,intent(in) :: r0,r1,z0,z1,z_shift
194  double precision,intent(in) :: cm0,d_iso0,dz,dx
195 ! ..
196 ! .. Dummy Arrays ..
197  double precision,dimension(r0:r1-1,z0-1:z1-1),intent(in) :: ez_array
198  double precision,dimension(r0:r1-1,z0-1:z1-1),intent(out) :: efield
199  double precision,dimension(z0-1:z1-1),intent(out) :: d_iso
200 ! ..
201 ! .. Local Variables ..
202  integer :: iz
203 ! ..
204 ! .. Local arrays ..
205  double precision,dimension(z0-1:z1-1) :: back_dens_inv,cm
206 ! ..
207 !-----------------------------------------------
208 
209  ! If the densities are varying...
210  if (sprite_module) then
211  do iz = z0-1,z1+1
212  ! spr_density_at(iz) = (iz-z_shift + half) * dz
213  ! back_dens_inv = ONE / spr_density_at
214  back_dens_inv(iz) = one / ((iz-z_shift - half) * dz)
215  end do
216  cm = cm0 * back_dens_inv
217  d_iso = d_iso0 * back_dens_inv
218  else
219  back_dens_inv = one
220  cm = cm0
221  d_iso = d_iso0
222  end if
223 
224  d_iso = d_iso / dx
225 
226  ! efield = cm * ez_array
227  do iz = z0-1,z1-1
228  efield(:,iz) = cm(iz) * ez_array(:,iz)
229  end do
230 
231  end subroutine mat_efield_z
232 
233  subroutine mat_cdr_f_ad_r(data_array,efield,d_iso,mat_r,r0,r1,z0,z1)
234 ! ..
235  implicit none
236 ! ..
237 ! .. Parameters ..
238  double precision,parameter :: zero = 0.0d0,half = 0.5d0, &
239  one = 1.0d0,three = 3.0d0, &
240  verysmall = 1.0e-20
241 ! ..
242 ! .. Dummy Variables ..
243  integer,intent(in) :: r0,r1,z0,z1
244 ! ..
245 ! .. Dummy Arrays ..
246  double precision,dimension(z1-z0),intent(in) :: d_iso
247  double precision,dimension(r1-r0+1,z1-z0),intent(in) :: efield
248  double precision,dimension(r1-r0+4,z1-z0),intent(in) :: data_array
249  double precision,dimension(r1-r0+1,z1-z0),intent(out) :: mat_r
250 ! ..
251 ! .. Local Variables ..
252  integer :: iz,length
253 ! ..
254 ! .. Local arrays ..
255  double precision,dimension(r1-r0+1,z1-z0) :: aux,psi_p
256  double precision,dimension(r1-r0+3,z1-z0) :: mat_diff
257 ! ..
258 !-----------------------------------------------
259  length = r1 - r0 + 1
260  mat_diff(1:length+2,:) = data_array(2:length+3,:) - &
261  data_array(1:length+2,:)
262 
263  where (abs(mat_diff(2:length+1,:)) > verysmall)
264  where (efield>=zero)
265  ! p_min =data_array(ir-k1,iz-k2)
266  ! p =data_array(ir,iz)
267  ! p_plus =data_array(ir+k1,iz+k2)
268 
269  psi_p = mat_diff(1:length,:) / mat_diff(2:length+1,:)
270  elsewhere (efield<zero)
271  ! p =data_array(ir,iz)
272  ! p_plus =data_array(ir+ k1,iz+ k2)
273  ! p_2plus =data_array(ir+2*k1,iz+2*k2)
274 
275  psi_p = mat_diff(3:length+2,:) / mat_diff(2:length+1,:)
276  end where
277  where (psi_p <= zero)
278  ! psi_p = psi_MN (psi_p)
279  psi_p = zero
280  elsewhere (psi_p >= 4)
281  psi_p = one
282  elsewhere (psi_p >= 0.4)
283  psi_p = (one + half * psi_p) / three
284  end where
285  elsewhere ! abs(mat_diff(2:length+1,:)) < verysmall
286  psi_p = -one
287  end where
288 
289  ! sigmasigma is (sigma_{i+1,j} - sigma_{i,j})*/
290  ! sigmasigma = p - p_plus
291  aux = psi_p * mat_diff(2:length+1,:)
292  ! aux = p + psi_p * sigmasigma
293  where (efield>=zero)
294  aux = data_array(2:length+1,:) + aux
295  elsewhere
296  aux = data_array(3:length+2,:) - aux
297  end where
298 
299  ! aux = efield * (p + psi_p * sigmasigma)
300  aux = aux * efield
301  do iz = 1,z1-z0
302  ! mat_r = d_iso * sigmasigma
303  mat_r(1:length,iz) = d_iso(iz) * mat_diff(2:length+1,iz)
304  end do
305  ! mat_r = efield * (p + psi_p * sigmasigma) + d_iso * sigmasigma
306  mat_r = aux - mat_r
307 
308  end subroutine mat_cdr_f_ad_r
309 
310  subroutine mat_cdr_f_ad_z(data_array,efield,d_iso,mat_z,r0,r1,z0,z1)
311 ! ..
312  implicit none
313 ! ..
314 ! .. Parameters ..
315  double precision,parameter :: zero = 0.0d0,half = 0.5d0, &
316  one = 1.0d0,three = 3.0d0, &
317  verysmall = 1.0e-20
318 ! ..
319 ! .. Dummy Variables ..
320  integer,intent(in) :: r0,r1,z0,z1
321 ! ..
322 ! .. Dummy Arrays ..
323  double precision,dimension(z1-z0+1),intent(in) :: d_iso
324  double precision,dimension(r0:r1-1,z0-1:z1-1),intent(in) :: efield
325  double precision,dimension(r1-r0,z1-z0+4),intent(in) :: data_array
326  double precision,dimension(r1-r0,z1-z0+1),intent(out) :: mat_z
327 ! ..
328 ! .. Local Variables ..
329  integer :: iz,length
330 ! ..
331 ! .. Local Arrays ..
332  double precision,dimension(r1-r0,z1-z0+1) :: aux,psi_p,pij
333  double precision,dimension(r1-r0,z1-z0+3) :: mat_diff
334 !-----------------------------------------------
335  length = z1 - z0 + 1
336  mat_diff(:,1:length+2) = data_array(:,2:length+3) - &
337  data_array(:,1:length+2)
338 
339  where (abs(mat_diff(:,2:length+1)) > verysmall)
340  ! psi_p = (p - p_min) / (p_plus- p)
341  where (efield>zero)
342  ! p_min =data_array(ir-k1,iz-k2)
343  ! p =data_array(ir,iz)
344  ! p_plus =data_array(ir+k1,iz+k2)
345 
346  psi_p = mat_diff(:,1:length) / mat_diff(:,2:length+1)
347  elsewhere ! efield=<ZERO
348  ! p =data_array(ir,iz)
349  ! p_plus =data_array(ir+ k1,iz+ k2)
350  ! p_2plus =data_array(ir+2*k1,iz+2*k2)
351 
352  psi_p = mat_diff(:,3:length+2) / mat_diff(:,2:length+1)
353  end where
354  ! psi_p = psi_MN (psi_p)
355  pij = psi_p
356  where ((psi_p > 0.4) .and. (psi_p < 4))
357  pij = (one + half * psi_p) / three
358  elsewhere (psi_p >= 4)
359  pij = one
360  else where (psi_p <= zero)
361  pij = zero
362  end where
363  elsewhere ! abs(mat_diff(:,2:length+1)) < verysmall
364  pij = one
365  end where
366 
367  ! sigmasigma is (sigma_{i+1,j} - sigma_{i,j})*/
368  ! sigmasigma = p - p_plus
369  aux = pij * mat_diff(:,2:length+1)
370  ! aux = p + pij * sigmasigma
371  where (efield>=0)
372  aux = data_array(:,2:length+1) + aux
373  elsewhere
374  aux = data_array(:,3:length+2) - aux
375  end where
376 
377  ! aux = efield * (p + pij * sigmasigma)
378  aux = aux * efield
379  do iz = 1,length
380  ! mat_z = d_iso * sigmasigma
381  mat_z(:,iz) = d_iso(iz) * mat_diff(:,iz+1)
382  end do
383  ! mat_z = efield * (p + pij * sigmasigma) + d_iso * sigmasigma
384  mat_z = aux - mat_z
385 
386  end subroutine mat_cdr_f_ad_z
387 
388  subroutine make_vecs_er (er_min,er,r0,r1,dr)
389 ! ..
390  implicit none
391 !
392 ! .. Parameters ..
393  double precision,parameter :: half = 0.5d0,one = 1.0d0
394 ! ..
395 ! .. Dummy Variables ..
396  integer,intent(in) :: r0,r1
397  double precision,intent(in) :: dr
398 ! ..
399 ! .. Dummy Arrays ..
400  double precision,dimension(r0:r1-1),intent(out) :: er_min
401  double precision,dimension(r0+1:r1),intent(out) :: er
402 ! ..
403 ! .. Local Variables ..
404  integer :: ir
405  double precision :: r_inv
406 ! ..
407 ! .............................................................
408  do ir = r0,r1-1
409  r_inv = one / ((ir+ half) * dr)
410  er_min(ir) = ir * r_inv
411  er(ir+1) = (ir+1) * r_inv
412  end do
413 
414  end subroutine make_vecs_er
Generated on Fri Jul 4 2014 15:40:19 for Streamer fluid modeling - An overview of ARCoS by   doxygen 1.8.2