MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_amr_fct.t
Go to the documentation of this file.
1 module mod_amr_fct
2  implicit none
3  private
4 
5  type facealloc
6  double precision, dimension(:^D&), pointer :: face
7  end type facealloc
8 
9  type fake_neighbors
10  integer :: igrid
11  integer :: ipe
12  end type fake_neighbors
13 
14  type(facealloc), dimension(:,:,:), allocatable, public :: pface
15 
16  type(fake_neighbors), dimension(:^D&,:,:), allocatable, public :: fine_neighbors
17 
18  integer, dimension(:,:^D&,:), allocatable, public :: old_neighbor
19 
20  integer :: itag, isend, irecv
21  integer :: nrecv, nsend, ibuf_recv, ibuf_send, ibuf_send_next
22  integer, dimension(^ND) :: isize
23  integer, dimension(:), allocatable :: recvrequest, sendrequest
24  integer, dimension(:,:), allocatable :: recvstatus, sendstatus
25  double precision, allocatable :: recvbuffer(:), sendbuffer(:)
26 
27  public :: store_faces
28  public :: comm_faces
29  public :: end_comm_faces
30  public :: deallocatebfaces
31  public :: old_neighbors
32  public :: prolong_2nd_stg
33  public :: already_fine
34 
35 contains
36  !> This subroutine performs a 2nd order prolongation for a staggered field F,
37  !> preserving the divergence of the coarse cell.
38  !> This is useful for preserving DivF=0.
39  !> If DivF=f(x), a different algorithm must be used.
40  subroutine prolong_2nd_stg(sCo,sFi,ixCo^Lin,ixFi^Lin,dxCo^D,xComin^D,dxFi^D,xFimin^D,ghost,fine_^Lin)
42  use mod_physics
43 
44  logical, intent(in) :: ghost
45  integer, intent(in) :: ixco^lin, ixfi^lin
46  double precision, intent(in) :: dxco^d, xcomin^d, dxfi^d, xfimin^d
47  type(state), intent(in) :: sco
48  type(state), intent(inout) :: sfi
49 
50  logical, optional :: fine_^lin
51  logical :: fine_^l
52 
53  double precision :: eta^d, invdxco^d
54  integer :: ixco^l,ixfi^l
55  integer :: idim1,idim2,ix^de,idim3,ixfis^l,ixgs^l,ixcos^l,ixfisc^l
56  integer :: ixcosv^l(1:ndim),ixfisv^l(1:ndim)
57  integer :: hxcos^l,jxcos^l,ixcose^l,ixfise^l,hxfisc^l,jxfisc^l,ipxfisc^l,ixcosc^l,imxfisc^l,jpxfisc^l,jmxfisc^l,hpxfisc^l
58  integer :: hxfi^l,jxfi^l,hijxfi^l,hjixfi^l,hjjxfi^l
59  integer :: iihxfi^l,iijxfi^l,ijhxfi^l,ijjxfi^l,ihixfi^l,ijixfi^l,ihjxfi^l
60  integer :: jihxfi^l,jijxfi^l,jjhxfi^l,jjjxfi^l,jhixfi^l,jjixfi^l,jhjxfi^l
61  double precision :: bfluxco(sco%ixgs^s,nws),bfluxfi(sfi%ixgs^s,nws)
62  double precision :: slopes(sco%ixgs^s,ndim),b_energy_change(ixg^t)
63  {^ifthreed
64  ! Directional bias, see pdf
65  ! These arrays should have ranges ixO^S
66  ! Since they are still not known, we give ranges ixG^T,
67  ! even though it is wasteful
68  double precision :: sigmau(ixg^t),sigmad(ixg^t),sigma(ixg^t,1:ndim), alpha(ixg^t,1:ndim)
69  ! Auxiliary arrays for magnetic fluxes
70  double precision :: f1(ixg^t),f2(ixg^t),f3(ixg^t),f4(ixg^t)
71  }
72 
73  {^ifoned
74  call mpistop("CT prolongation not implemented in 1D. But CT is not needed.")
75  }
76 
77  {^nooned
78  ! Note on the indices:
79  ! ixCo Cells where
80  ! divergence-preserving prolongation will be applied.
81  ! ixFi Fine cells which correspond to that extent.
82  ! ixCoE For 'expanded', faces that need to be used to calculate
83  ! the slopes for the prolongation
84  ! ixCos
85  ! ixFis
86  ! ixCosV And
87  ! ixFisV For 'variable', since the ranges depend on the direction of the faces
88  ! ixCosC For 'copy',
89  ! ixFisC For 'copy', to fill the information from the fine grid needed in the internal faces prolongation step.
90  ! At the end, ixFisC is used to copy the assign the magnetic
91  ! fields components to the state structure.
92 
93  fine_^l=.false.;
94 
95  {if(present(fine_min^din)) fine_min^d=fine_min^din;}
96  {if(present(fine_max^din)) fine_max^d=fine_max^din;}
97 
98  ! When filling ghost cells, ixFi^L are given.
99  ! When refining the block, ixCo^L are given.
100  if (ghost) then
101  ixfi^l=ixfi^lin;
102  {ixcomin^d=int((ixfimin^d+nghostcells+1)/2);}
103  {ixcomax^d=int((ixfimax^d+nghostcells+1)/2);}
104  else
105  ixco^l=ixco^lin;
106  ixfi^l=ixm^ll;
107  end if
108 
109  ! Expanded range for staggered variables
110  ixcosmin^d=ixcomin^d-1;
111  ixcosmax^d=ixcomax^d;
112 
113  ixfismin^d=ixfimin^d-1;
114  ixfismax^d=ixfimax^d;
115 
116 
117  associate(wcos=>sco%ws, wfis=>sfi%ws,wco=>sco%w, wfi=>sfi%w)
118  ! Assemble general indices
119  ixgsmin^d=sfi%ixGsmin^d;
120  ixgsmax^d=sfi%ixGsmax^d;
121 
122  do idim1=1,ndim
123  ixcosvmin^d(idim1)=ixcomin^d-kr(^d,idim1);
124  ixcosvmax^d(idim1)=ixcomax^d;
125  ixfisvmin^d(idim1)=ixfimin^d-kr(^d,idim1);
126  ixfisvmax^d(idim1)=ixfimax^d;
127  end do
128 
129  ! Initialize auxiliary arrays at zero
130  bfluxco = zero
131  bfluxfi = zero
132  slopes = zero
133 
134 
135  invdxco^d=1.d0/dxco^d;
136 
137  ! Fill coarse magnetic flux array
138  do idim1=1,ndim
139  ! Fill information from parts already at the fine level
140  ! First set up indices
141  if (ghost) then
142  ixfisemin^d=max(1-kr(^d,idim1),ixfisvmin^d(idim1)-2*(1-kr(^d,idim1)));
143  ixfisemax^d=min(ixgsmax^d,ixfisvmax^d(idim1)+2*(1-kr(^d,idim1)));
144  ixcose^l=int((ixfise^l+nghostcells+1)/2);
145  else
146  ixcose^l=ixcosv^l(idim1)^ladd(1-kr(idim1,^d));
147  ixfise^l=ixfisv^l(idim1)^ladd2*(1-kr(idim1,^d));
148  end if
149  ! Convert fine fields to fluxes
150  bfluxfi(ixfise^s,idim1)=wfis(ixfise^s,idim1)*sfi%surfaceC(ixfise^s,idim1)
151  {^iftwod
152  idim2=1+mod(idim1,2)
153  }
154  {^ifthreed
155  idim2=1+mod(idim1,3)
156  idim3=1+mod(idim1+1,3)
157  }
158  bfluxco(ixcose^s,idim1) = zero
159  ! Add fine fluxes sharing the same fine face
160  {do ix^de=0,1\}
161  ixfisc^l=ixfise^l+ix2*kr(idim2,^d);
162  {^ifthreed
163  ixfisc^l=ixfisc^l+ix3*kr(idim3,^d);
164  }
165  bfluxco(ixcose^s,idim1)=bfluxco(ixcose^s,idim1)+bfluxfi(ixfiscmin^d:ixfiscmax^d:2,idim1)
166  {end do^DE&\}
167  end do
168 
169  ! Omit indices for already refined face, if any
170  {
171  if (fine_min^d) then
172  ixcosvmin^d(^d)=ixcosvmin^d(^d)+1
173  ixfisvmin^d(^d)=ixfisvmin^d(^d)+2
174  end if
175  if (fine_max^d) then
176  ixcosvmax^d(^d)=ixcosvmax^d(^d)-1
177  ixfisvmax^d(^d)=ixfisvmax^d(^d)-2
178  end if
179  \}
180 
181  do idim1=1,ndim
182  ixcose^l=ixcosv^l(idim1);
183  ! Omit part already refined
184  {
185  if (^d/=idim1) then
186  if ((.not.fine_min^d).or.(.not.ghost)) then
187  ixcosemin^d=ixcosvmin^d(idim1)-1
188  end if
189  if ((.not.fine_max^d).or.(.not.ghost)) then
190  ixcosemax^d=ixcosvmax^d(idim1)+1
191  end if
192  end if
193  \}
194  ! Fill coarse flux array from coarse field
195  bfluxco(ixcose^s,idim1)=wcos(ixcose^s,idim1)*sco%surfaceC(ixcose^s,idim1)
196  end do
197  ! Finished filling coarse flux array
198 
199  ! Distribute coarse fluxes among fine fluxes
200  ! There are too many loops here, perhaps optimize later
201  do idim1=1,ndim
202  ixcos^l=ixcosv^l(idim1);
203  do idim2=1,ndim
204  if(idim1==idim2) cycle
205  ! Calculate slope in direction idim2
206  ! Set up indices
207  jxcos^l=ixcos^l+kr(idim2,^d);
208  hxcos^l=ixcos^l-kr(idim2,^d);
209  slopes(ixcos^s,idim2)=0.125d0*(bfluxco(jxcos^s,idim1)-bfluxco(hxcos^s,idim1))
210  {^iftwod
211  do ix2=0,1
212  ixfiscmin^d=ixfisvmin^d(idim1)+ix2*kr(^d,idim2);
213  ixfiscmax^d=ixfisvmax^d(idim1)+ix2*kr(^d,idim2);
214  bfluxfi(ixfiscmin^d:ixfiscmax^d:2,idim1)=half*(bfluxco(ixcos^s,idim1)&
215  +(2*ix2-1)*slopes(ixcos^s,idim2))
216  end do
217  }
218  {^ifthreed
219  do idim3=1,ndim
220  ! Calculate slope in direction idim3
221  ! Set up indices
222  jxcos^l=ixcos^l+kr(idim3,^d);
223  hxcos^l=ixcos^l-kr(idim3,^d);
224  slopes(ixcos^s,idim3)=0.125d0*(bfluxco(jxcos^s,idim1)-bfluxco(hxcos^s,idim1))
225  if(lvc(idim1,idim2,idim3)<1) cycle
226  do ix2=0,1
227  do ix3=0,1
228  {ixfiscmin^d=ixfisvmin^d(idim1)+ix2*kr(^d,idim2)+ix3*kr(^d,idim3);}
229  {ixfiscmax^d=ixfisvmax^d(idim1)+ix2*kr(^d,idim2)+ix3*kr(^d,idim3);}
230  bfluxfi(ixfiscmin^d:ixfiscmax^d:2,idim1)=quarter*bfluxco(ixcos^s,idim1)&
231  +quarter*(2*ix2-1)*slopes(ixcos^s,idim2)&
232  +quarter*(2*ix3-1)*slopes(ixcos^s,idim3)
233  end do
234  end do
235  end do
236  }
237  end do
238  end do
239 
240  ! Calculate interior fine fluxes
241  {^ifthreed
242  ! Directional bias for nonlinear prolongation
243  do idim1=1,ndim
244  do idim2=1,ndim
245  do idim3=1,ndim
246  if (lvc(idim1,idim2,idim3)<1) cycle
247  ! Set up indices
248  hxfi^l=ixfi^l-kr(idim1,^d);
249  jxfi^l=ixfi^l+kr(idim1,^d);
250 
251  hijxfi^l=hxfi^l+kr(idim3,^d);
252  hjixfi^l=hxfi^l+kr(idim2,^d);
253  hjjxfi^l=hijxfi^l+kr(idim2,^d);
254 
255  iihxfi^l=ixfi^l-kr(idim3,^d);
256  iijxfi^l=ixfi^l+kr(idim3,^d);
257  ihixfi^l=ixfi^l-kr(idim2,^d);
258  ihjxfi^l=ihixfi^l+kr(idim3,^d);
259  ijixfi^l=ixfi^l+kr(idim2,^d);
260  ijhxfi^l=ijixfi^l-kr(idim3,^d);
261  ijjxfi^l=ijixfi^l+kr(idim3,^d);
262 
263  jihxfi^l=jxfi^l-kr(idim3,^d);
264  jijxfi^l=jxfi^l+kr(idim3,^d);
265  jhixfi^l=jxfi^l-kr(idim2,^d);
266  jhjxfi^l=jhixfi^l+kr(idim3,^d);
267  jjixfi^l=jxfi^l+kr(idim2,^d);
268  jjhxfi^l=jjixfi^l-kr(idim3,^d);
269  jjjxfi^l=jjixfi^l+kr(idim3,^d);
270 
271  sigmau(ixco^s)=&
272  abs(bfluxfi(jhixfimin^d:jhixfimax^d:2,idim2))&
273  + abs(bfluxfi(jhjxfimin^d:jhjxfimax^d:2,idim2))&
274  + abs(bfluxfi(jjixfimin^d:jjixfimax^d:2,idim2))&
275  + abs(bfluxfi(jjjxfimin^d:jjjxfimax^d:2,idim2))&
276  + abs(bfluxfi(jihxfimin^d:jihxfimax^d:2,idim3))&
277  + abs(bfluxfi(jijxfimin^d:jijxfimax^d:2,idim3))&
278  + abs(bfluxfi(jjhxfimin^d:jjhxfimax^d:2,idim3))&
279  + abs(bfluxfi(jjjxfimin^d:jjjxfimax^d:2,idim3))
280 
281  sigmad(ixco^s)=&
282  abs(bfluxfi(ihixfimin^d:ihixfimax^d:2,idim2))&
283  + abs(bfluxfi(ihjxfimin^d:ihjxfimax^d:2,idim2))&
284  + abs(bfluxfi(ijixfimin^d:ijixfimax^d:2,idim2))&
285  + abs(bfluxfi(ijjxfimin^d:ijjxfimax^d:2,idim2))&
286  + abs(bfluxfi(iihxfimin^d:iihxfimax^d:2,idim3))&
287  + abs(bfluxfi(iijxfimin^d:iijxfimax^d:2,idim3))&
288  + abs(bfluxfi(ijhxfimin^d:ijhxfimax^d:2,idim3))&
289  + abs(bfluxfi(ijjxfimin^d:ijjxfimax^d:2,idim3))
290 
291  sigma(ixco^s,idim1)=sigmau(ixco^s)+sigmad(ixco^s)
292  where(sigma(ixco^s,idim1)/=zero)
293  sigma(ixco^s,idim1)=abs(sigmau(ixco^s)-sigmad(ixco^s))/sigma(ixco^s,idim1)
294  elsewhere
295  sigma(ixco^s,idim1)=zero
296  end where
297 
298  end do
299  end do
300  end do
301 
302  ! Directional bias for Toth-Roe prolongation
303  do idim1=1,ndim
304  do idim2=1,ndim
305  do idim3=1,ndim
306  if (lvc(idim1,idim2,idim3)<1) cycle
307  ! Nonlinear
308  ! alpha(ixCo^S,idim1)=sigma(ixCo^S,idim2)-sigma(ixCo^S,idim3)
309  ! Homogeneous
310  ! alpha(ixCo^S,idim1)=0.d0
311  ! Toth-Roe
312  alpha(ixco^s,idim1)=(sco%dx(ixco^s,idim2)-sco%dx(ixco^s,idim3))/(sco%dx(ixco^s,idim2)+sco%dx(ixco^s,idim3))
313  end do
314  end do
315  end do
316  }
317 
318  do idim1=1,ndim
319  do idim2=1,ndim
320  {^iftwod
321  if (idim1==idim2) cycle
322  ixfiscmin^d=ixfismin^d+1;
323  ixfiscmax^d=ixfismax^d-1;
324  jxfisc^l=ixfisc^l+kr(idim1,^d);
325  hxfisc^l=ixfisc^l-kr(idim1,^d);
326  ipxfisc^l=ixfisc^l+kr(idim2,^d);
327  imxfisc^l=ixfisc^l-kr(idim2,^d);
328  jpxfisc^l=jxfisc^l+kr(idim2,^d);
329  jmxfisc^l=jxfisc^l-kr(idim2,^d);
330  hpxfisc^l=hxfisc^l+kr(idim2,^d);
331 
332  bfluxfi(ixfiscmin^d:ixfiscmax^d:2,idim1)=&
333  half*(bfluxfi(jxfiscmin^d:jxfiscmax^d:2,idim1)&
334  +bfluxfi(hxfiscmin^d:hxfiscmax^d:2,idim1))&
335  -quarter*(bfluxfi(ipxfiscmin^d:ipxfiscmax^d:2,idim2)&
336  -bfluxfi(jpxfiscmin^d:jpxfiscmax^d:2,idim2)&
337  -bfluxfi(imxfiscmin^d:imxfiscmax^d:2,idim2)&
338  +bfluxfi(jmxfiscmin^d:jmxfiscmax^d:2,idim2))
339 
340  bfluxfi(ipxfiscmin^d:ipxfiscmax^d:2,idim1)=&
341  half*(bfluxfi(jpxfiscmin^d:jpxfiscmax^d:2,idim1)&
342  +bfluxfi(hpxfiscmin^d:hpxfiscmax^d:2,idim1))&
343  -quarter*(bfluxfi(ipxfiscmin^d:ipxfiscmax^d:2,idim2)&
344  -bfluxfi(jpxfiscmin^d:jpxfiscmax^d:2,idim2)&
345  -bfluxfi(imxfiscmin^d:imxfiscmax^d:2,idim2)&
346  +bfluxfi(jmxfiscmin^d:jmxfiscmax^d:2,idim2))
347  }
348  {^ifthreed
349  do idim3=1,ndim
350  if (lvc(idim1,idim2,idim3)<1) cycle
351  ! Set up indices
352  hxfi^l=ixfi^l-kr(idim1,^d);
353  jxfi^l=ixfi^l+kr(idim1,^d);
354 
355  hijxfi^l=hxfi^l+kr(idim3,^d);
356  hjixfi^l=hxfi^l+kr(idim2,^d);
357  hjjxfi^l=hijxfi^l+kr(idim2,^d);
358 
359  iihxfi^l=ixfi^l-kr(idim3,^d);
360  iijxfi^l=ixfi^l+kr(idim3,^d);
361  ihixfi^l=ixfi^l-kr(idim2,^d);
362  ihjxfi^l=ihixfi^l+kr(idim3,^d);
363  ijixfi^l=ixfi^l+kr(idim2,^d);
364  ijhxfi^l=ijixfi^l-kr(idim3,^d);
365  ijjxfi^l=ijixfi^l+kr(idim3,^d);
366 
367  jihxfi^l=jxfi^l-kr(idim3,^d);
368  jijxfi^l=jxfi^l+kr(idim3,^d);
369  jhixfi^l=jxfi^l-kr(idim2,^d);
370  jhjxfi^l=jhixfi^l+kr(idim3,^d);
371  jjixfi^l=jxfi^l+kr(idim2,^d);
372  jjhxfi^l=jjixfi^l-kr(idim3,^d);
373  jjjxfi^l=jjixfi^l+kr(idim3,^d);
374 
375  ! Prolongation formulas
376  f1(ixco^s)=bfluxfi(ihixfimin^d:ihixfimax^d:2,idim2)&
377  -bfluxfi(jhixfimin^d:jhixfimax^d:2,idim2)&
378  -bfluxfi(ijixfimin^d:ijixfimax^d:2,idim2)&
379  +bfluxfi(jjixfimin^d:jjixfimax^d:2,idim2)
380 
381  f2(ixco^s)=bfluxfi(ihjxfimin^d:ihjxfimax^d:2,idim2)&
382  -bfluxfi(jhjxfimin^d:jhjxfimax^d:2,idim2)&
383  -bfluxfi(ijjxfimin^d:ijjxfimax^d:2,idim2)&
384  +bfluxfi(jjjxfimin^d:jjjxfimax^d:2,idim2)
385 
386  f3(ixco^s)=bfluxfi(iihxfimin^d:iihxfimax^d:2,idim3)&
387  -bfluxfi(jihxfimin^d:jihxfimax^d:2,idim3)&
388  -bfluxfi(iijxfimin^d:iijxfimax^d:2,idim3)&
389  +bfluxfi(jijxfimin^d:jijxfimax^d:2,idim3)
390 
391  f4(ixco^s)=bfluxfi(ijhxfimin^d:ijhxfimax^d:2,idim3)&
392  -bfluxfi(jjhxfimin^d:jjhxfimax^d:2,idim3)&
393  -bfluxfi(ijjxfimin^d:ijjxfimax^d:2,idim3)&
394  +bfluxfi(jjjxfimin^d:jjjxfimax^d:2,idim3)
395 
396  bfluxfi(ixfimin^d:ixfimax^d:2,idim1)=&
397  half*(bfluxfi(hxfimin^d:hxfimax^d:2,idim1)+bfluxfi(jxfimin^d:jxfimax^d:2,idim1))&
398  +6.25d-2*((3.d0+alpha(ixco^s,idim2))*f1(ixco^s)&
399  +(1.d0-alpha(ixco^s,idim2))*f2(ixco^s)&
400  +(3.d0-alpha(ixco^s,idim3))*f3(ixco^s)&
401  +(1.d0+alpha(ixco^s,idim3))*f4(ixco^s))
402 
403  bfluxfi(ijixfimin^d:ijixfimax^d:2,idim1)=&
404  half*(bfluxfi(hjixfimin^d:hjixfimax^d:2,idim1)+bfluxfi(jjixfimin^d:jjixfimax^d:2,idim1))&
405  +6.25d-2*((3.d0+alpha(ixco^s,idim2))*f1(ixco^s)&
406  +(1.d0-alpha(ixco^s,idim2))*f2(ixco^s)&
407  +(1.d0+alpha(ixco^s,idim3))*f3(ixco^s)&
408  +(3.d0-alpha(ixco^s,idim3))*f4(ixco^s))
409 
410  bfluxfi(iijxfimin^d:iijxfimax^d:2,idim1)=&
411  half*(bfluxfi(hijxfimin^d:hijxfimax^d:2,idim1)+bfluxfi(jijxfimin^d:jijxfimax^d:2,idim1))&
412  +6.25d-2*((1.d0-alpha(ixco^s,idim2))*f1(ixco^s)&
413  +(3.d0+alpha(ixco^s,idim2))*f2(ixco^s)&
414  +(3.d0-alpha(ixco^s,idim3))*f3(ixco^s)&
415  +(1.d0+alpha(ixco^s,idim3))*f4(ixco^s))
416 
417  bfluxfi(ijjxfimin^d:ijjxfimax^d:2,idim1)=&
418  half*(bfluxfi(hjjxfimin^d:hjjxfimax^d:2,idim1)+bfluxfi(jjjxfimin^d:jjjxfimax^d:2,idim1))&
419  +6.25d-2*((1.d0-alpha(ixco^s,idim2))*f1(ixco^s)&
420  +(3.d0+alpha(ixco^s,idim2))*f2(ixco^s)&
421  +(1.d0+alpha(ixco^s,idim3))*f3(ixco^s)&
422  +(3.d0-alpha(ixco^s,idim3))*f4(ixco^s))
423  end do
424  }
425  end do
426  end do
427 
428  ! Go back to magnetic fields
429  do idim1=1,ndim
430  ixfiscmax^d=ixfimax^d;
431  ixfiscmin^d=ixfimin^d-kr(^d,idim1);
432  where(sfi%surfaceC(ixfisc^s,idim1)/=zero)
433  wfis(ixfisc^s,idim1)=bfluxfi(ixfisc^s,idim1)/sfi%surfaceC(ixfisc^s,idim1)
434  elsewhere
435  wfis(ixfisc^s,idim1)=zero
436  end where
437  end do
438 
439  if(phys_total_energy.and. .not.prolongprimitive) then
440  b_energy_change(ixfi^s)=0.5d0*sum(wfi(ixfi^s,iw_mag(:))**2,dim=ndim+1)
441  end if
442  call phys_face_to_center(ixfi^l,sfi)
443  if(phys_total_energy.and. .not.prolongprimitive) then
444  b_energy_change(ixfi^s)=0.5d0*sum(wfi(ixfi^s,iw_mag(:))**2,dim=ndim+1)-&
445  b_energy_change(ixfi^s)
446  wfi(ixfi^s,iw_e)=wfi(ixfi^s,iw_e)+b_energy_change(ixfi^s)
447  end if
448 
449  end associate
450 
451  } ! END NOONED
452  end subroutine prolong_2nd_stg
453 
454  !> To achive consistency and thus conservation of divergence,
455  !> when refining a block we take into account the faces of the
456  !> already fine neighbours, if any. This routine stores them.
457  subroutine store_faces
458  use mod_forest, only: refine
460  integer :: igrid, iigrid, idims, iside, ineighbor, ipe_neighbor
461  integer :: nx^d, i^d, ic^d, inc^d
462 
463  if (npe>1) then
464  nsend_fc=0
465  nrecv_fc=0
466  end if
467 
468  ! Size of the block face
469  nx^d=ixmhi^d-ixmlo^d+1;
470 
471  do iigrid=1,igridstail; igrid=igrids(iigrid);
472  ! Check whether it is necessary to store any block face, i.e.
473  ! if any coarser neighbour is going to be refined
474  {do iside=1,2
475  i^dd=kr(^dd,^d)*(2*iside-3);
476  if (neighbor_pole(i^dd,igrid)/=0) cycle
477  if (neighbor_type(i^dd,igrid)==neighbor_coarse) then
478  ineighbor =neighbor(1,i^dd,igrid)
479  ipe_neighbor=neighbor(2,i^dd,igrid)
480  if (refine(ineighbor,ipe_neighbor)) then
481  allocate(pface(iside,^d,igrid)%face(1^d%1:nx^dd))
482  !! Store the faces
483  if (iside==1) then !! left side
484  pface(iside,^d,igrid)%face(1^d%1:nx^dd)=&
485  ps(igrid)%ws(ixmlo^d-1^d%ixM^t,^d)
486  else !! right side
487  pface(iside,^d,igrid)%face(1^d%1:nx^dd)=&
488  ps(igrid)%ws(ixmhi^d^d%ixM^t,^d)
489  end if
490  if (ipe_neighbor/=mype) nsend_fc(^d)=nsend_fc(^d)+1
491  end if
492  end if
493  end do\}
494 
495  ! If a grid is going to be refined,
496  ! remember what are its neighbours.
497  if (refine(igrid,mype)) then
498  ! Initialize neighbour array
499  fine_neighbors(:^d&,:,igrid)%igrid=-1
500  fine_neighbors(:^d&,:,igrid)%ipe=-1
501  do idims=1,ndim
502  do iside=1,2
503  i^d=kr(^d,idims)*(2*iside-3);
504  if (neighbor_pole(i^d,igrid)/=0) cycle
505  if (neighbor_type(i^d,igrid)==neighbor_fine) then
506  {do ic^db=1+int((1+i^db)/2),2-int((1-i^db)/2)
507  inc^db=ic^db+i^db\}
508  ineighbor=neighbor_child(1,inc^d,igrid)
509  ipe_neighbor=neighbor_child(2,inc^d,igrid)
510 
511  fine_neighbors(ic^d,idims,igrid)%igrid= ineighbor
512  fine_neighbors(ic^d,idims,igrid)%ipe=ipe_neighbor
513 
514  if (ipe_neighbor/=mype) nrecv_fc(idims)=nrecv_fc(idims)+1
515  {end do\}
516  end if
517  end do
518  end do
519  end if
520 
521  end do
522 
523  end subroutine store_faces
524 
525  !> When refining a coarse block with fine neighbours, it is necessary
526  !> prolong consistently with the already fine faces.
527  !> This routine takes care of the communication of such faces.
528  subroutine comm_faces
529  use mod_forest, only: refine
531 
532  integer :: iigrid,igrid,ineighbor,ipe_neighbor
533  integer :: idims,iside,i^d,ic^d,inc^d,nx^d
534  integer :: recvsize, sendsize
535 
536  ! Communicate the block faces to achieve consistency when refining
537  ! Initialize communication structures
538 
539  nrecv=0
540  nsend=0
541  recvsize=0
542  sendsize=0
543 
544  do idims=1,ndim
545  select case (idims)
546  {case (^d)
547  nrecv=nrecv+nrecv_fc(^d)
548  nsend=nsend+nsend_fc(^d)
549  nx^d=1^d%nx^dd=ixmhi^dd-ixmlo^dd+1;
550  isize(^d)={nx^dd*}
551  recvsize=recvsize+nrecv_fc(^d)*isize(^d)
552  sendsize=sendsize+nsend_fc(^d)*isize(^d)
553  \}
554  end select
555  end do
556 
557  if (nrecv>0) then
558  ! Allocate receive buffer
559  allocate(recvbuffer(recvsize),recvstatus(mpi_status_size,nrecv), &
560  recvrequest(nrecv))
561  recvrequest=mpi_request_null
562  ibuf_recv=1
563  irecv=0
564 
565  ! Receive
566  do iigrid=1,igridstail; igrid=igrids(iigrid);
567  if (refine(igrid,mype)) then
568  {do ic^db=1,2\}
569  ! Only one of the sides will be necessary,
570  ! so we do the loop only over dimensions, instead of
571  ! over dimensions and sizes as in the routines
572  ! old_neighbors and already_fine.
573  do idims=1,ndim
574  ipe_neighbor=fine_neighbors(ic^d,idims,igrid)%ipe
575  ineighbor =fine_neighbors(ic^d,idims,igrid)%igrid
576  if (ineighbor>0.and.ipe_neighbor/=mype) then
577  {if (idims==^d) iside=ic^d\}
578  !!! Check indices
579  i^d=kr(^d,idims)*(2*iside-3);
580  if (neighbor_pole(i^d,igrid)/=0) cycle
581  inc^d=ic^d+i^d;
582  irecv=irecv+1
583  itag=4**^nd*(igrid-1)+{inc^d*4**(^d-1)+}
584 
585  call mpi_irecv(recvbuffer(ibuf_recv),isize(idims), &
586  mpi_double_precision,ipe_neighbor,itag, &
587  icomm,recvrequest(irecv),ierrmpi)
588  ibuf_recv=ibuf_recv+isize(idims)
589  end if
590  end do
591  {end do\}
592  end if
593  end do
594 
595  end if
596 
597  if (nsend>0) then
598  ! Allocate send buffer
599  allocate(sendbuffer(sendsize),sendstatus(mpi_status_size,nsend),sendrequest(nsend))
600  sendrequest=mpi_request_null
601  isend=0
602  ibuf_send=1
603  do iigrid=1,igridstail; igrid=igrids(iigrid);
604  ! Check whether it is necessary to store any block face, i.e.
605  ! if any coarser neighbour is going to be refined
606  {do iside=1,2
607  i^dd=kr(^dd,^d)*(2*iside-3);
608  ! When there is a pole, faces are always zero and this is not necessary
609  if (neighbor_pole(i^dd,igrid)/=0) cycle
610  if (neighbor_type(i^dd,igrid)==neighbor_coarse) then
611  ineighbor =neighbor(1,i^dd,igrid)
612  ipe_neighbor=neighbor(2,i^dd,igrid)
613  if (refine(ineighbor,ipe_neighbor)) then
614  if (ipe_neighbor/=mype) then
615  ic^dd=1+modulo(node(pig^dd_,igrid)-1,2);
616  inc^d=-2*i^d+ic^d^d%inc^dd=ic^dd;
617  itag=4**^nd*(ineighbor-1)+{inc^dd*4**(^dd-1)+}
618  isend=isend+1
619  ibuf_send_next=ibuf_send+isize(^d)
620  sendbuffer(ibuf_send:ibuf_send_next-1)=&
621  reshape(pface(iside,^d,igrid)%face,(/isize(^d)/))
622  call mpi_isend(sendbuffer(ibuf_send),isize(^d), &
623  mpi_double_precision,ipe_neighbor,itag, &
624  icomm,sendrequest(isend),ierrmpi)
625  ibuf_send=ibuf_send_next
626  end if
627  end if
628  end if
629  end do\}
630  end do
631  end if
632 
633  ! Waitalls
634  if (nrecv>0) then
635  call mpi_waitall(nrecv,recvrequest,recvstatus,ierrmpi)
636  deallocate(recvstatus,recvrequest)
637  ibuf_recv=1
638  end if
639 
640  if (nsend>0) then
641  call mpi_waitall(nsend,sendrequest,sendstatus,ierrmpi)
642  deallocate(sendbuffer,sendstatus,sendrequest)
643  end if
644 
645  end subroutine comm_faces
646 
647  subroutine end_comm_faces
649  ! Deallocate receive buffer
650  if (nrecv>0) deallocate(recvbuffer)
651  end subroutine end_comm_faces
652 
653  subroutine deallocatebfaces
655  integer :: igrid, iigrid, iside
656 
657  do iigrid=1,igridstail; igrid=igrids(iigrid);
658  {do iside=1,2
659  if (associated(pface(iside,^d,igrid)%face)) then
660  deallocate(pface(iside,^d,igrid)%face)
661  end if
662  end do\}
663  end do
664 
665  end subroutine deallocatebfaces
666 
667  subroutine old_neighbors(child_igrid,child_ipe,igrid,ipe)
669  integer, dimension(2^D&), intent(in) :: child_igrid, child_ipe
670  integer, intent(in) :: igrid, ipe
671  integer :: iside, i^d, ic^d
672 
673  {do ic^db=1,2\}
674  old_neighbor(:,:^d&,child_igrid(ic^d))=-1
675  {do iside=1,2
676  if (ic^d==iside) then
677  i^dd=kr(^dd,^d)*(2*iside-3);
678  old_neighbor(1,i^dd,child_igrid(ic^dd))=&
679  fine_neighbors(ic^dd,^d,igrid)%igrid
680  old_neighbor(2,i^dd,child_igrid(ic^dd))=&
681  fine_neighbors(ic^dd,^d,igrid)%ipe
682  end if
683  end do\}
684  {end do\}
685 
686  end subroutine old_neighbors
687 
688  !> This routine fills the fine faces before prolonging.
689  !> It is the face equivalent of fix_conserve
690  subroutine already_fine(sFi,ichild,fine_^L)
691  use mod_forest
693  type(tree_node_ptr) :: tree
694  type(state) :: sfi
695  integer, intent(in) :: ichild
696  logical :: fine_^l
697 
698  integer :: ineighbor,ipe_neighbor,ibufnext
699  integer :: iside,iotherside,i^d,nx^d
700 
701  ! Size of the block face
702  nx^d=ixmhi^d-ixmlo^d+1;
703 
704  ! Initialise everything to zero and false
705  fine_min^d=.false.;
706  fine_max^d=.false.;
707  sfi%ws=zero
708 
709  {^nooned
710  ! This face communication is not needed in 1D
711  {do iside=1,2
712  i^dd=kr(^dd,^d)*(2*iside-3);
713  ! This is not necessary at the pole.
714  ! We are inside a loop over the children, so the grid index is ichild
715  if (neighbor_pole(i^dd,ichild)/=0) cycle
716 
717  ! Get old ipe and igrid of neighbour from array fake_neighbor
718  ! Then plug it into the structure pfaces and get the faces
719  ineighbor =old_neighbor(1,i^dd,ichild)
720  ipe_neighbor=old_neighbor(2,i^dd,ichild)
721 
722  iotherside=3-iside
723  if (ineighbor>0) then
724  if (iside==1) then ! left side
725  if (ipe_neighbor==mype) then
726  sfi%ws(ixmlo^d-1^d%ixM^t,^d)=pface(iotherside,^d,ineighbor)%face(1^d%1:nx^dd)
727  else
728  ibufnext=ibuf_recv+isize(^d)
729  sfi%ws(ixmlo^d-1^d%ixM^t,^d)=reshape(&
730  source=recvbuffer(ibuf_recv:ibufnext-1),&
731  shape=shape(sfi%ws(ixmlo^d-1^d%ixM^t,^d)))
732  ibuf_recv=ibufnext
733  end if
734  fine_min^d=.true.
735  else ! right side
736  if (ipe_neighbor==mype) then
737  sfi%ws(ixmhi^d^d%ixM^t,^d)=pface(iotherside,^d,ineighbor)%face(1^d%1:nx^dd)
738  else
739  ibufnext=ibuf_recv+isize(^d)
740  sfi%ws(ixmhi^d^d%ixM^t,^d)=reshape(&
741  source=recvbuffer(ibuf_recv:ibufnext-1),&
742  shape=shape(sfi%ws(ixmhi^d^d%ixM^t,^d)))
743  ibuf_recv=ibufnext
744  end if
745  fine_max^d=.true.
746  end if
747  end if
748  end do\}
749  }
750  end subroutine already_fine
751 
752 end module mod_amr_fct
type(fake_neighbors), dimension(:^d &,:,:), allocatable, public fine_neighbors
Definition: mod_amr_fct.t:16
subroutine, public old_neighbors(child_igrid, child_ipe, igrid, ipe)
Definition: mod_amr_fct.t:668
integer, dimension(:,:^d &,:), allocatable, public old_neighbor
Definition: mod_amr_fct.t:18
subroutine, public store_faces
To achive consistency and thus conservation of divergence, when refining a block we take into account...
Definition: mod_amr_fct.t:458
subroutine, public comm_faces
When refining a coarse block with fine neighbours, it is necessary prolong consistently with the alre...
Definition: mod_amr_fct.t:529
subroutine, public already_fine(sFi, ichild, fine_L)
This routine fills the fine faces before prolonging. It is the face equivalent of fix_conserve.
Definition: mod_amr_fct.t:691
type(facealloc), dimension(:,:,:), allocatable, public pface
Definition: mod_amr_fct.t:14
subroutine, public end_comm_faces
Definition: mod_amr_fct.t:648
subroutine, public prolong_2nd_stg(sCo, sFi, ixCoLin, ixFiLin, dxCoD, xCominD, dxFiD, xFiminD, ghost, fine_Lin)
This subroutine performs a 2nd order prolongation for a staggered field F, preserving the divergence ...
Definition: mod_amr_fct.t:41
subroutine, public deallocatebfaces
Definition: mod_amr_fct.t:654
Module with basic grid data structures.
Definition: mod_forest.t:2
logical, dimension(:,:), allocatable, save refine
Definition: mod_forest.t:70
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, parameter ndim
Number of spatial dimensions for grid variables.
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
integer, dimension(:), allocatable, parameter d
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
integer npe
The number of MPI tasks.
integer nghostcells
Number of ghost cells surrounding a grid.
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition: mod_physics.t:4
Pointer to a tree_node.
Definition: mod_forest.t:6