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