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 validate_mod
15  USE stagger_mod
16  implicit none
17  private
18  type, public, extends(mhd_t):: solver_t
19  contains
20  procedure:: init
21  procedure:: update
22  procedure, nopass:: cast2solver
23  procedure:: p2u
24  procedure:: u2p
25  procedure:: e2e_th
26  procedure:: e_th2e
27  !procedure:: e2s
28  !procedure:: s2e
29  !procedure:: e2ss
30  !procedure:: ss2e
31  procedure:: log_density => void
32  procedure:: log_pressure => void
33  procedure:: gas_pressure_ngz
34  !procedure:: cell_gas_pressure
35  procedure:: gas_temperature
36  procedure:: gas_temperature_ngz
37  !procedure:: cell_gas_temperature
38  procedure:: velocity_magnitude => void
39  procedure:: magnetic_field_magnitude => void
40  procedure:: grav_potential => void
41  procedure:: apply_heating
42  procedure:: gas_velocity_vector
43  procedure:: gas_velocity_scalar
44  procedure:: compression_magnitude
45  procedure:: vorticity_magnitude
46  end type
47  type(solver_t), public:: solver
48 CONTAINS
49 
50 !=======================================================================
51 !> Organize calls to the extras and the hydro solver
52 !=======================================================================
53 SUBROUTINE init (self)
54  class(solver_t) :: self
55  !.....................................................................
56  call self%mhd_t%init
57  call self%extras_t%init
58  call validate%init
59 END SUBROUTINE init
60 
61 !=======================================================================
62 !> Organize calls to the extras and the hydro solver
63 !=======================================================================
64 SUBROUTINE update (self)
65  class(solver_t) :: self
66  associate(d=>self%mem(:,:,:,self%idx%d,self%it,1))
67  !.....................................................................
68  call self%extras_t%pre_update
69  call validate%check (self, d, 'before update')
70  call self%mhd_t%update
71  call validate%check (self, d, ' after update')
72  call self%extras_t%post_update
73  end associate
74 END SUBROUTINE update
75 
76 !===============================================================================
77 !> Cast a generic task_t to patch_t
78 !===============================================================================
79 FUNCTION cast2solver (task) RESULT(solver)
80  class(task_t), pointer:: task
81  class(solver_t), pointer:: solver
82  !.............................................................................
83  select type (task)
84  class is (solver_t)
85  solver => task
86  class default
87  nullify(solver)
88  call io%abort ('patch_t%cast: failed to cast a task to patch_t')
89  end select
90 END FUNCTION cast2solver
91 
92 !=======================================================================
93 !> Get velocities from momenta
94 !=======================================================================
95 SUBROUTINE p2u(self, U, it)
96  class(solver_t) :: self
97  real, dimension(:,:,:,:), pointer:: u
98  integer:: it
99  !.....................................................................
100  real, dimension(:,:,:,:), pointer:: dd
101  integer:: i
102  !---------------------------------------------------------------------
103  associate(d => self%mem(:,:,:,self%idx%d, it,1), &
104  p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
105  do i=1,3
106  u(:,:,:,i) = p(:,:,:,i)/d
107  end do
108  end associate
109 END SUBROUTINE p2u
110 
111 !=======================================================================
112 !> Put velocities to momenta
113 !=======================================================================
114 SUBROUTINE u2p(self, U, it)
115  class(solver_t) :: self
116  real, dimension(:,:,:,:), pointer:: u
117  integer:: it
118  !.....................................................................
119  real, dimension(:,:,:,:), pointer:: dd
120  integer:: i
121  !---------------------------------------------------------------------
122  associate(d => self%mem(:,:,:,self%idx%d ,it,1), &
123  p => self%mem(:,:,:,self%idx%px:self%idx%pz,it,1))
124  do i=1,3
125  p(:,:,:,i) = u(:,:,:,i)*d
126  end do
127  end associate
128 END SUBROUTINE u2p
129 
130 !=======================================================================
131 !> Get thermal energy per unit mass from entropy
132 !=======================================================================
133 SUBROUTINE e2e_th(self, E_th, it)
134  class(solver_t) :: self
135  integer:: it
136  !.....................................................................
137  real, dimension(:,:,:), pointer:: e_th
138  real:: g1
139  !---------------------------------------------------------------------
140  associate(d => self%mem(:,:,:,self%idx%d,it,1), &
141  e => self%mem(:,:,:,self%idx%e,it,1))
142  e_th = e/d
143  end associate
144 END SUBROUTINE e2e_th
145 
146 !=======================================================================
147 !> Put thermal energy to entropy per unit volume
148 !=======================================================================
149 SUBROUTINE e_th2e(self, E_th, it)
150  class(solver_t) :: self
151  integer:: it
152  !.....................................................................
153  real, dimension(:,:,:), pointer:: e_th
154  real:: g1
155  !---------------------------------------------------------------------
156  associate(d => self%mem(:,:,:,self%idx%d,it,1), &
157  e => self%mem(:,:,:,self%idx%e,it,1))
158  e = d*e_th
159  end associate
160 END SUBROUTINE e_th2e
161 
162 !===============================================================================
163 FUNCTION up(f,i) RESULT (g)
164  real, dimension(:,:,:), intent(in):: f
165  real, dimension(size(f,1),size(f,2),size(f,3)):: g
166  integer:: i
167  !.............................................................................
168  g = 0.5*(cshift(f,1,i)+f)
169 END FUNCTION
170 
171 !===============================================================================
172 FUNCTION gas_pressure_ngz (self) RESULT (pg)
173  class(solver_t):: self
174  real, dimension(self%n(1),self%n(2),self%n(3)):: pg
175  integer :: l(3), u(3)
176  !-----------------------------------------------------------------------------
177  l = self%li; u = self%ui
178  associate(d => self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%d ,self%it,1), &
179  px => self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%px,self%it,1), &
180  py => self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%py,self%it,1), &
181  pz => self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%pz,self%it,1), &
182  bx => self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%bx,self%it,1), &
183  by => self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%by,self%it,1), &
184  bz => self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%bz,self%it,1), &
185  e => self%mem(l(1):u(1),l(2):u(2),l(3):u(3),self%idx%e ,self%it,1))
186  if (self%gamma==1.0) then
187  pg = d(l(1):u(1),l(2):u(2),l(3):u(3))
188  else
189  pg = (self%gamma-1.0) &
190  *(e(l(1):u(1),l(2):u(2),l(3):u(3)) -0.5*( &
191  (px(l(1):u(1),l(2):u(2),l(3):u(3))**2 + &
192  py(l(1):u(1),l(2):u(2),l(3):u(3))**2 + &
193  pz(l(1):u(1),l(2):u(2),l(3):u(3))**2) &
194  /d(l(1):u(1),l(2):u(2),l(3):u(3)) + &
195  up(bx(l(1):u(1),l(2):u(2),l(3):u(3)),1)**2 + &
196  up(by(l(1):u(1),l(2):u(2),l(3):u(3)),2)**2 + &
197  up(bz(l(1):u(1),l(2):u(2),l(3):u(3)),3)**2))
198  end if
199  end associate
200 END FUNCTION gas_pressure_ngz
201 
202 !===============================================================================
203 !> temperature from ideal gas law, taking into account that the temperature is
204 !> normalised by k_B and mu (so that specific gas constant = 1)
205 !===============================================================================
206 FUNCTION gas_temperature (self, lnd, ss) RESULT (tmp)
207  class(solver_t):: self
208  real, dimension(:,:,:), pointer:: lnd, ss
209  real, dimension(self%gn(1),self%gn(2),self%gn(3)):: tmp
210  !-----------------------------------------------------------------------------
211  associate(d => self%mem(:,:,:, 1,self%it,1))
212  tmp = self%gas_pressure()/d
213  if (io%verbose > 2) then
214  if (any(tmp<0)) then
215  error stop "solver_t%gas_temperature: T<0"
216  end if
217  end if
218  end associate
219 END FUNCTION gas_temperature
220 
221 !===============================================================================
222 FUNCTION gas_temperature_ngz (self) RESULT (tmp)
223  class(solver_t):: self
224  real, dimension(self%n(1),self%n(2),self%n(3)):: tmp
225  integer :: l(3), u(3)
226  !---------------------------------------------------------------------
227  l = self%li; u = self%ui
228  associate(d => self%mem(l(1):u(1),l(2):u(2),l(3):u(3), 1,self%it,1))
229  tmp = self%gas_pressure_ngz()/d
230 
231  if (io%verbose > 2) then
232  if (any(tmp<0)) then
233  error stop "solver_t%gas_temperature: T<0"
234  end if
235  end if
236  end associate
237 END FUNCTION gas_temperature_ngz
238 
239 !===============================================================================
240 SUBROUTINE log_density (self, v)
241  class(solver_t):: self
242  real, dimension(:,:,:), pointer:: v
243  !-----------------------------------------------------------------------------
244  v = log(self%mem(:,:,:,self%id,self%it,1))
245 END SUBROUTINE
246 
247 !===============================================================================
248 SUBROUTINE log_pressure (self, lnd, ss, v)
249  class(solver_t):: self
250  real, dimension(:,:,:), pointer:: lnd, ss, v
251  !-----------------------------------------------------------------------------
252  v = log(self%gas_pressure())
253 END SUBROUTINE
254 
255 !===============================================================================
256 SUBROUTINE velocity_magnitude (self, v)
257  class(solver_t):: self
258  real, dimension(:,:,:), pointer:: v, d
259  real, dimension(:,:,:,:), pointer:: p, u
260  integer:: i
261  !-----------------------------------------------------------------------------
262  d => self%mem(:,:,:,self%idx%d,self%it,1)
263  p => self%mem(:,:,:,self%idx%px:self%idx%px,self%it,1)
264  v = 0.0
265  do i=1,3
266  v(:,:,:) = v(:,:,:) + (p(:,:,:,i)/d)**2
267  end do
268  v = sqrt(v)
269 END SUBROUTINE
270 
271 !===============================================================================
272 SUBROUTINE magnetic_field_magnitude (self, v)
273  class(solver_t):: self
274  real, dimension(:,:,:), pointer:: v
275  real, dimension(:,:,:,:), pointer:: b
276  integer:: i
277  !-----------------------------------------------------------------------------
278  associate(bx => self%mem(:,:,:,self%idx%bx,self%it,1), &
279  by => self%mem(:,:,:,self%idx%by,self%it,1), &
280  bz => self%mem(:,:,:,self%idx%bz,self%it,1))
281  v = up(bx,1)**2 + up(by,2)**2 + up(bz,3)**2
282  v = sqrt(v)
283  end associate
284 END SUBROUTINE
285 
286 !===============================================================================
287 SUBROUTINE void (self, v)
288  class(solver_t):: self
289  real, dimension(:,:,:), pointer:: v
290  !-----------------------------------------------------------------------------
291  v = 0.0
292 END SUBROUTINE
293 
294 !===============================================================================
295 SUBROUTINE apply_heating (self, q)
296  class(solver_t):: self
297  real, dimension(:,:,:):: q
298 END SUBROUTINE
299 
300 !===============================================================================
301 FUNCTION gas_velocity_vector (self) RESULT (v)
302  class(solver_t):: self
303  real(kind=KindScalarVar), dimension(:,:,:), pointer :: d
304  real(kind=KindScalarVar), dimension(:,:,:,:), pointer:: p
305  real(kind=KindScalarVar), dimension(self%gn(1),self%gn(2),self%gn(3),3):: v
306  integer:: i
307  !-----------------------------------------------------------------------------
308  d => self%mem(:,:,:,self%idx%d,self%it,1)
309  p => self%mem(:,:,:,self%idx%px:self%idx%pz,self%it,1)
310  do i=1,3
311  v(:,:,:,i) = p(:,:,:,i)/d
312  end do
313 END FUNCTION gas_velocity_vector
314 
315 !===============================================================================
316 FUNCTION gas_velocity_scalar (self, idir) RESULT (v)
317  class(solver_t):: self
318  real(kind=KindScalarVar), dimension(:,:,:), pointer :: d
319  real(kind=KindScalarVar), dimension(self%gn(1),self%gn(2),self%gn(3)):: v
320  integer:: idir
321  !-----------------------------------------------------------------------------
322  d => self%mem(:,:,:,self%idx%d,self%it,1)
323  select case (idir)
324  case (1)
325  v = self%mem(:,:,:,self%idx%px,self%it,1) / d
326  case (2)
327  v = self%mem(:,:,:,self%idx%py,self%it,1) / d
328  case (3)
329  v = self%mem(:,:,:,self%idx%pz,self%it,1) / d
330  case default
331  error stop "solver_mod::gas_velocity_scalar:: invalid value of idir"
332  end select
333 END FUNCTION gas_velocity_scalar
334 
335 !===============================================================================
336 !> Compute the compression
337 !===============================================================================
338 SUBROUTINE compression_magnitude (self, w)
339  class(solver_t):: self
340  real(kind=KindScalarVar), dimension(:,:,:):: w
341  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: vx, vy, vz
342  integer:: i
343  integer, save:: itimer=0
344  !-----------------------------------------------------------------------------
345  call trace%begin ('solver_t%compression_magnitude', itimer=itimer)
346  allocate (vx(self%gn(1),self%gn(2),self%gn(3)))
347  allocate (vy(self%gn(1),self%gn(2),self%gn(3)))
348  allocate (vz(self%gn(1),self%gn(2),self%gn(3)))
349  vx = self%mem(:,:,:,self%idx%px,self%it,1)/self%mem(:,:,:,self%idx%d,self%it,1)
350  vy = self%mem(:,:,:,self%idx%py,self%it,1)/self%mem(:,:,:,self%idx%d,self%it,1)
351  vz = self%mem(:,:,:,self%idx%pz,self%it,1)/self%mem(:,:,:,self%idx%d,self%it,1)
352  w = max(- stagger%ddx(self%ds(1),vx) &
353  - stagger%ddy(self%ds(2),vy) &
354  - stagger%ddz(self%ds(3),vz), 0.0)
355  deallocate (vx, vy, vz)
356  call trace%end (itimer)
357 END SUBROUTINE compression_magnitude
358 
359 !===============================================================================
360 !> Compute the vorticity
361 !===============================================================================
362 SUBROUTINE vorticity_magnitude (self, w)
363  class(solver_t):: self
364  real(kind=KindScalarVar), dimension(:,:,:):: w
365  real(kind=KindScalarVar), dimension(:,:,:), allocatable:: vx, vy, vz
366  integer:: i
367  integer, save:: itimer=0
368  !-----------------------------------------------------------------------------
369  call trace%begin ('solver_t%vorticity_magnitude', itimer=itimer)
370  allocate (vx(self%gn(1),self%gn(2),self%gn(3)))
371  allocate (vy(self%gn(1),self%gn(2),self%gn(3)))
372  allocate (vz(self%gn(1),self%gn(2),self%gn(3)))
373  vx = self%mem(:,:,:,self%idx%px,self%it,1)/self%mem(:,:,:,self%idx%d,self%it,1)
374  vy = self%mem(:,:,:,self%idx%py,self%it,1)/self%mem(:,:,:,self%idx%d,self%it,1)
375  vz = self%mem(:,:,:,self%idx%pz,self%it,1)/self%mem(:,:,:,self%idx%d,self%it,1)
376  w = sqrt((stagger%ddy(self%ds(2),vz)-stagger%ddz(self%ds(3),vy))**2 &
377  + (stagger%ddz(self%ds(3),vx)-stagger%ddx(self%ds(1),vz))**2 &
378  + (stagger%ddx(self%ds(1),vy)-stagger%ddy(self%ds(2),vx))**2)
379  deallocate (vx, vy, vz)
380  call trace%end (itimer)
381 END SUBROUTINE vorticity_magnitude
382 
383 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
6th order stagger operators, with self-test procedure
Definition: stagger_2nd.f90:4
Template module for tasks.
Definition: task_mod.f90:4