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