DISPATCH
solver_mod.f90
1 !===============================================================================
2 !> The solver_t data type extends the patch_t with a generic guard zone dnload
3 !> procedure, which covers the need for all simple, Cartesian mesh based
4 !> solvers, and eliminates the need for a dnload procedure in each one of them.
5 !> It also adds restart functionality when initializing patches.
6 !===============================================================================
7 MODULE solver_mod
8  USE io_mod
9  USE trace_mod
10  USE task_mod
11  USE patch_mod
12  USE kinds_mod
13  USE mhd_mod
14  USE vector_mod
15  USE vector_ops
16  USE validate_mod
17  implicit none
18  private
19  type, public, extends(mhd_t):: solver_t
20  contains
21  procedure:: init
22  procedure:: update
23  procedure, nopass:: cast2solver
24  procedure:: p2u
25  procedure:: u2p
26  procedure:: e2e_th
27  procedure:: e_th2e
28  procedure:: e2s
29  procedure:: s2e
30  procedure:: e2ss
31  procedure:: ss2e
32  procedure:: log_density => void
33  procedure:: log_pressure => void
34  !procedure:: gas_pressure_ngz
35  !procedure:: cell_gas_pressure
36  !procedure:: gas_pressure
37  !procedure:: gas_temperature
38  !procedure:: gas_temperature_ngz
39  !procedure:: cell_gas_temperature
40  procedure:: velocity_magnitude => void
41  procedure:: magnetic_field_magnitude => void
42  procedure:: grav_potential => void
43  procedure:: apply_heating
44  procedure:: compression_magnitude
45  procedure:: vorticity_magnitude
46  procedure:: gas_velocity_vector
47  procedure:: gas_velocity_scalar
48  end type
49  type(solver_t), public:: solver
50 CONTAINS
51 
52 !===============================================================================
53 !> Organize calls to the extras and the hydro solver.
54 !===============================================================================
55 SUBROUTINE init (self)
56  class(solver_t) :: self
57  !.............................................................................
58  call self%mhd_t%pre_init ! calls self%mhd_t%pre_init
59  call self%mhd_t%init ! calls self%gpatch_t%init
60  call validate%init
61 END SUBROUTINE init
62 
63 !===============================================================================
64 !> Cast a generic task_t to patch_t
65 !===============================================================================
66 FUNCTION cast2solver (task) RESULT(solver)
67  class(task_t), pointer:: task
68  class(solver_t), pointer:: solver
69  !.............................................................................
70  select type (task)
71  class is (solver_t)
72  solver => task
73  class default
74  nullify(solver)
75  call io%abort ('solver_t%cast2solver: failed to cast a task to solver_t')
76  end select
77 END FUNCTION cast2solver
78 
79 !===============================================================================
80 !> Organize calls to the extras and the hydro solver
81 !===============================================================================
82 SUBROUTINE update (self)
83  class(solver_t) :: self
84  associate(d=>self%mem(:,:,:,self%idx%d,self%it,1))
85  !.............................................................................
86  call trace%begin ('solver_t%update')
87  call self%mhd_t%pre_update
88  call validate%check (self, d, 'before update')
89  call self%mhd_t%update
90  call validate%check (self, d, ' after update')
91  call self%mhd_t%post_update
92  end associate
93  call trace%end()
94 END SUBROUTINE update
95 
96 !===============================================================================
97 !> Get velocities from momenta
98 !===============================================================================
99 SUBROUTINE p2u(self, U, it)
100  class(solver_t) :: self
101  real, dimension(:,:,:,:), pointer:: u
102  integer:: it
103  !.............................................................................
104  real, dimension(:,:,:,:), allocatable:: dd
105  !-----------------------------------------------------------------------------
106  associate(d => self%mem(:,:,:,self%idx%d, it,1), &
107  p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
108  allocate(dd(size(d,1),size(d,2),size(d,3),3))
109  if (self%kind(1:13) == 'stagger2e_pic') then
110  dd=up(log(d))
111  u = p/exp(dd)
112  else
113  dd=down(log(d))
114  u = p/exp(dd)
115  end if
116  deallocate(dd)
117  end associate
118 END SUBROUTINE p2u
119 
120 !=======================================================================
121 !> Put velocities to momenta
122 !=======================================================================
123 SUBROUTINE u2p(self, U, it)
124  class(solver_t) :: self
125  real, dimension(:,:,:,:), pointer:: u
126  integer:: it
127  !.............................................................................
128  real, dimension(:,:,:,:), allocatable:: dd
129  !-----------------------------------------------------------------------------
130  associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
131  p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
132  allocate(dd(size(d,1),size(d,2),size(d,3),3))
133  if (self%kind(1:13) == 'stagger2e_pic') then
134  dd=up(log(d))
135  p = u*exp(dd)
136  else
137  dd=down(log(d))
138  p = u*exp(dd)
139  end if
140  deallocate(dd)
141  end associate
142 END SUBROUTINE u2p
143 
144 !===============================================================================
145 !> Get thermal energy per unit mass from entropy
146 !===============================================================================
147 SUBROUTINE e2e_th(self, E_th, it)
148  class(solver_t) :: self
149  integer:: it
150  !.............................................................................
151  real, dimension(:,:,:), pointer:: e_th
152  real:: g1
153  !-----------------------------------------------------------------------------
154  associate(d => self%mem(:,:,:,self%idx%d,it,1), &
155  s => self%mem(:,:,:,self%idx%s,it,1))
156  g1 = self%gamma-1.0
157  e_th = d**g1*exp(s*g1/d)/g1
158  end associate
159 END SUBROUTINE e2e_th
160 
161 !===============================================================================
162 !> Put thermal energy to entropy per unit volume
163 !===============================================================================
164 SUBROUTINE e_th2e(self, E_th, it)
165  class(solver_t) :: self
166  integer:: it
167  !.............................................................................
168  real, dimension(:,:,:), pointer:: e_th
169  real:: g1
170  !-----------------------------------------------------------------------------
171  associate(d => self%mem(:,:,:,self%idx%d,it,1), &
172  s => self%mem(:,:,:,self%idx%s,it,1))
173  g1 = self%gamma-1.0
174  s = d*log(e_th*g1/d**g1)/g1
175  end associate
176 END SUBROUTINE e_th2e
177 
178 !===============================================================================
179 !> Get entropy per unit volume from energy per unit volume
180 !===============================================================================
181 SUBROUTINE e2s(self, s, it)
182  class(solver_t) :: self
183  integer:: it
184  !.............................................................................
185  real, dimension(:,:,:):: s
186  real:: g1
187  !-----------------------------------------------------------------------------
188  associate(d => self%mem(:,:,:,self%idx%d,it,1), &
189  e => self%mem(:,:,:,self%idx%s,it,1))
190  g1 = self%gamma-1.0
191  s = d*log(e*g1/d**self%gamma)/g1
192  end associate
193 END SUBROUTINE e2s
194 
195 !===============================================================================
196 !> Get entropy per unit mass from energy per unit volume
197 !===============================================================================
198 SUBROUTINE e2ss(self, ss, it)
199  class(solver_t) :: self
200  integer:: it
201  !.............................................................................
202  real, dimension(:,:,:):: ss
203  real:: g1
204  !-----------------------------------------------------------------------------
205  associate(d => self%mem(:,:,:,self%idx%d,it,1), &
206  e => self%mem(:,:,:,self%idx%s,it,1))
207  g1 = self%gamma-1.0
208  ss = log(e*g1/d**self%gamma)/g1
209  end associate
210 END SUBROUTINE e2ss
211 
212 !===============================================================================
213 !> Get thermal energy per unit volume from entropy per unit volume
214 !===============================================================================
215 SUBROUTINE s2e(self, s, it)
216  class(solver_t) :: self
217  integer:: it
218  !.............................................................................
219  real, dimension(:,:,:):: s
220  real:: g1
221  !-----------------------------------------------------------------------------
222  associate(d => self%mem(:,:,:,self%idx%d,it,1), &
223  e => self%mem(:,:,:,self%idx%s,it,1))
224  g1 = self%gamma-1.0
225  e = d**self%gamma*exp(s/d*g1)/g1
226  end associate
227 END SUBROUTINE s2e
228 
229 !===============================================================================
230 !> Get thermal energy per unit volume from entropy per unit mass
231 !===============================================================================
232 SUBROUTINE ss2e(self, ss, it)
233  class(solver_t) :: self
234  integer:: it
235  !.............................................................................
236  real, dimension(:,:,:):: ss
237  real:: g1
238  !-----------------------------------------------------------------------------
239  associate(d => self%mem(:,:,:,self%idx%d,it,1), &
240  e => self%mem(:,:,:,self%idx%s,it,1))
241  g1 = self%gamma-1.0
242  e = d**self%gamma*exp(ss*g1)/g1
243  end associate
244 END SUBROUTINE ss2e
245 
246 !===============================================================================
247 !> temperature from ideal gas law, taking into account that the temperature is
248 !> normalised by k_B and mu (so that specific gas constant = 1)
249 !===============================================================================
250 !FUNCTION gas_pressure (self, lnd, ss) RESULT (pg)
251 ! class(solver_t):: self
252 ! real, dimension(:,:,:), pointer:: lnd, ss, d, s
253 ! optional:: lnd, ss
254 ! real, dimension(self%gn(1),self%gn(2),self%gn(3)):: pg
255 ! !-----------------------------------------------------------------------------
256 ! if (present(lnd)) then
257 ! pg = exp(lnd*self%gamma)*exp(ss*(self%gamma-1d0))
258 ! else
259 ! d => self%mem(:,:,:,self%idx%d,self%it,1)
260 ! s => self%mem(:,:,:,self%idx%s,self%it,1)
261 ! pg = d**self%gamma*exp(s/d*(self%gamma-1d0))
262 ! end if
263 !END FUNCTION gas_pressure
264 
265 !===============================================================================
266 !> temperature from ideal gas law, taking into account that the temperature is
267 !> normalised by k_B and mu (so that specific gas constant = 1)
268 !===============================================================================
269 !FUNCTION gas_temperature (self, lnd, ss) RESULT (tmp)
270 ! class(solver_t):: self
271 ! real, dimension(:,:,:), pointer:: lnd, ss
272 ! real, dimension(self%gn(1),self%gn(2),self%gn(3)):: tmp
273 ! !-----------------------------------------------------------------------------
274 ! associate (d => self%mem(:,:,:, 1,self%it,1))
275 ! tmp = self%gas_pressure()/d
276 ! if (io%verbose > 2) then
277 ! if (any(tmp<0)) then
278 ! call io%abort ("solver_t%gas_temperature: T<0")
279 ! end if
280 ! end if
281 ! end associate
282 !END FUNCTION gas_temperature
283 
284 !===============================================================================
285 SUBROUTINE log_density (self, v)
286  class(solver_t):: self
287  real, dimension(:,:,:), pointer:: v
288  !-----------------------------------------------------------------------------
289  v = log(self%mem(:,:,:,self%id,self%it,1))
290 END SUBROUTINE
291 
292 !===============================================================================
293 SUBROUTINE log_pressure (self, lnd, ss, v)
294  class(solver_t):: self
295  real, dimension(:,:,:), pointer:: lnd, ss, v
296  !-----------------------------------------------------------------------------
297  v = log(self%gas_pressure())
298 END SUBROUTINE
299 
300 !===============================================================================
301 SUBROUTINE velocity_magnitude (self, v)
302  class(solver_t):: self
303  real, dimension(:,:,:), pointer:: v
304  real, dimension(:,:,:,:), pointer:: p
305  !-----------------------------------------------------------------------------
306  p => self%mem(:,:,:,self%idx%px:self%idx%px,self%it,1)
307  if (self%kind(1:13) == 'stagger2e_pic') then
308  v = norm(p/exp(up(log(self%mem(:,:,:,self%id,self%it,1)))))
309  else
310  v = norm(p/exp(down(log(self%mem(:,:,:,self%id,self%it,1)))))
311  end if
312 END SUBROUTINE
313 
314 !===============================================================================
315 SUBROUTINE magnetic_field_magnitude (self, v)
316  class(solver_t):: self
317  real, dimension(:,:,:), pointer:: v
318  real, dimension(:,:,:,:), pointer:: b
319  !-----------------------------------------------------------------------------
320  b => self%mem(:,:,:,self%idx%bx:self%idx%bx,self%it,1)
321  v = norm(b)
322 END SUBROUTINE
323 
324 !===============================================================================
325 SUBROUTINE void (self, v)
326  class(solver_t):: self
327  real, dimension(:,:,:), pointer:: v
328  !-----------------------------------------------------------------------------
329  v = 0.0
330 END SUBROUTINE
331 
332 !===============================================================================
333 SUBROUTINE apply_heating (self, q)
334  class(solver_t):: self
335  real, dimension(:,:,:):: q
336 END SUBROUTINE
337 
338 !===============================================================================
339 FUNCTION gas_velocity_vector (self) RESULT (v)
340  USE scalar_mod
341  class(solver_t):: self
342  real(kind=KindScalarVar), dimension(:,:,:), pointer :: d
343  real(kind=KindScalarVar), dimension(:,:,:), allocatable :: lnd
344  real(kind=KindScalarVar), dimension(:,:,:,:), pointer:: p
345  real(kind=KindScalarVar), dimension(:,:,:,:), allocatable:: ld, dd
346  real(kind=KindScalarVar), dimension(self%gn(1),self%gn(2),self%gn(3),3):: v
347  !-----------------------------------------------------------------------------
348  d => self%mem(:,:,:,self%idx%d,self%it,1)
349  p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%it,1)
350  call allocate_vectors_a (self%gn, ld, dd)
351  call allocate_scalars_a (self%gn, lnd)
352  if (self%kind(1:13) == 'stagger2e_pic') then
353  lnd = log(d)
354  ld = up(lnd)
355  dd = exp(ld)
356  v = p/dd
357  else
358  lnd = log(d)
359  ld = down(lnd)
360  dd = exp(ld)
361  v = p/dd
362  end if
363  call deallocate_vectors_a (ld, dd)
364  call deallocate_scalars_a (lnd)
365  nullify (d,p)
366 END FUNCTION gas_velocity_vector
367 
368 !===============================================================================
369 FUNCTION gas_velocity_scalar (self, idir) RESULT (v)
370  USE scalar_mod
371  class(solver_t):: self
372  integer:: idir
373  real(kind=KindScalarVar), dimension(:,:,:), pointer:: d
374  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: lnd, ld, dd
375  real(kind=KindScalarVar), dimension(self%gn(1),self%gn(2),self%gn(3)):: v
376  !-----------------------------------------------------------------------------
377  d => self%mem(:,:,:,self%idx%d,self%it,1)
378  call allocate_scalars_a (self%gn, lnd, ld, dd)
379  lnd = log(d)
380  if (self%kind(1:13) == 'stagger2e_pic') then
381  select case (idir)
382  case (1)
383  ld = xup(lnd)
384  dd = exp(ld)
385  v = self%mem(:,:,:,self%idx%px,self%it,1) / dd
386  case (2)
387  ld = yup(lnd)
388  dd = exp(ld)
389  v = self%mem(:,:,:,self%idx%py,self%it,1) / dd
390  case (3)
391  ld = zup(lnd)
392  dd = exp(ld)
393  v = self%mem(:,:,:,self%idx%pz,self%it,1) / dd
394  case default
395  call io%abort ("solver_mod::gas_velocity_scalar:: invalid value of idir")
396  end select
397  else
398  select case (idir)
399  case (1)
400  ld = xdn(lnd)
401  dd = exp(ld)
402  v = self%mem(:,:,:,self%idx%px,self%it,1) / dd
403  case (2)
404  ld = ydn(lnd)
405  dd = exp(ld)
406  v = self%mem(:,:,:,self%idx%py,self%it,1) / dd
407  case (3)
408  ld = zdn(lnd)
409  dd = exp(ld)
410  v = self%mem(:,:,:,self%idx%pz,self%it,1) / dd
411  case default
412  call io%abort ("solver_mod::gas_velocity_scalar:: invalid value of idir")
413  end select
414  end if
415  call deallocate_scalars_a (lnd, ld, dd)
416  nullify (d)
417 END FUNCTION gas_velocity_scalar
418 
419 !===============================================================================
420 !> Compute the compression, centered in cells, given that momenta are dn-staggered
421 !> in stagger2, except for stagger2e_pic, where it up-staggered
422 !===============================================================================
423 SUBROUTINE compression_magnitude (self, w)
424  class(solver_t):: self
425  real(kind=KindScalarVar), dimension(:,:,:):: w
426  real(kind=KindScalarVar), dimension(:,:,:,:), pointer:: d
427  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: vx, vy, vz
428  integer:: i
429  integer, save:: itimer=0
430  !-----------------------------------------------------------------------------
431  call trace%begin ('solver_t%compression_magnitude', itimer=itimer)
432  allocate (vx(self%gn(1),self%gn(2),self%gn(3)))
433  allocate (vy(self%gn(1),self%gn(2),self%gn(3)))
434  allocate (vz(self%gn(1),self%gn(2),self%gn(3)))
435  allocate ( d(self%gn(1),self%gn(2),self%gn(3),3))
436  if (self%kind(1:13) == 'stagger2e_pic') then
437  d = up(log(self%mem(:,:,:,self%idx%d,self%it,1)))
438  else
439  d = down(log(self%mem(:,:,:,self%idx%d,self%it,1)))
440  end if
441  vx = self%mem(:,:,:,self%idx%px,self%it,1)/exp(d(:,:,:,1))
442  vy = self%mem(:,:,:,self%idx%py,self%it,1)/exp(d(:,:,:,2))
443  vz = self%mem(:,:,:,self%idx%pz,self%it,1)/exp(d(:,:,:,3))
444  if (self%kind(1:13) == 'stagger2e_pic') then
445  w = max(- ddxdn(self%ds,vx) &
446  - ddydn(self%ds,vy) &
447  - ddzdn(self%ds,vz), 0.0)
448  else
449  w = max(- ddxup(self%ds,vx) &
450  - ddyup(self%ds,vy) &
451  - ddzup(self%ds,vz), 0.0)
452  end if
453  deallocate (vx, vy, vz, d)
454  call trace%end (itimer)
455 END SUBROUTINE compression_magnitude
456 
457 !===============================================================================
458 !> Compute the vorticity
459 !===============================================================================
460 SUBROUTINE vorticity_magnitude (self, w)
461  class(solver_t):: self
462  real(kind=KindScalarVar), dimension(:,:,:):: w
463  real(kind=KindScalarVar), dimension(:,:,:,:), pointer:: d
464  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: vx, vy, vz
465  integer:: i
466  integer, save:: itimer=0
467  !-----------------------------------------------------------------------------
468  call trace%begin ('solver_t%vorticity_magnitude', itimer=itimer)
469  allocate (vx(self%gn(1),self%gn(2),self%gn(3)))
470  allocate (vy(self%gn(1),self%gn(2),self%gn(3)))
471  allocate (vz(self%gn(1),self%gn(2),self%gn(3)))
472  allocate ( d(self%gn(1),self%gn(2),self%gn(3),3))
473  if (self%kind(1:13) == 'stagger2e_pic') then
474  d = up(log(self%mem(:,:,:,self%idx%d,self%it,1)))
475  else
476  d = down(log(self%mem(:,:,:,self%idx%d,self%it,1)))
477  end if
478  vx = self%mem(:,:,:,self%idx%px,self%it,1)/exp(d(:,:,:,1))
479  vy = self%mem(:,:,:,self%idx%py,self%it,1)/exp(d(:,:,:,2))
480  vz = self%mem(:,:,:,self%idx%pz,self%it,1)/exp(d(:,:,:,3))
481  if (self%kind(1:13) == 'stagger2e_pic') then
482  w = sqrt(ydn(zdn(ddyup(self%ds,vz)-ddzup(self%ds,vy)))**2 &
483  + zdn(xdn(ddzup(self%ds,vx)-ddxup(self%ds,vz)))**2 &
484  + xdn(ydn(ddxup(self%ds,vy)-ddyup(self%ds,vx)))**2)
485  else
486  w = sqrt(yup(zup(ddydn(self%ds,vz)-ddzdn(self%ds,vy)))**2 &
487  + zup(xup(ddzdn(self%ds,vx)-ddxdn(self%ds,vz)))**2 &
488  + xup(yup(ddxdn(self%ds,vy)-ddydn(self%ds,vx)))**2)
489  end if
490  deallocate (vx, vy, vz)
491  call trace%end (itimer)
492 END SUBROUTINE vorticity_magnitude
493 
494 END MODULE solver_mod
Generic validation module. The general idea is to be able to compare two runs at critical points in t...
Template module for patches, which adds pointers to memory and mesh, and number of dimensions and var...
Definition: patch_mod.f90:6
RAMSES Godunov solvers, use of guard zones; specifically in HLLD.
Definition: mhd_mod.f90:16
Definition: io_mod.f90:4
This module contains all experiment specific information necessary to solve the heat diffusion proble...
Definition: solver_mod.f90:5
Template module for tasks.
Definition: task_mod.f90:4