MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_amr_solution_node.t
Go to the documentation of this file.
2  use mod_comm_lib, only: mpistop
3 
4  implicit none
5  private
6 
7  public :: getnode, putnode
8  public :: alloc_node, alloc_state
9  public :: dealloc_node
10 
11 contains
12 
13 
14  !> Get first available igrid on processor ipe
15  integer function getnode(ipe)
16  use mod_forest, only: igrid_inuse
18 
19  integer, intent(in) :: ipe
20  integer :: igrid, igrid_available
21 
22  igrid_available=0
23 
24  do igrid=1,max_blocks
25  if (igrid_inuse(igrid,ipe)) cycle
26 
27  igrid_available=igrid
28  exit
29  end do
30 
31  if (igrid_available == 0) then
32  getnode = -1
33  print *, "Current maximum number of grid blocks:", max_blocks
34  call mpistop("Insufficient grid blocks; increase max_blocks in meshlist")
35  else
36  getnode=igrid_available
37  igrid_inuse(igrid,ipe)=.true.
38  end if
39 
40  if (ipe==mype) then
41  ! initialize nodal block
42  node(1:nodehi,getnode) = 0
43  rnode(1:rnodehi,getnode) = zero
44  end if
45 
46  end function getnode
47 
48  ! put igrid on processor ipe to be not in use
49  subroutine putnode(igrid,ipe)
50  use mod_forest
51  implicit none
52 
53  integer, intent(in) :: igrid, ipe
54 
55  igrid_inuse(igrid,ipe)=.false.
56 
57  end subroutine putnode
58 
59  !> allocate arrays on igrid node
60  subroutine alloc_node(igrid)
61  use mod_forest
63  use mod_geometry
66  use mod_b0, only: set_b0_grid
67 
68  integer, intent(in) :: igrid
69 
70  integer :: level, ig^d, ign^d, ixcog^l, ix, i^d
71  integer :: imin, imax, index, igco^d, ixshift, offset, ifirst
72  integer :: icase, ixgext^l
73  double precision :: dx^d, summeddx, sizeuniformpart^d
74  double precision :: xext(ixglo^d-1:ixghi^d+1,1:ndim)
75  double precision :: delx_ext(ixglo1-1:ixghi1+1)
76  double precision :: exp_factor_ext(ixglo1-1:ixghi1+1),del_exp_factor_ext(ixglo1-1:ixghi1+1),exp_factor_primitive_ext(ixglo1-1:ixghi1+1)
77  double precision :: xc(ixglo1:ixghi1),delxc(ixglo1:ixghi1)
78  double precision :: exp_factor_coarse(ixglo1:ixghi1),del_exp_factor_coarse(ixglo1:ixghi1),exp_factor_primitive_coarse(ixglo1:ixghi1)
79 
80  ixcogmin^d=1;
81  ixcogmax^d=(ixghi^d-2*nghostcells)/2+2*nghostcells;
82 
83  icase=mod(nghostcells,2)
84  if(stagger_grid) icase=1
85  select case(icase)
86  case(0)
87  ixgext^l=ixg^ll;
88  case(1)
89  ! for ghost cell related prolongations, we need
90  ! an extra layer with known volumes and dx-intervals
91  ! in case the number of ghost cells is odd
92  ixgext^l=ixg^ll^ladd1;
93  case default
94  call mpistop("no such case")
95  end select
96 
97  ! set level information
98  level=igrid_to_node(igrid,mype)%node%level
99 
100  if(.not. allocated(ps(igrid)%w)) then
101 
102  ! allocate arrays for solution and space
103  call alloc_state(igrid, ps(igrid), ixg^ll, ixgext^l, .true.)
104  ! allocate arrays for one level coarser solution
105  call alloc_state_coarse(igrid, psc(igrid), ixcog^l, ixcog^l)
106  if(.not.convert) then
107  ! allocate arrays for temp solution 1
108  call alloc_state(igrid, ps1(igrid), ixg^ll, ixgext^l, .false.)
109 
110  ! allocate temporary solution space
111  select case (t_integrator)
113  call alloc_state(igrid, ps2(igrid), ixg^ll, ixgext^l, .false.)
114  case(rk3_bt,rk4,ssprk5,imex_cb3a)
115  call alloc_state(igrid, ps2(igrid), ixg^ll, ixgext^l, .false.)
116  call alloc_state(igrid, ps3(igrid), ixg^ll, ixgext^l, .false.)
117  case(imex_ars3,imex_232)
118  call alloc_state(igrid, ps2(igrid), ixg^ll, ixgext^l, .false.)
119  call alloc_state(igrid, ps3(igrid), ixg^ll, ixgext^l, .false.)
120  call alloc_state(igrid, ps4(igrid), ixg^ll, ixgext^l, .false.)
121  end select
122  end if
123 
124  end if
125 
126  ! avoid dividing by zero rho in skipped corner ghostcells when phys_req_diagonal=F
127  ps(igrid)%w(:^d&,1)=1.d0
128  ps(igrid)%level=level
129  psc(igrid)%level=level-1
130  ! avoid dividing by zero rho in skipped corner ghostcells when phys_req_diagonal=F
131  psc(igrid)%w(:^d&,1)=1.d0
132  if(phys_trac) ps(igrid)%special_values=0.d0
133  if(.not.convert) then
134  ps1(igrid)%level=level
135  select case (t_integrator)
137  ps2(igrid)%level=level
138  case(rk3_bt,rk4,ssprk5,imex_cb3a)
139  ps2(igrid)%level=level
140  ps3(igrid)%level=level
141  case(imex_ars3,imex_232)
142  ps2(igrid)%level=level
143  ps3(igrid)%level=level
144  ps4(igrid)%level=level
145  end select
146  end if
147 
148  ! block pointer to current block
149  block=>ps(igrid)
150 
151  ig^d=igrid_to_node(igrid,mype)%node%ig^d;
152 
153  node(plevel_,igrid)=level
154  ^d&node(pig^d_,igrid)=ig^d\
155 
156  ! set dx information
157  ^d&rnode(rpdx^d_,igrid)=dx(^d,level)\
158  dxlevel(:)=dx(:,level)
159 
160  ! uniform cartesian case as well as all unstretched coordinates
161  ! determine the minimal and maximal corners
162  ^d&rnode(rpxmin^d_,igrid)=xprobmin^d+dble(ig^d-1)*dg^d(level)\
163  ^d&rnode(rpxmax^d_,igrid)=xprobmin^d+dble(ig^d)*dg^d(level)\
164  {if(rnode(rpxmax^d_,igrid)>xprobmax^d) rnode(rpxmax^d_,igrid)=xprobmax^d\}
165 
166  ^d&dx^d=rnode(rpdx^d_,igrid)\
167  {do ix=ixglo^d,ixmhi^d-nghostcells
168  ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmin^d_,igrid)+(dble(ix-nghostcells)-half)*dx^d
169  end do\}
170  ! update overlap cells of neighboring blocks in the same way to get the same values
171  {do ix=ixmhi^d-nghostcells+1,ixghi^d
172  ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmax^d_,igrid)+(dble(ix-ixmhi^d)-half)*dx^d
173  end do\}
174 
175  ^d&dx^d=2.0d0*rnode(rpdx^d_,igrid)\
176  {do ix=ixcogmin^d,ixcogmax^d
177  psc(igrid)%x(ix^d%ixCoG^s,^d)=rnode(rpxmin^d_,igrid)+(dble(ix-nghostcells)-half)*dx^d
178  end do\}
179 
180  ^d&ps(igrid)%dx(ixgext^s,^d)=rnode(rpdx^d_,igrid);
181  ^d&psc(igrid)%dx(ixcog^s,^d)=2.0d0*rnode(rpdx^d_,igrid);
182  ^d&dx^d=rnode(rpdx^d_,igrid)\
183  {do ix=ixgextmin^d,ixgextmax^d
184  xext(ix^d%ixGext^s,^d)=rnode(rpxmin^d_,igrid)+(dble(ix-nghostcells)-half)*dx^d
185  end do\}
186 
187  if(any(stretched_dim)) then
188  {if(stretch_type(^d) == stretch_uni)then
189  imin=(ig^d-1)*block_nx^d
190  imax=ig^d*block_nx^d
191  rnode(rpxmin^d_,igrid)=xprobmin^d+dxfirst_1mq(level,^d) &
192  *(1.0d0-qstretch(level,^d)**imin)
193  rnode(rpxmax^d_,igrid)=xprobmin^d+dxfirst_1mq(level,^d) &
194  *(1.0d0-qstretch(level,^d)**imax)
195  ! fix possible out of bound due to precision
196  if(rnode(rpxmax^d_,igrid)>xprobmax^d) rnode(rpxmax^d_,igrid)=xprobmax^d
197  ixshift=(ig^d-1)*block_nx^d-nghostcells
198  do ix=ixgextmin^d,ixgextmax^d
199  index=ixshift+ix
200  ps(igrid)%dx(ix^d%ixGext^s,^d)=dxfirst(level,^d)*qstretch(level,^d)**(index-1)
201  enddo
202  igco^d=(ig^d-1)/2
203  ixshift=igco^d*block_nx^d+(1-modulo(ig^d,2))*block_nx^d/2-nghostcells
204  do ix=ixcogmin^d,ixcogmax^d
205  index=ixshift+ix
206  psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxfirst(level-1,^d)*qstretch(level-1,^d)**(index-1)
207  psc(igrid)%x(ix^d%ixCoG^s,^d)=xprobmin^d+dxfirst_1mq(level-1,^d)&
208  *(1.0d0-qstretch(level-1,^d)**(index-1))&
209  + 0.5d0*dxfirst(level-1,^d)*qstretch(level-1,^d)**(index-1)
210  end do
211  ! now that dx and grid boundaries are known: fill cell centers
212  ifirst=nghostcells+1
213  ! first fill the mesh
214  summeddx=0.0d0
215  do ix=ixmlo^d,ixmhi^d
216  ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmin^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
217  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
218  enddo
219  ! then ghost cells to left
220  summeddx=0.0d0
221  do ix=nghostcells,1,-1
222  ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmin^d_,igrid)-summeddx-0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
223  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
224  enddo
225  ! then ghost cells to right
226  summeddx=0.0d0
227  do ix=ixghi^d-nghostcells+1,ixghi^d
228  ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmax^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
229  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
230  enddo
231  select case(icase)
232  case(0)
233  ! if even number of ghost cells: xext is just copy of local x
234  xext(ixgext^s,^d)=ps(igrid)%x(ixgext^s,^d)
235  case(1)
236  ! if uneven number of ghost cells: extra layer left/right
237  summeddx=0.0d0
238  do ix=ixmlo^d,ixmhi^d
239  xext(ix^d%ixGext^s,^d)=rnode(rpxmin^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
240  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
241  enddo
242  ! then ghost cells to left
243  summeddx=0.0d0
244  do ix=nghostcells,ixgextmin^d,-1
245  xext(ix^d%ixGext^s,^d)=rnode(rpxmin^d_,igrid)-summeddx-0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
246  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
247  enddo
248  ! then ghost cells to right
249  summeddx=0.0d0
250  do ix=ixghi^d-nghostcells+1,ixgextmax^d
251  xext(ix^d%ixGext^s,^d)=rnode(rpxmax^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
252  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
253  enddo
254  case default
255  call mpistop("no such case")
256  end select
257  endif\}
258  {if(stretch_type(^d) == stretch_symm)then
259  ! here we distinguish three kinds of grid blocks
260  ! depending on their ig-index, set per level
261  ! the first n_stretchedblocks/2 will stretch to the left
262  ! the middle ntotal-n_stretchedblocks will be uniform
263  ! the last n_stretchedblocks/2 will stretch to the right
264  if(ig^d<=nstretchedblocks(level,^d)/2)then
265  ! stretch to the left
266  offset=block_nx^d*nstretchedblocks(level,^d)/2
267  imin=(ig^d-1)*block_nx^d
268  imax=ig^d*block_nx^d
269  rnode(rpxmin^d_,igrid)=xprobmin^d+xstretch^d-dxfirst_1mq(level,^d) &
270  *(1.0d0-qstretch(level,^d)**(offset-imin))
271  rnode(rpxmax^d_,igrid)=xprobmin^d+xstretch^d-dxfirst_1mq(level,^d) &
272  *(1.0d0-qstretch(level,^d)**(offset-imax))
273  ! fix possible out of bound due to precision
274  if(rnode(rpxmin^d_,igrid)<xprobmin^d) rnode(rpxmin^d_,igrid)=xprobmin^d
275  ixshift=(ig^d-1)*block_nx^d-nghostcells
276  do ix=ixgextmin^d,ixgextmax^d
277  index=ixshift+ix
278  ps(igrid)%dx(ix^d%ixGext^s,^d)=dxfirst(level,^d)*qstretch(level,^d)**(offset-index)
279  enddo
280  ixshift=(nstretchedblocks(level,^d)/2-ig^d)*(block_nx^d/2)+block_nx^d/2+nghostcells
281  do ix=ixcogmin^d,ixcogmax^d
282  index=ixshift-ix
283  psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxfirst(level-1,^d)*qstretch(level-1,^d)**index
284  enddo
285  ! last block: to modify ghost cells!!!
286  if(ig^d==nstretchedblocks(level,^d)/2)then
287  if(ng^d(level)==nstretchedblocks(level,^d))then
288  ! if middle blocks do not exist then use symmetry
289  do ix=ixghi^d-nghostcells+1,ixgextmax^d
290  ps(igrid)%dx(ix^d%ixGext^s,^d)= &
291  ps(igrid)%dx(2*(ixghi^d-nghostcells)+1-ix^d%ixGext^s,^d)
292  enddo
293  do ix=ixcogmax^d-nghostcells+1,ixcogmax^d
294  psc(igrid)%dx(ix^d%ixCoG^s,^d)= &
295  psc(igrid)%dx(2*(ixcogmax^d-nghostcells)+1-ix^d%ixCoG^s,^d)
296  enddo
297  else
298  ! if middle blocks exist then use same as middle blocks:
299  do ix=ixghi^d-nghostcells+1,ixgextmax^d
300  ps(igrid)%dx(ix^d%ixGext^s,^d)=dxmid(level,^d)
301  enddo
302  do ix=ixcogmax^d-nghostcells+1,ixcogmax^d
303  psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxmid(level-1,^d)
304  enddo
305  endif
306  endif
307  ! first block: make ghost cells symmetric (to allow periodicity)
308  if(ig^d==1)then
309  do ix=ixgextmin^d,nghostcells
310  ps(igrid)%dx(ix^d%ixGext^s,^d)=ps(igrid)%dx(2*nghostcells+1-ix^d%ixGext^s,^d)
311  enddo
312  do ix=1,nghostcells
313  psc(igrid)%dx(ix^d%ixCoG^s,^d)=psc(igrid)%dx(2*nghostcells+1-ix^d%ixCoG^s,^d)
314  enddo
315  endif
316  else
317  if(ig^d<=ng^d(level)-nstretchedblocks(level,^d)/2) then
318  ! keep uniform
319  ps(igrid)%dx(ixgext^s,^d)=dxmid(level,^d)
320  psc(igrid)%dx(ixcog^s,^d)=dxmid(level-1,^d)
321  rnode(rpxmin^d_,igrid)=xprobmin^d+xstretch^d+(ig^d-nstretchedblocks(level,^d)/2-1)*block_nx^d*dxmid(level,^d)
322  rnode(rpxmax^d_,igrid)=xprobmin^d+xstretch^d+(ig^d-nstretchedblocks(level,^d)/2) *block_nx^d*dxmid(level,^d)
323  ! first and last block: to modify the ghost cells!!!
324  if(ig^d==nstretchedblocks(level,^d)/2+1)then
325  do ix=ixgextmin^d,nghostcells
326  ps(igrid)%dx(ix^d%ixGext^s,^d)=dxfirst(level,^d)*qstretch(level,^d)**(nghostcells-ix)
327  enddo
328  do ix=1,nghostcells
329  psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxfirst(level-1,^d)*qstretch(level-1,^d)**(nghostcells-ix)
330  enddo
331  endif
332  if(ig^d==ng^d(level)-nstretchedblocks(level,^d))then
333  do ix=ixghi^d-nghostcells+1,ixgextmax^d
334  ps(igrid)%dx(ix^d%ixGext^s,^d)=dxfirst(level,^d)*qstretch(level,^d)**(ix-block_nx^d-nghostcells-1)
335  enddo
336  do ix=ixcogmax^d-nghostcells+1,ixcogmax^d
337  psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxfirst(level-1,^d)*qstretch(level-1,^d)**(ix-ixcogmax^d+nghostcells-1)
338  enddo
339  endif
340  else
341  ! stretch to the right
342  offset=block_nx^d*(ng^d(level)-nstretchedblocks(level,^d)/2)
343  sizeuniformpart^d=dxmid(1,^d)*(domain_nx^d-nstretchedblocks_baselevel(^d)*block_nx^d)
344  imin=(ig^d-1)*block_nx^d-offset
345  imax=ig^d*block_nx^d-offset
346  rnode(rpxmin^d_,igrid)=xprobmin^d+xstretch^d+sizeuniformpart^d+dxfirst_1mq(level,^d) &
347  *(1.0d0-qstretch(level,^d)**imin)
348  rnode(rpxmax^d_,igrid)=xprobmin^d+xstretch^d+sizeuniformpart^d+dxfirst_1mq(level,^d) &
349  *(1.0d0-qstretch(level,^d)**imax)
350  ! fix possible out of bound due to precision
351  if(rnode(rpxmax^d_,igrid)>xprobmax^d) rnode(rpxmax^d_,igrid)=xprobmax^d
352  ixshift=(ig^d-1)*block_nx^d-nghostcells-offset
353  do ix=ixgextmin^d,ixgextmax^d
354  index=ixshift+ix
355  ps(igrid)%dx(ix^d%ixGext^s,^d)=dxfirst(level,^d)*qstretch(level,^d)**(index-1)
356  enddo
357  ixshift=(ig^d+nstretchedblocks(level,^d)/2-ng^d(level)-1)*(block_nx^d/2)-nghostcells
358  do ix=ixcogmin^d,ixcogmax^d
359  index=ixshift+ix
360  psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxfirst(level-1,^d)*qstretch(level-1,^d)**(index-1)
361  enddo
362  ! first block: modify ghost cells!!!
363  if(ig^d==ng^d(level)-nstretchedblocks(level,^d)+1)then
364  if(ng^d(level)==nstretchedblocks(level,^d))then
365  ! if middle blocks do not exist then use symmetry
366  do ix=ixgextmin^d,nghostcells
367  ps(igrid)%dx(ix^d%ixGext^s,^d)=ps(igrid)%dx(2*nghostcells+1-ix^d%ixGext^s,^d)
368  enddo
369  do ix=1,nghostcells
370  psc(igrid)%dx(ix^d%ixCoG^s,^d)=psc(igrid)%dx(2*nghostcells+1-ix^d%ixCoG^s,^d)
371  enddo
372  else
373  ! if middle blocks exist then use same as middle blocks:
374  do ix=ixgextmin^d,nghostcells
375  ps(igrid)%dx(ix^d%ixGext^s,^d)=dxmid(level,^d)
376  enddo
377  do ix=1,nghostcells
378  psc(igrid)%dx(ix^d%ixCoG^s,^d)=dxmid(level-1,^d)
379  enddo
380  endif
381  endif
382  ! last block: make ghost cells symmetric (to allow periodicity)
383  if(ig^d==ng^d(level))then
384  do ix=ixghi^d-nghostcells+1,ixgextmax^d
385  ps(igrid)%dx(ix^d%ixGext^s,^d)=ps(igrid)%dx(2*(ixghi^d-nghostcells)+1-ix^d%ixGext^s,^d)
386  enddo
387  do ix=ixcogmax^d-nghostcells+1,ixcogmax^d
388  psc(igrid)%dx(ix^d%ixCoG^s,^d)=psc(igrid)%dx(2*(ixcogmax^d-nghostcells)+1-ix^d%ixCoG^s,^d)
389  enddo
390  endif
391  endif
392  endif
393  ! now that dx and grid boundaries are known: fill cell centers
394  ifirst=nghostcells+1
395  ! first fill the mesh
396  summeddx=0.0d0
397  do ix=ixmlo^d,ixmhi^d
398  ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmin^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
399  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
400  enddo
401  ! then ghost cells to left
402  summeddx=0.0d0
403  do ix=nghostcells,1,-1
404  ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmin^d_,igrid)-summeddx-0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
405  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
406  enddo
407  ! then ghost cells to right
408  summeddx=0.0d0
409  do ix=ixghi^d-nghostcells+1,ixghi^d
410  ps(igrid)%x(ix^d%ixG^t,^d)=rnode(rpxmax^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixG^t,^d)
411  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
412  enddo
413  ! and next for the coarse representation
414  ! first fill the mesh
415  summeddx=0.0d0
416  do ix=nghostcells+1,ixcogmax^d-nghostcells
417  psc(igrid)%x(ix^d%ixCoG^s,^d)=rnode(rpxmin^d_,igrid)+summeddx+0.5d0*psc(igrid)%dx(ix^d%ixCoG^s,^d)
418  summeddx=summeddx+psc(igrid)%dx(ix^d%ifirst,^d)
419  enddo
420  ! then ghost cells to left
421  summeddx=0.0d0
422  do ix=nghostcells,1,-1
423  psc(igrid)%x(ix^d%ixCoG^s,^d)=rnode(rpxmin^d_,igrid)-summeddx-0.5d0*psc(igrid)%dx(ix^d%ixCoG^s,^d)
424  summeddx=summeddx+psc(igrid)%dx(ix^d%ifirst,^d)
425  enddo
426  ! then ghost cells to right
427  summeddx=0.0d0
428  do ix=ixcogmax^d-nghostcells+1,ixcogmax^d
429  psc(igrid)%x(ix^d%ixCoG^s,^d)=rnode(rpxmax^d_,igrid)+summeddx+0.5d0*psc(igrid)%dx(ix^d%ixCoG^s,^d)
430  summeddx=summeddx+psc(igrid)%dx(ix^d%ifirst,^d)
431  enddo
432  select case(icase)
433  case(0)
434  ! if even number of ghost cells: xext is just copy of local x
435  xext(ixgext^s,^d)=ps(igrid)%x(ixgext^s,^d)
436  case(1)
437  ! if uneven number of ghost cells: extra layer left/right
438  summeddx=0.0d0
439  do ix=ixmlo^d,ixmhi^d
440  xext(ix^d%ixGext^s,^d)=rnode(rpxmin^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
441  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
442  enddo
443  ! then ghost cells to left
444  summeddx=0.0d0
445  do ix=nghostcells,ixgextmin^d,-1
446  xext(ix^d%ixGext^s,^d)=rnode(rpxmin^d_,igrid)-summeddx-0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
447  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
448  enddo
449  ! then ghost cells to right
450  summeddx=0.0d0
451  do ix=ixghi^d-nghostcells+1,ixgextmax^d
452  xext(ix^d%ixGext^s,^d)=rnode(rpxmax^d_,igrid)+summeddx+0.5d0*ps(igrid)%dx(ix^d%ixGext^s,^d)
453  summeddx=summeddx+ps(igrid)%dx(ix^d%ifirst,^d)
454  enddo
455  case default
456  call mpistop("no such case")
457  end select
458  endif\}
459  endif
460 
461  ! calculate area of cell surfaces for standard block
462  call get_surface_area(ps(igrid),ixg^ll)
463  ! calculate area of cell surfaces for coarser representative block
464  call get_surface_area(psc(igrid),ixcog^l)
465  ! calculate volume and distance of cells
466  ps(igrid)%dsC=1.d0
467  select case (coordinate)
468  case (cartesian)
469  ps(igrid)%dvolume(ixgext^s)= {^d&rnode(rpdx^d_,igrid)|*}
470  ps(igrid)%ds(ixgext^s,1:ndim)=ps(igrid)%dx(ixgext^s,1:ndim)
471  ps(igrid)%dsC(ixgext^s,1:ndim)=ps(igrid)%dx(ixgext^s,1:ndim)
472  psc(igrid)%dvolume(ixcog^s)= {^d&2.d0*rnode(rpdx^d_,igrid)|*}
473  psc(igrid)%ds(ixcog^s,1:ndim)=psc(igrid)%dx(ixcog^s,1:ndim)
474  case (cartesian_stretched)
475  ps(igrid)%dvolume(ixgext^s)= {^d&ps(igrid)%dx(ixgext^s,^d)|*}
476  ps(igrid)%ds(ixgext^s,1:ndim)=ps(igrid)%dx(ixgext^s,1:ndim)
477  ps(igrid)%dsC(ixgext^s,1:ndim)=ps(igrid)%dx(ixgext^s,1:ndim)
478  psc(igrid)%dvolume(ixcog^s)= {^d&psc(igrid)%dx(ixcog^s,^d)|*}
479  psc(igrid)%ds(ixcog^s,1:ndim)=psc(igrid)%dx(ixcog^s,1:ndim)
480  case (cartesian_expansion)
481  {^ifoned
482  delx_ext(ixgext^s) = ps(igrid)%dx(ixgext^s,1)
483  if(associated(usr_set_surface)) call usr_set_surface(ixgext^l,xext(ixgext^s,1),delx_ext(ixgext^s),exp_factor_ext(ixgext^s),del_exp_factor_ext(ixgext^s),exp_factor_primitive_ext(ixgext^s))
484  ps(igrid)%dvolume(ixgext^s)= exp_factor_primitive_ext(ixgext^s)
485  ps(igrid)%ds(ixgext^s,1)=ps(igrid)%dx(ixgext^s,1)
486  ps(igrid)%dsC(ixgext^s,1)=ps(igrid)%dx(ixgext^s,1)
487  xc(ixcog^s) = psc(igrid)%x(ixcog^s,1)
488  delxc(ixcog^s) = psc(igrid)%dx(ixcog^s,1)
489  if(associated(usr_set_surface)) call usr_set_surface(ixcog^l,xc(ixcog^s),delxc(ixcog^s),exp_factor_coarse(ixcog^s),del_exp_factor_coarse(ixcog^s),exp_factor_primitive_coarse(ixcog^s))
490  psc(igrid)%dvolume(ixcog^s)= exp_factor_primitive_coarse(ixcog^s)
491  psc(igrid)%ds(ixcog^s,1)=psc(igrid)%dx(ixcog^s,1)
492  }
493  case (spherical)
494  ps(igrid)%dvolume(ixgext^s)=(xext(ixgext^s,1)**2 &
495  +ps(igrid)%dx(ixgext^s,1)**2/12.0d0)*&
496  ps(igrid)%dx(ixgext^s,1){^nooned &
497  *two*dabs(dsin(xext(ixgext^s,2))) &
498  *dsin(half*ps(igrid)%dx(ixgext^s,2))}{^ifthreed*ps(igrid)%dx(ixgext^s,3)}
499  psc(igrid)%dvolume(ixcog^s)=(psc(igrid)%x(ixcog^s,1)**2 &
500  +psc(igrid)%dx(ixcog^s,1)**2/12.0d0)*&
501  psc(igrid)%dx(ixcog^s,1){^nooned &
502  *two*dabs(dsin(psc(igrid)%x(ixcog^s,2))) &
503  *dsin(half*psc(igrid)%dx(ixcog^s,2))}{^ifthreed*psc(igrid)%dx(ixcog^s,3)}
504  ps(igrid)%ds(ixgext^s,1)=ps(igrid)%dx(ixgext^s,1)
505  {^nooned ps(igrid)%ds(ixgext^s,2)=xext(ixgext^s,1)*ps(igrid)%dx(ixgext^s,2)}
506  {^ifthreed ps(igrid)%ds(ixgext^s,3)=xext(ixgext^s,1)*dsin(xext(ixgext^s,2))*&
507  ps(igrid)%dx(ixgext^s,3)}
508  ps(igrid)%dsC(ixgext^s,1)=ps(igrid)%dx(ixgext^s,1)
509  {^nooned ps(igrid)%dsC(ixgext^s,2)=(xext(ixgext^s,1)+half*ps(igrid)%dx(ixgext^s,1))*&
510  ps(igrid)%dx(ixgext^s,2)
511  if(ndir>ndim) then
512  ps(igrid)%dsC(ixgext^s,3)=(xext(ixgext^s,1)+half*ps(igrid)%dx(ixgext^s,1))*&
513  dsin(xext(ixgext^s,2)+half*ps(igrid)%dx(ixgext^s,2))
514  end if
515  }
516  {^ifthreed ps(igrid)%dsC(ixgext^s,3)=(xext(ixgext^s,1)+half*ps(igrid)%dx(ixgext^s,1))*&
517  dsin(xext(ixgext^s,2)+half*ps(igrid)%dx(ixgext^s,2))*&
518  ps(igrid)%dx(ixgext^s,3)}
519  case (cylindrical)
520  ps(igrid)%dvolume(ixgext^s)=dabs(xext(ixgext^s,1)) &
521  *ps(igrid)%dx(ixgext^s,1){^de&*ps(igrid)%dx(ixgext^s,^de) }
522  psc(igrid)%dvolume(ixcog^s)=dabs(psc(igrid)%x(ixcog^s,1)) &
523  *psc(igrid)%dx(ixcog^s,1){^de&*psc(igrid)%dx(ixcog^s,^de) }
524  ps(igrid)%ds(ixgext^s,r_)=ps(igrid)%dx(ixgext^s,r_)
525  ps(igrid)%dsC(ixgext^s,r_)=ps(igrid)%dx(ixgext^s,r_)
526  if(z_>0.and.z_<=ndim) then
527  ps(igrid)%ds(ixgext^s,z_)=ps(igrid)%dx(ixgext^s,z_)
528  ps(igrid)%dsC(ixgext^s,z_)=ps(igrid)%dx(ixgext^s,z_)
529  if(phi_>z_.and.ndir>ndim) then
530  ps(igrid)%dsC(ixgext^s,phi_)=xext(ixgext^s,1)+half*ps(igrid)%dx(ixgext^s,1)
531  end if
532  end if
533  if(phi_>0.and.phi_<=ndim) then
534  ps(igrid)%ds(ixgext^s,phi_)=xext(ixgext^s,1)*ps(igrid)%dx(ixgext^s,phi_)
535  ps(igrid)%dsC(ixgext^s,phi_)=(xext(ixgext^s,1)+&
536  half*ps(igrid)%dx(ixgext^s,1))*ps(igrid)%dx(ixgext^s,phi_)
537  if(z_>phi_.and.ndir>ndim) ps(igrid)%dsC(ixgext^s,z_)=1.d0
538  end if
539  case default
540  call mpistop("Sorry, coordinate unknown")
541  end select
542 
543  ! initialize background non-evolving solution
544  if (b0field) call set_b0_grid(igrid)
545  if (number_equi_vars>0) call phys_set_equi_vars(igrid)
546 
547  ! find the blocks on the boundaries
548  ps(igrid)%is_physical_boundary=.false.
549  {
550  do i^d=-1,1
551  if(i^d==0) cycle
552  ign^d=ig^d+i^d
553  ! blocks at periodic boundary have neighbors in the physical domain
554  ! thus threated at internal blocks with no physical boundary
555  if (periodb(^d)) ign^d=1+modulo(ign^d-1,ng^d(level))
556  if (ign^d > ng^d(level)) then
557  if(phi_ > 0 .and. poleb(2,^d)) then
558  ! if at a pole, the boundary is not physical boundary
559  ps(igrid)%is_physical_boundary(2*^d)=.false.
560  else
561  ps(igrid)%is_physical_boundary(2*^d)=.true.
562  end if
563  else if (ign^d < 1) then
564  if(phi_ > 0 .and. poleb(1,^d)) then
565  ! if at a pole, the boundary is not physical boundary
566  ps(igrid)%is_physical_boundary(2*^d-1)=.false.
567  else
568  ps(igrid)%is_physical_boundary(2*^d-1)=.true.
569  end if
570  end if
571  end do
572  \}
573  if(any(ps(igrid)%is_physical_boundary)) then
574  phyboundblock(igrid)=.true.
575  else
576  phyboundblock(igrid)=.false.
577  end if
578 
579  end subroutine alloc_node
580 
581  !> allocate memory to physical state of igrid node
582  subroutine alloc_state(igrid, s, ixG^L, ixGext^L, alloc_once_for_ps)
584  type(state) :: s
585  integer, intent(in) :: igrid, ixg^l, ixgext^l
586  logical, intent(in) :: alloc_once_for_ps
587  integer :: ixgs^l
588 
589  allocate(s%w(ixg^s,1:nw))
590  s%igrid=igrid
591  s%w=0.d0
592  s%ixG^l=ixg^l;
593  {^d& ixgsmin^d = ixgmin^d-1; ixgsmax^d = ixgmax^d|;}
594  if(stagger_grid) then
595  allocate(s%ws(ixgs^s,1:nws))
596  s%ws=0.d0
597  if(record_electric_field) then
598  allocate(s%we(ixgs^s,1:nws))
599  s%we=0.d0
600  end if
601  s%ixGs^l=ixgs^l;
602  end if
603  if(alloc_once_for_ps) then
604  ! allocate extra variables for ps state
605  if(nw_extra>0) allocate(s%wextra(ixg^s,1:nw_extra))
606  ! allocate coordinates
607  allocate(s%x(ixg^s,1:ndim))
608  allocate(s%dx(ixgext^s,1:ndim), &
609  s%ds(ixgext^s,1:ndim),s%dsC(ixgext^s,1:3))
610  allocate(s%dvolume(ixgext^s))
611  allocate(s%surfaceC(ixgs^s,1:ndim), &
612  s%surface(ixg^s,1:ndim))
613  ! allocate physical boundary flag
614  allocate(s%is_physical_boundary(2*ndim))
615  if(local_timestep) then
616  allocate(s%dt(ixg^s))
617  endif
618 
619  if(b0field) then
620  allocate(s%B0(ixg^s,1:ndir,0:ndim))
621  allocate(s%J0(ixg^s,7-2*ndir:3))
622  end if
623  if(number_equi_vars > 0) then
624  allocate(s%equi_vars(ixg^s,1:number_equi_vars,0:ndim))
625  endif
626 
627  ! allocate space for special values for each block state
628  if(phys_trac) then
629  ! special_values(1) Tcoff local
630  ! special_values(2) Tmax local
631  ! special_values(3:2+ndim) Bdir local
632  allocate(s%special_values(ndim+2))
633  end if
634  else
635  ! share common info from ps states to save memory
636  if(nw_extra>0) s%wextra=>ps(igrid)%wextra
637  s%x=>ps(igrid)%x
638  s%dx=>ps(igrid)%dx
639  s%ds=>ps(igrid)%ds
640  s%dsC=>ps(igrid)%dsC
641  s%dvolume=>ps(igrid)%dvolume
642  s%surfaceC=>ps(igrid)%surfaceC
643  s%surface=>ps(igrid)%surface
644  s%is_physical_boundary=>ps(igrid)%is_physical_boundary
645  if(b0field) then
646  s%B0=>ps(igrid)%B0
647  s%J0=>ps(igrid)%J0
648  end if
649  if(number_equi_vars > 0) then
650  s%equi_vars=>ps(igrid)%equi_vars
651  endif
652  if(phys_trac) s%special_values=>ps(igrid)%special_values
653  end if
654  end subroutine alloc_state
655 
656  !> allocate memory to one-level coarser physical state of igrid node
657  subroutine alloc_state_coarse(igrid, s, ixG^L, ixGext^L)
659  type(state) :: s
660  integer, intent(in) :: igrid, ixG^L, ixGext^L
661  integer :: ixGs^L
662 
663  allocate(s%w(ixg^s,1:nw))
664  s%igrid=igrid
665  s%w=0.d0
666  s%ixG^l=ixg^l;
667  {^d& ixgsmin^d = ixgmin^d-1; ixgsmax^d = ixgmax^d|;}
668  if(stagger_grid) then
669  allocate(s%ws(ixgs^s,1:nws))
670  s%ws=0.d0
671  s%ixGs^l=ixgs^l;
672  end if
673  if(b0fieldalloccoarse) then
674  allocate(s%B0(ixg^s,1:ndir,0:ndim))
675  end if
676  ! allocate coordinates
677  allocate(s%x(ixg^s,1:ndim))
678  allocate(s%dx(ixgext^s,1:ndim), &
679  s%ds(ixgext^s,1:ndim),s%dsC(ixgext^s,1:3))
680  allocate(s%dvolume(ixgext^s))
681  allocate(s%surfaceC(ixgs^s,1:ndim), &
682  s%surface(ixg^s,1:ndim))
683  ! allocate physical boundary flag
684  allocate(s%is_physical_boundary(2*ndim))
685  end subroutine alloc_state_coarse
686 
687  subroutine dealloc_state(igrid, s,dealloc_x)
689  integer, intent(in) :: igrid
690  type(state) :: s
691  logical, intent(in) :: dealloc_x
692 
693  deallocate(s%w)
694  if(stagger_grid) then
695  deallocate(s%ws)
696  if(record_electric_field) deallocate(s%we)
697  end if
698  if(dealloc_x) then
699  if(nw_extra>0) deallocate(s%wextra)
700  ! deallocate coordinates
701  deallocate(s%x)
702  deallocate(s%dx,s%ds,s%dsC)
703  deallocate(s%dvolume)
704  deallocate(s%surfaceC,s%surface)
705  deallocate(s%is_physical_boundary)
706  if(b0field) then
707  deallocate(s%B0)
708  deallocate(s%J0)
709  end if
710  if(number_equi_vars > 0) then
711  deallocate(s%equi_vars)
712  end if
713  else
714  nullify(s%x,s%dx,s%ds,s%dsC,s%dvolume,s%surfaceC,s%surface)
715  nullify(s%is_physical_boundary)
716  if(b0field) nullify(s%B0,s%J0)
717  if(number_equi_vars > 0) then
718  nullify(s%equi_vars)
719  end if
720  if(nw_extra>0) nullify(s%wextra)
721  end if
722  end subroutine dealloc_state
723 
724  subroutine dealloc_state_coarse(igrid, s)
726  integer, intent(in) :: igrid
727  type(state) :: s
728 
729  deallocate(s%w)
730  if(stagger_grid) then
731  deallocate(s%ws)
732  end if
733  if(b0fieldalloccoarse) then
734  deallocate(s%B0)
735  end if
736  ! deallocate coordinates
737  deallocate(s%x)
738  deallocate(s%dx,s%ds,s%dsC)
739  deallocate(s%dvolume)
740  deallocate(s%surfaceC,s%surface)
741  deallocate(s%is_physical_boundary)
742  end subroutine dealloc_state_coarse
743 
744  subroutine dealloc_node(igrid)
746 
747  integer, intent(in) :: igrid
748 
749  if (igrid==0) then
750  call mpistop("trying to delete a non-existing grid in dealloc_node")
751  end if
752 
753  call dealloc_state(igrid, ps(igrid),.true.)
754  call dealloc_state_coarse(igrid, psc(igrid))
755  if(.not.convert) then
756  call dealloc_state(igrid, ps1(igrid),.false.)
757  ! deallocate temporary solution space
758  select case (t_integrator)
760  call dealloc_state(igrid, ps2(igrid),.false.)
761  case(rk3_bt,rk4,ssprk5,imex_cb3a)
762  call dealloc_state(igrid, ps2(igrid),.false.)
763  call dealloc_state(igrid, ps3(igrid),.false.)
764  case(imex_ars3,imex_232)
765  call dealloc_state(igrid, ps2(igrid),.false.)
766  call dealloc_state(igrid, ps3(igrid),.false.)
767  call dealloc_state(igrid, ps4(igrid),.false.)
768  end select
769  end if
770 
771  end subroutine dealloc_node
772 
773 end module mod_amr_solution_node
subroutine alloc_state_coarse(igrid, s, ixGL, ixGextL)
allocate memory to one-level coarser physical state of igrid node
subroutine, public putnode(igrid, ipe)
subroutine dealloc_state_coarse(igrid, s)
subroutine, public dealloc_node(igrid)
subroutine, public alloc_node(igrid)
allocate arrays on igrid node
subroutine, public alloc_state(igrid, s, ixGL, ixGextL, alloc_once_for_ps)
allocate memory to physical state of igrid node
integer function, public getnode(ipe)
Get first available igrid on processor ipe.
subroutine dealloc_state(igrid, s, dealloc_x)
Definition: mod_b0.t:1
subroutine, public set_b0_grid(igrid)
Definition: mod_b0.t:9
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
Module with basic grid data structures.
Definition: mod_forest.t:2
logical, dimension(:,:), allocatable, save igrid_inuse
Definition: mod_forest.t:70
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
Definition: mod_forest.t:32
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
integer ixghi
Upper index of grid block arrays.
integer, parameter imex_232
double precision, dimension(:), allocatable dg
extent of grid blocks in domain per dimension, in array over levels
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical stagger_grid
True for using stagger grid.
integer, parameter imex_ars3
integer, parameter imex_trapezoidal
integer mype
The rank of the current MPI task.
integer, parameter plevel_
integer, dimension(:), allocatable, parameter d
logical local_timestep
each cell has its own timestep or not
integer, parameter nodehi
grid hierarchy info (level and grid indices)
integer ndir
Number of spatial dimensions (components) for vector variables.
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
logical phys_trac
Use TRAC (Johnston 2019 ApJL, 873, L22) for MHD or 1D HD.
logical convert
If true and restart_from_file is given, convert snapshots to other file formats.
integer, parameter rnodehi
grid location info (corner coordinates and grid spacing)
integer max_blocks
The maximum number of grid blocks in a processor.
integer, parameter imex_cb3a
logical record_electric_field
True for record electric field.
integer, parameter imex_222
integer t_integrator
time integrator method
integer, dimension(:,:), allocatable node
double precision, dimension(ndim) dxlevel
integer, parameter imex_midpoint
integer, parameter ixglo
Lower index of grid block arrays (always 1)
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
procedure(sub_set_equi_vars), pointer phys_set_equi_vars
Definition: mod_physics.t:90
Module with all the methods that users can customize in AMRVAC.
procedure(set_surface), pointer usr_set_surface