MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code
mod_fix_conserve.t
Go to the documentation of this file.
1 !> Module for flux conservation near refinement boundaries
3  implicit none
4  private
5 
6  type fluxalloc
7  double precision, dimension(:^D&,:), pointer:: flux => null()
8  double precision, dimension(:^D&,:), pointer:: edge => null()
9  end type fluxalloc
10  !> store flux to fix conservation
11  type(fluxalloc), dimension(:,:,:), allocatable, public :: pflux
12 
13  integer, save :: nrecv, nsend
14  double precision, allocatable, save :: recvbuffer(:), sendbuffer(:)
15  integer, dimension(:), allocatable :: fc_recvreq, fc_sendreq
16  integer, dimension(:,:), allocatable :: fc_recvstat, fc_sendstat
17  integer, dimension(^ND), save :: isize
18  integer :: ibuf, ibuf_send
19  ! ct for corner total
20  integer, save :: nrecv_ct, nsend_ct
21  ! buffer for corner coarse
22  double precision, allocatable, save :: recvbuffer_cc(:), sendbuffer_cc(:)
23  integer, dimension(:), allocatable :: cc_recvreq, cc_sendreq
24  integer, dimension(:,:), allocatable :: cc_recvstat, cc_sendstat
25  integer, dimension(^ND), save :: isize_stg
26  integer :: ibuf_cc, ibuf_cc_send
27  integer :: itag, itag_cc, isend, isend_cc, irecv, irecv_cc
28 
29  public :: init_comm_fix_conserve
30  public :: allocatebflux
31  public :: deallocatebflux
32  public :: sendflux
33  public :: recvflux
34  public :: store_flux
35  public :: store_edge
36  public :: fix_conserve
37  public :: fix_edges
38 
39  contains
40 
41  subroutine init_comm_fix_conserve(idim^LIM,nwfluxin)
43 
44  integer, intent(in) :: idim^lim,nwfluxin
45 
46  integer :: iigrid, igrid, idims, iside, i^d, nxco^d
47  integer :: ic^d, inc^d, ipe_neighbor
48  integer :: recvsize, sendsize
49  integer :: recvsize_cc, sendsize_cc
50 
51  nsend = 0
52  nrecv = 0
53  recvsize = 0
54  sendsize = 0
55  if(stagger_grid) then
56  ! Special communication for diagonal 'coarse corners'
57  ! nrecv/send_cc (for 'coarse corners' is a dim=ndim-1 array which
58  ! stores the faces that must be communicated in each direction.
59  ! nrecv/send_ct (for 'corners total' is the total number of
60  ! necessary communications. These special cases have their own
61  ! send and receive buffers (send/recvbuffer_cc), their tags, etc.
62  nsend_ct=0
63  nrecv_ct=0
64  recvsize_cc=0
65  sendsize_cc=0
66  end if
67 
68  do idims= idim^lim
69  select case (idims)
70  {case (^d)
71  nrecv=nrecv+nrecv_fc(^d)
72  nsend=nsend+nsend_fc(^d)
73  nxco^d=1^d%nxCo^dd=ixghi^dd/2-nghostcells;
74  isize(^d)={nxco^dd*}*(nwfluxin)
75  recvsize=recvsize+nrecv_fc(^d)*isize(^d)
76  sendsize=sendsize+nsend_fc(^d)*isize(^d)
77  if(stagger_grid) then
78  ! This does not consider the 'coarse corner' case
79  nxco^d=1^d%nxCo^dd=ixghi^dd/2-nghostcells+1;
80  isize_stg(^d)={nxco^dd*}*(^nd-1)
81  ! the whole size is used (cell centered and staggered)
82  isize(^d)=isize(^d)+isize_stg(^d)
83  recvsize=recvsize+nrecv_fc(^d)*isize_stg(^d)
84  sendsize=sendsize+nsend_fc(^d)*isize_stg(^d)
85  ! Coarse corner case
86  nrecv_ct=nrecv_ct+nrecv_cc(^d)
87  nsend_ct=nsend_ct+nsend_cc(^d)
88  recvsize_cc=recvsize_cc+nrecv_cc(^d)*isize_stg(^d)
89  sendsize_cc=sendsize_cc+nsend_cc(^d)*isize_stg(^d)
90  end if
91  \}
92  end select
93  end do
94 
95  ! Reallocate buffers when size differs
96  if (allocated(recvbuffer)) then
97  if (recvsize /= size(recvbuffer)) then
98  deallocate(recvbuffer)
99  allocate(recvbuffer(recvsize))
100  end if
101  else
102  allocate(recvbuffer(recvsize))
103  end if
104 
105  if (allocated(fc_recvreq)) then
106  if (nrecv /= size(fc_recvreq)) then
107  deallocate(fc_recvreq, fc_recvstat)
108  allocate(fc_recvstat(mpi_status_size,nrecv), fc_recvreq(nrecv))
109  end if
110  else
111  allocate(fc_recvstat(mpi_status_size,nrecv), fc_recvreq(nrecv))
112  end if
113 
114  if (allocated(fc_sendreq)) then
115  if (nsend /= size(fc_sendreq)) then
116  deallocate(fc_sendreq, fc_sendstat)
117  allocate(fc_sendstat(mpi_status_size,nsend), fc_sendreq(nsend))
118  end if
119  else
120  allocate(fc_sendstat(mpi_status_size,nsend), fc_sendreq(nsend))
121  end if
122 
123  if(stagger_grid) then
124 
125  if (allocated(sendbuffer)) then
126  if (sendsize /= size(sendbuffer)) then
127  deallocate(sendbuffer)
128  allocate(sendbuffer(sendsize))
129  end if
130  else
131  allocate(sendbuffer(sendsize))
132  end if
133 
134  if (allocated(recvbuffer_cc)) then
135  if (recvsize_cc /= size(recvbuffer_cc)) then
136  deallocate(recvbuffer_cc)
137  allocate(recvbuffer_cc(recvsize_cc))
138  end if
139  else
140  allocate(recvbuffer_cc(recvsize_cc))
141  end if
142 
143  if (allocated(cc_recvreq)) then
144  if (nrecv_ct /= size(cc_recvreq)) then
145  deallocate(cc_recvreq, cc_recvstat)
146  allocate(cc_recvstat(mpi_status_size,nrecv_ct), cc_recvreq(nrecv_ct))
147  end if
148  else
149  allocate(cc_recvstat(mpi_status_size,nrecv_ct), cc_recvreq(nrecv_ct))
150  end if
151 
152  if (allocated(sendbuffer_cc)) then
153  if (sendsize_cc /= size(sendbuffer_cc)) then
154  deallocate(sendbuffer_cc)
155  allocate(sendbuffer_cc(sendsize_cc))
156  end if
157  else
158  allocate(sendbuffer_cc(sendsize_cc))
159  end if
160 
161  if (allocated(cc_sendreq)) then
162  if (nsend_ct /= size(cc_sendreq)) then
163  deallocate(cc_sendreq, cc_sendstat)
164  allocate(cc_sendstat(mpi_status_size,nsend_ct), cc_sendreq(nsend_ct))
165  end if
166  else
167  allocate(cc_sendstat(mpi_status_size,nsend_ct), cc_sendreq(nsend_ct))
168  end if
169  end if
170 
171  end subroutine init_comm_fix_conserve
172 
173  subroutine recvflux(idim^LIM)
175 
176  integer, intent(in) :: idim^lim
177 
178  integer :: iigrid, igrid, idims, iside, i^d, nxco^d
179  integer :: ic^d, inc^d, ipe_neighbor
180  integer :: pi^d,mi^d,ph^d,mh^d,idir
181 
182  if (nrecv>0) then
183  fc_recvreq=mpi_request_null
184  ibuf=1
185  irecv=0
186 
187  do iigrid=1,igridstail; igrid=igrids(iigrid);
188  do idims= idim^lim
189  do iside=1,2
190  i^d=kr(^d,idims)*(2*iside-3);
191 
192  if (neighbor_pole(i^d,igrid)/=0) cycle
193 
194  if (neighbor_type(i^d,igrid)/=4) cycle
195  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
196  inc^db=2*i^db+ic^db\}
197  ipe_neighbor=neighbor_child(2,inc^d,igrid)
198  if (ipe_neighbor/=mype) then
199  irecv=irecv+1
200  itag=4**^nd*(igrid-1)+{inc^d*4**(^d-1)+}
201  call mpi_irecv(recvbuffer(ibuf),isize(idims), &
202  mpi_double_precision,ipe_neighbor,itag, &
203  icomm,fc_recvreq(irecv),ierrmpi)
204  ibuf=ibuf+isize(idims)
205  end if
206  {end do\}
207  end do
208  end do
209  end do
210  end if
211 
212  if(stagger_grid) then
213  ! receive corners
214  if (nrecv_ct>0) then
215  cc_recvreq=mpi_request_null
216  ibuf_cc=1
217  irecv_cc=0
218 
219  do iigrid=1,igridstail; igrid=igrids(iigrid);
220  do idims= idim^lim
221  do iside=1,2
222  i^d=kr(^d,idims)*(2*iside-3);
223  ! Check if there are special corners
224  ! (Coarse block diagonal to a fine block)
225  ! If there are, receive.
226  ! Tags are calculated in the same way as for
227  ! normal fluxes, but should not overlap because
228  ! inc^D are different
229  if (neighbor_type(i^d,igrid)==3) then
230  do idir=idims+1,ndim
231  pi^d=i^d+kr(idir,^d);
232  mi^d=i^d-kr(idir,^d);
233  ph^d=pi^d-kr(idims,^d)*(2*iside-3);
234  mh^d=mi^d-kr(idims,^d)*(2*iside-3);
235 
236  if (neighbor_type(pi^d,igrid)==4.and.&
237  neighbor_type(ph^d,igrid)==3.and.&
238  neighbor_pole(pi^d,igrid)==0) then
239  ! Loop on children (several in 3D)
240  {do ic^db=1+int((1-pi^db)/2),2-int((1+pi^db)/2)
241  inc^db=2*pi^db+ic^db\}
242  ipe_neighbor=neighbor_child(2,inc^d,igrid)
243  if (mype/=ipe_neighbor) then
244  irecv_cc=irecv_cc+1
245  itag_cc=4**^nd*(igrid-1)+{inc^d*4**(^d-1)+}
246  call mpi_irecv(recvbuffer_cc(ibuf_cc),isize_stg(idims),&
247  mpi_double_precision,ipe_neighbor,itag_cc,&
248  icomm,cc_recvreq(irecv_cc),ierrmpi)
249  ibuf_cc=ibuf_cc+isize_stg(idims)
250  end if
251  {end do\}
252  end if
253 
254  if (neighbor_type(mi^d,igrid)==4.and.&
255  neighbor_type(mh^d,igrid)==3.and.&
256  neighbor_pole(mi^d,igrid)==0) then
257  ! Loop on children (several in 3D)
258  {do ic^db=1+int((1-mi^db)/2),2-int((1+mi^db)/2)
259  inc^db=2*mi^db+ic^db\}
260  ipe_neighbor=neighbor_child(2,inc^d,igrid)
261  if (mype/=ipe_neighbor) then
262  irecv_cc=irecv_cc+1
263  itag_cc=4**^nd*(igrid-1)+{inc^d*4**(^d-1)+}
264  call mpi_irecv(recvbuffer_cc(ibuf_cc),isize_stg(idims),&
265  mpi_double_precision,ipe_neighbor,itag_cc,&
266  icomm,cc_recvreq(irecv_cc),ierrmpi)
267  ibuf_cc=ibuf_cc+isize_stg(idims)
268  end if
269  {end do\}
270  end if
271  end do
272  end if
273  end do
274  end do
275  end do
276  end if
277  end if ! end if stagger grid
278 
279  end subroutine recvflux
280 
281  subroutine sendflux(idim^LIM)
283 
284  integer, intent(in) :: idim^lim
285 
286  integer :: idims, iside, i^d, ic^d, inc^d, ix^d, ixco^d, nxco^d, iw
287  integer :: ineighbor, ipe_neighbor, igrid, iigrid, ibuf_send_next
288  integer :: idir, ibuf_cc_send_next, pi^d, ph^d, mi^d, mh^d
289 
290  fc_sendreq = mpi_request_null
291  isend = 0
292  if(stagger_grid) then
293  ibuf_send = 1
294  cc_sendreq=mpi_request_null
295  isend_cc=0
296  ibuf_cc_send=1
297  end if
298 
299  do iigrid=1,igridstail; igrid=igrids(iigrid);
300  do idims = idim^lim
301  select case (idims)
302  {case (^d)
303  do iside=1,2
304  i^dd=kr(^dd,^d)*(2*iside-3);
305 
306  if (neighbor_pole(i^dd,igrid)/=0) cycle
307 
308  if (neighbor_type(i^dd,igrid)==neighbor_coarse) then
309  ! send flux to coarser neighbor
310  ineighbor=neighbor(1,i^dd,igrid)
311  ipe_neighbor=neighbor(2,i^dd,igrid)
312  if (ipe_neighbor/=mype) then
313  ic^dd=1+modulo(node(pig^dd_,igrid)-1,2);
314  inc^d=-2*i^d+ic^d^d%inc^dd=ic^dd;
315  itag=4**^nd*(ineighbor-1)+{inc^dd*4**(^dd-1)+}
316  isend=isend+1
317 
318  if(stagger_grid) then
319  ibuf_send_next=ibuf_send+isize(^d)
320  sendbuffer(ibuf_send:ibuf_send_next-isize_stg(^d)-1)=&
321  reshape(pflux(iside,^d,igrid)%flux,(/isize(^d)-isize_stg(^d)/))
322 
323  sendbuffer(ibuf_send_next-isize_stg(^d):ibuf_send_next-1)=&
324  reshape(pflux(iside,^d,igrid)%edge,(/isize_stg(^d)/))
325  call mpi_isend(sendbuffer(ibuf_send),isize(^d), &
326  mpi_double_precision,ipe_neighbor,itag, &
327  icomm,fc_sendreq(isend),ierrmpi)
328  ibuf_send=ibuf_send_next
329  else
330  call mpi_isend(pflux(iside,^d,igrid)%flux,isize(^d), &
331  mpi_double_precision,ipe_neighbor,itag, &
332  icomm,fc_sendreq(isend),ierrmpi)
333  end if
334  end if
335 
336  if(stagger_grid) then
337  ! If we are in a fine block surrounded by coarse blocks
338  do idir=idims+1,ndim
339  pi^dd=i^dd+kr(idir,^dd);
340  mi^dd=i^dd-kr(idir,^dd);
341  ph^dd=pi^dd-kr(idims,^dd)*(2*iside-3);
342  mh^dd=mi^dd-kr(idims,^dd)*(2*iside-3);
343 
344  if (neighbor_type(pi^dd,igrid)==2.and.&
345  neighbor_type(ph^dd,igrid)==2.and.&
346  mype/=neighbor(2,pi^dd,igrid).and.&
347  neighbor_pole(pi^dd,igrid)==0) then
348  ! Get relative position in the grid for tags
349  ineighbor=neighbor(1,pi^dd,igrid)
350  ipe_neighbor=neighbor(2,pi^dd,igrid)
351  ic^dd=1+modulo(node(pig^dd_,igrid)-1,2);
352  inc^dd=-2*pi^dd+ic^dd;
353  itag_cc=4**^nd*(ineighbor-1)+{inc^dd*4**(^dd-1)+}
354  ! Reshape to buffer and send
355  isend_cc=isend_cc+1
356  ibuf_cc_send_next=ibuf_cc_send+isize_stg(^d)
357  sendbuffer_cc(ibuf_cc_send:ibuf_cc_send_next-1)=&
358  reshape(pflux(iside,^d,igrid)%edge,shape=(/isize_stg(^d)/))
359  call mpi_isend(sendbuffer_cc(ibuf_cc_send),isize_stg(^d),&
360  mpi_double_precision,ipe_neighbor,itag_cc,&
361  icomm,cc_sendreq(isend_cc),ierrmpi)
362  ibuf_cc_send=ibuf_cc_send_next
363  end if
364 
365  if (neighbor_type(mi^dd,igrid)==2.and.&
366  neighbor_type(mh^dd,igrid)==2.and.&
367  mype/=neighbor(2,mi^dd,igrid).and.&
368  neighbor_pole(mi^dd,igrid)==0) then
369  ! Get relative position in the grid for tags
370  ineighbor=neighbor(1,mi^dd,igrid)
371  ipe_neighbor=neighbor(2,mi^dd,igrid)
372  ic^dd=1+modulo(node(pig^dd_,igrid)-1,2);
373  inc^dd=-2*pi^dd+ic^dd;
374  inc^dd=-2*mi^dd+ic^dd;
375  itag_cc=4**^nd*(ineighbor-1)+{inc^dd*4**(^dd-1)+}
376  ! Reshape to buffer and send
377  isend_cc=isend_cc+1
378  ibuf_cc_send_next=ibuf_cc_send+isize_stg(^d)
379  sendbuffer_cc(ibuf_cc_send:ibuf_cc_send_next-1)=&
380  reshape(pflux(iside,^d,igrid)%edge,shape=(/isize_stg(^d)/))
381  call mpi_isend(sendbuffer_cc(ibuf_cc_send),isize_stg(^d),&
382  mpi_double_precision,ipe_neighbor,itag_cc,&
383  icomm,cc_sendreq(isend_cc),ierrmpi)
384  ibuf_cc_send=ibuf_cc_send_next
385  end if
386  end do
387  end if ! end if stagger grid
388 
389  end if
390  end do\}
391  end select
392  end do
393  end do
394  end subroutine sendflux
395 
396  subroutine allocatebflux
398 
399  integer :: iigrid, igrid, iside, i^d, nx^d, nxco^d
400  integer :: idir,idim,pi^d, mi^d, ph^d, mh^d ! To detect corners
401 
402  nx^d=ixmhi^d-ixmlo^d+1;
403  nxco^d=nx^d/2;
404 
405  do iigrid=1,igridstail; igrid=igrids(iigrid);
406  ! For every grid,
407  ! arrays for the fluxes are allocated for every face direction(^D)
408  ! and every side (1=left, 2=right)
409  {do iside=1,2
410  i^dd=kr(^dd,^d)*(2*iside-3);
411 
412  if (neighbor_pole(i^dd,igrid)/=0) cycle
413 
414  select case (neighbor_type(i^dd,igrid))
415  case(neighbor_fine)
416  allocate(pflux(iside,^d,igrid)%flux(1^d%1:nx^dd,1:nwflux))
417  if(stagger_grid) allocate(pflux(iside,^d,igrid)%edge(1^d%0:nx^dd,1:ndim-1))
418  case(neighbor_coarse)
419  allocate(pflux(iside,^d,igrid)%flux(1^d%1:nxco^dd,1:nwflux))
420  if(stagger_grid) allocate(pflux(iside,^d,igrid)%edge(1^d%0:nxco^dd,1:ndim-1))
421  case(neighbor_sibling)
422  if(stagger_grid) then
423  idim=^d
424  do idir=idim+1,ndim
425  !do idir=min(idim+1,ndim),ndim
426  pi^dd=i^dd+kr(idir,^dd);
427  mi^dd=i^dd-kr(idir,^dd);
428  ph^dd=pi^dd-kr(^d,^dd)*(2*iside-3);
429  mh^dd=mi^dd-kr(^d,^dd)*(2*iside-3);
430  if ((neighbor_type(pi^dd,igrid)==4&
431  .and.neighbor_type(ph^dd,igrid)==3)&
432  .or.(neighbor_type(mi^dd,igrid)==4&
433  .and.neighbor_type(mh^dd,igrid)==3)) then
434  allocate(pflux(iside,^d,igrid)%edge(1^d%0:nx^dd,1:ndim-1))
435  exit
436  end if
437  end do
438  end if
439  end select
440  end do\}
441  end do
442 
443  end subroutine allocatebflux
444 
445  subroutine deallocatebflux
447 
448  integer :: iigrid, igrid, iside
449 
450  do iigrid=1,igridstail; igrid=igrids(iigrid);
451  {do iside=1,2
452  if (associated(pflux(iside,^d,igrid)%flux)) then
453  deallocate(pflux(iside,^d,igrid)%flux)
454  nullify(pflux(iside,^d,igrid)%flux)
455  end if
456  if (associated(pflux(iside,^d,igrid)%edge)) then
457  deallocate(pflux(iside,^d,igrid)%edge)
458  nullify(pflux(iside,^d,igrid)%edge)
459  end if
460  end do\}
461  end do
462 
463  end subroutine deallocatebflux
464 
465  subroutine fix_conserve(psb,idim^LIM,nw0,nwfluxin)
467 
468  integer, intent(in) :: idim^lim, nw0, nwfluxin
469  type(state) :: psb(max_blocks)
470 
471  integer :: iigrid, igrid, idims, iside, iotherside, i^d, ic^d, inc^d, ix^l
472  integer :: nxco^d, iw, ix, ipe_neighbor, ineighbor, nbuf, ibufnext, nw1
473  double precision :: cofiratio
474 
475  nw1=nw0-1+nwfluxin
476  if (slab_uniform) then
477  ! The flux is divided by volume of fine cell. We need, however,
478  ! to divide by volume of coarse cell => muliply by volume ratio
479  cofiratio=one/dble(2**ndim)
480  end if
481 
482  if (nrecv>0) then
483  call mpi_waitall(nrecv,fc_recvreq,fc_recvstat,ierrmpi)
484  ibuf=1
485  end if
486 
487  nxco^d=(ixmhi^d-ixmlo^d+1)/2;
488 
489  ! for all grids: perform flux update at Coarse-Fine interfaces
490  do iigrid=1,igridstail; igrid=igrids(iigrid);
491  do idims= idim^lim
492  select case (idims)
493  {case (^d)
494  do iside=1,2
495  i^dd=kr(^dd,^d)*(2*iside-3);
496 
497  if (neighbor_pole(i^dd,igrid)/=0) cycle
498 
499  if (neighbor_type(i^dd,igrid)/=4) cycle
500 
501  ! opedit: skip over active/passive interface since flux for passive ones is
502  ! not computed, keep the buffer counter up to date:
503  if (.not.neighbor_active(i^dd,igrid).or.&
504  .not.neighbor_active(0^dd&,igrid) ) then
505  {do ic^ddb=1+int((1-i^ddb)/2),2-int((1+i^ddb)/2)
506  inc^ddb=2*i^ddb+ic^ddb\}
507  ipe_neighbor=neighbor_child(2,inc^dd,igrid)
508  if (ipe_neighbor/=mype) then
509  ibufnext=ibuf+isize(^d)
510  ibuf=ibufnext
511  end if
512  {end do\}
513  cycle
514  end if
515  !
516 
517  select case (iside)
518  case (1)
519  ix=ixmlo^d
520  case (2)
521  ix=ixmhi^d
522  end select
523 
524  ! remove coarse flux
525  if (slab_uniform) then
526  psb(igrid)%w(ix^d%ixM^t,nw0:nw1) &
527  = psb(igrid)%w(ix^d%ixM^t,nw0:nw1) &
528  -pflux(iside,^d,igrid)%flux(1^d%:^dd&,1:nwfluxin)
529  else
530  do iw=nw0,nw1
531  psb(igrid)%w(ix^d%ixM^t,iw)=psb(igrid)%w(ix^d%ixM^t,iw)&
532  -pflux(iside,^d,igrid)%flux(1^d%:^dd&,iw-nw0+1) &
533  /ps(igrid)%dvolume(ix^d%ixM^t)
534  end do
535  end if
536 
537 
538  ! add fine flux
539  {do ic^ddb=1+int((1-i^ddb)/2),2-int((1+i^ddb)/2)
540  inc^ddb=2*i^ddb+ic^ddb\}
541  ineighbor=neighbor_child(1,inc^dd,igrid)
542  ipe_neighbor=neighbor_child(2,inc^dd,igrid)
543  ixmin^d=ix^d%ixmin^dd=ixmlo^dd+(ic^dd-1)*nxco^dd;
544  ixmax^d=ix^d%ixmax^dd=ixmin^dd-1+nxco^dd;
545  if (ipe_neighbor==mype) then
546  iotherside=3-iside
547  if (slab_uniform) then
548  psb(igrid)%w(ix^s,nw0:nw1) &
549  = psb(igrid)%w(ix^s,nw0:nw1) &
550  + pflux(iotherside,^d,ineighbor)%flux(:^dd&,1:nwfluxin)&
551  * cofiratio
552  else
553  do iw=nw0,nw1
554  psb(igrid)%w(ix^s,iw)=psb(igrid)%w(ix^s,iw) &
555  +pflux(iotherside,^d,ineighbor)%flux(:^dd&,iw-nw0+1) &
556  /ps(igrid)%dvolume(ix^s)
557  end do
558  end if
559  else
560  if (slab_uniform) then
561  ibufnext=ibuf+isize(^d)
562  if(stagger_grid) ibufnext=ibufnext-isize_stg(^d)
563  psb(igrid)%w(ix^s,nw0:nw1) &
564  = psb(igrid)%w(ix^s,nw0:nw1)+cofiratio &
565  *reshape(source=recvbuffer(ibuf:ibufnext-1), &
566  shape=shape(psb(igrid)%w(ix^s,nw0:nw1)))
567  ibuf=ibuf+isize(^d)
568  else
569  ibufnext=ibuf+isize(^d)
570  if(stagger_grid) then
571  nbuf=(isize(^d)-isize_stg(^d))/nwfluxin
572  else
573  nbuf=isize(^d)/nwfluxin
574  end if
575  do iw=nw0,nw1
576  psb(igrid)%w(ix^s,iw)=psb(igrid)%w(ix^s,iw) &
577  +reshape(source=recvbuffer(ibuf:ibufnext-1), &
578  shape=shape(psb(igrid)%w(ix^s,iw))) &
579  /ps(igrid)%dvolume(ix^s)
580  ibuf=ibuf+nbuf
581  end do
582  ibuf=ibufnext
583  end if
584  end if
585  {end do\}
586  end do\}
587  end select
588  end do
589  end do
590 
591  if (nsend>0) then
592  call mpi_waitall(nsend,fc_sendreq,fc_sendstat,ierrmpi)
593  end if
594 
595  end subroutine fix_conserve
596 
597  subroutine store_flux(igrid,fC,idim^LIM,nwfluxin)
599 
600  integer, intent(in) :: igrid, idim^lim, nwfluxin
601  double precision, intent(in) :: fc(ixg^t,1:nwfluxin,1:ndim)
602 
603  integer :: idims, iside, i^d, ic^d, inc^d, ix^d, ixco^d, nxco^d, iw
604 
605  do idims = idim^lim
606  select case (idims)
607  {case (^d)
608  do iside=1,2
609  i^dd=kr(^dd,^d)*(2*iside-3);
610 
611  if (neighbor_pole(i^dd,igrid)/=0) cycle
612 
613  select case (neighbor_type(i^dd,igrid))
614  case (neighbor_fine)
615  select case (iside)
616  case (1)
617  pflux(iside,^d,igrid)%flux(1^d%:^dd&,1:nwfluxin) = &
618  -fc(nghostcells^d%ixM^t,1:nwfluxin,^d)
619  case (2)
620  pflux(iside,^d,igrid)%flux(1^d%:^dd&,1:nwfluxin) = &
621  fc(ixmhi^d^d%ixM^t,1:nwfluxin,^d)
622  end select
623  case (neighbor_coarse)
624  nxco^d=1^d%nxCo^dd=ixghi^dd/2-nghostcells;
625  select case (iside)
626  case (1)
627  do iw=1,nwfluxin
628  {do ixco^ddb=1,nxco^ddb\}
629  ix^d=nghostcells^d%ix^dd=ixmlo^dd+2*(ixco^dd-1);
630  pflux(iside,^d,igrid)%flux(ixco^dd,iw) &
631  = {^noonedsum}(fc(ix^d^d%ix^dd:ix^dd+1,iw,^d))
632  {end do\}
633  end do
634  case (2)
635  do iw=1,nwfluxin
636  {do ixco^ddb=1,nxco^ddb\}
637  ix^d=ixmhi^d^d%ix^dd=ixmlo^dd+2*(ixco^dd-1);
638  pflux(iside,^d,igrid)%flux(ixco^dd,iw) &
639  =-{^noonedsum}(fc(ix^d^d%ix^dd:ix^dd+1,iw,^d))
640  {end do\}
641  end do
642  end select
643  end select
644  end do\}
645  end select
646  end do
647 
648  end subroutine store_flux
649 
650  subroutine store_edge(igrid,ixI^L,fE,idim^LIM)
652 
653  integer, intent(in) :: igrid, ixi^l, idim^lim
654  double precision, intent(in) :: fe(ixi^s,sdim:3)
655 
656  integer :: idims, idir, iside, i^d
657  integer :: pi^d, mi^d, ph^d, mh^d ! To detect corners
658  integer :: ixmc^l
659 
660  do idims = idim^lim !loop over face directions
661  !! Loop over block faces
662  do iside=1,2
663  i^d=kr(^d,idims)*(2*iside-3);
664  if (neighbor_pole(i^d,igrid)/=0) cycle
665  select case (neighbor_type(i^d,igrid))
666  case (neighbor_fine)
667  ! The neighbour is finer
668  ! Face direction, side (left or right), restrict ==ired?, fE
669  call flux_to_edge(igrid,ixi^l,idims,iside,.false.,fe)
670  case(neighbor_coarse)
671  ! The neighbour is coarser
672  call flux_to_edge(igrid,ixi^l,idims,iside,.true.,fe)
673  case(neighbor_sibling)
674  ! If the neighbour is at the same level,
675  ! check if there are corners
676  ! If there is any corner, store the fluxes from that side
677  do idir=idims+1,ndim
678  pi^d=i^d+kr(idir,^d);
679  mi^d=i^d-kr(idir,^d);
680  ph^d=pi^d-kr(idims,^d)*(2*iside-3);
681  mh^d=mi^d-kr(idims,^d)*(2*iside-3);
682  if (neighbor_type(pi^d,igrid)==4&
683  .and.neighbor_type(ph^d,igrid)==3) then
684  call flux_to_edge(igrid,ixi^l,idims,iside,.false.,fe)
685  end if
686  if (neighbor_type(mi^d,igrid)==4&
687  .and.neighbor_type(mh^d,igrid)==3) then
688  call flux_to_edge(igrid,ixi^l,idims,iside,.false.,fe)
689  end if
690  end do
691  end select
692  end do
693  end do
694 
695  end subroutine store_edge
696 
697  subroutine flux_to_edge(igrid,ixI^L,idims,iside,restrict,fE)
699 
700  integer :: igrid,ixI^L,idims,iside
701  logical :: restrict
702  double precision, intent(in) :: fE(ixI^S,sdim:3)
703 
704  integer :: idir1,idir2
705  integer :: ixE^L,ixF^L{^IFTHREED, jxF^L,}, nx^D,nxCo^D
706 
707  nx^d=ixmhi^d-ixmlo^d+1;
708  nxco^d=nx^d/2;
709  ! ixE are the indices on the 'edge' array.
710  ! ixF are the indices on the 'fE' array
711  ! jxF are indices advanced to perform the flux restriction (sum) in 3D
712  ! A line integral of the electric field on the coarse side
713  ! lies over two edges on the fine side. So, in 3D we restrict by summing
714  ! over two cells on the fine side.
715 
716  do idir1=1,ndim-1
717  {^ifthreed ! 3D: rotate indices among 1 and 2 to save space
718  idir2=mod(idir1+idims-1,3)+1}
719  {^iftwod ! Assign the only field component (3) to the only index (1)
720  idir2=3}
721 
722  if (restrict) then
723  ! Set up indices for restriction
724  ixfmin^d=ixmlo^d-1+kr(^d,idir2);
725  ixfmax^d=ixmhi^d-kr(^d,idir2);
726  {^ifthreed
727  jxf^l=ixf^l+kr(^d,idir2);}
728 
729  ixemin^d=0+kr(^d,idir2);
730  ixemax^d=nxco^d;
731  select case(idims)
732  {case(^d)
733  ixemin^d=1;ixemax^d=1;
734  select case(iside)
735  case(1)
736  ixfmax^d=ixfmin^d
737  {^ifthreedjxfmax^d=ixfmin^d}
738  case(2)
739  ixfmin^d=ixfmax^d
740  {^ifthreedjxfmin^d=ixfmax^d}
741  end select
742  \}
743  end select
744 
745  pflux(iside,idims,igrid)%edge(ixe^s,idir1)=&
746  fe(ixfmin^d:ixfmax^d:2,idir2){^ifthreed +&
747  fe(jxfmin^d:jxfmax^d:2,idir2)};
748 
749  else
750  ! Set up indices for copying
751  ixfmin^d=ixmlo^d-1+kr(^d,idir2);
752  ixfmax^d=ixmhi^d;
753  ixemin^d=0+kr(^d,idir2);
754  ixemax^d=nx^d;
755 
756  select case(idims)
757  {case(^d)
758  ixemin^d=1;ixemax^d=1;
759  select case(iside)
760  case(1)
761  ixfmax^d=ixfmin^d
762  case(2)
763  ixfmin^d=ixfmax^d
764  end select
765  \}
766  end select
767 
768  pflux(iside,idims,igrid)%edge(ixe^s,idir1)=fe(ixf^s,idir2)
769 
770  end if
771 
772  end do
773 
774  end subroutine flux_to_edge
775 
776  subroutine fix_edges(psuse,idim^LIM)
778 
779  type(state) :: psuse(max_blocks)
780  integer, intent(in) :: idim^lim
781 
782  integer :: iigrid, igrid, idims, iside, iotherside, i^d, ic^d, inc^d, ixmc^l
783  integer :: nbuf, ibufnext
784  integer :: ibufnext_cc
785  integer :: pi^d, mi^d, ph^d, mh^d ! To detect corners
786  integer :: ixe^l(1:3), ixte^l, ixf^l(1:ndim), ixfe^l(1:3)
787  integer :: nx^d, idir, ix, ipe_neighbor, ineighbor
788  logical :: pcorner(1:ndim),mcorner(1:ndim)
789 
790  if (nrecv_ct>0) then
791  call mpi_waitall(nrecv_ct,cc_recvreq,cc_recvstat,ierrmpi)
792  end if
793 
794  ! Initialise buffer counter again
795  ibuf=1
796  ibuf_cc=1
797  do iigrid=1,igridstail; igrid=igrids(iigrid);
798  do idims= idim^lim
799  do iside=1,2
800  i^d=kr(^d,idims)*(2*iside-3);
801  if (neighbor_pole(i^d,igrid)/=0) cycle
802  select case(neighbor_type(i^d,igrid))
803  case(neighbor_fine)
804  ! The first neighbour is finer
805  if (.not.neighbor_active(i^d,igrid).or.&
806  .not.neighbor_active(0^d&,igrid) ) then
807  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
808  inc^db=2*i^db+ic^db\}
809  ipe_neighbor=neighbor_child(2,inc^d,igrid)
810  !! When the neighbour is in a different process
811  if (ipe_neighbor/=mype) then
812  ibufnext=ibuf+isize(idims)
813  ibuf=ibufnext
814  end if
815  {end do\}
816  cycle
817  end if
818 
819  ! Check if there are corners
820  pcorner=.false.
821  mcorner=.false.
822  do idir=1,ndim
823  pi^d=i^d+kr(idir,^d);
824  mi^d=i^d-kr(idir,^d);
825  ph^d=pi^d-kr(idims,^d)*(2*iside-3);
826  mh^d=mi^d-kr(idims,^d)*(2*iside-3);
827  if (neighbor_type(ph^d,igrid)==neighbor_fine) pcorner(idir)=.true.
828  if (neighbor_type(mh^d,igrid)==neighbor_fine) mcorner(idir)=.true.
829  end do
830  ! Calculate indices range
831  call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.false.,.false.,0^d&,pcorner,mcorner)
832  ! Remove coarse part of circulation
833  call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iside,idims,igrid)%edge,idims,iside,.false.,psuse(igrid))
834  ! Add fine part of the circulation
835  {do ic^db=1+int((1-i^db)/2),2-int((1+i^db)/2)
836  inc^db=2*i^db+ic^db\}
837  ineighbor=neighbor_child(1,inc^d,igrid)
838  ipe_neighbor=neighbor_child(2,inc^d,igrid)
839  iotherside=3-iside
840  nx^d=(ixmhi^d-ixmlo^d+1)/2;
841  call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.true.,.false.,inc^d,pcorner,mcorner)
842  if (ipe_neighbor==mype) then
843  call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iotherside,idims,ineighbor)%edge,idims,iside,.true.,psuse(igrid))
844  else
845  ibufnext=ibuf+isize(idims)
846  call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,&
847  reshape(source=recvbuffer(ibufnext-isize_stg(idims):ibufnext-1),&
848  shape=(/ ixtemax^d-ixtemin^d+1 ,^nd-1 /)),&
849  idims,iside,.true.,psuse(igrid))
850  ibuf=ibufnext
851  end if
852  {end do\}
853 
854  case(neighbor_sibling)
855  ! The first neighbour is at the same level
856  ! Check if there are corners
857  do idir=idims+1,ndim
858  pcorner=.false.
859  mcorner=.false.
860  pi^d=i^d+kr(idir,^d);
861  mi^d=i^d-kr(idir,^d);
862  ph^d=pi^d-kr(idims,^d)*(2*iside-3);
863  mh^d=mi^d-kr(idims,^d)*(2*iside-3);
864  if (neighbor_type(pi^d,igrid)==neighbor_fine&
865  .and.neighbor_type(ph^d,igrid)==neighbor_sibling&
866  .and.neighbor_pole(pi^d,igrid)==0) then
867  pcorner(idir)=.true.
868  ! Remove coarse part
869  ! Set indices
870  call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.false.,.true.,0^d&,pcorner,mcorner)
871  ! Remove
872  call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iside,idims,igrid)%edge,idims,iside,.false.,psuse(igrid))
873  ! Add fine part
874  ! Find relative position of finer grid
875  {^ifthreed do ix=1,2}
876  inc^d=kr(idims,^d)*3*(iside-1)+3*kr(idir,^d){^ifthreed+kr(6-idir-idims,^d)*ix};
877  ineighbor=neighbor_child(1,inc^d,igrid)
878  ipe_neighbor=neighbor_child(2,inc^d,igrid)
879  iotherside=3-iside
880  ! Set indices
881  call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.true.,.true.,inc^d,pcorner,mcorner)
882  ! add
883  if (ipe_neighbor==mype) then
884  call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iotherside,idims,ineighbor)%edge,idims,iside,.true.,psuse(igrid))
885  else
886  ibufnext_cc=ibuf_cc+isize_stg(idims)
887  call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,&
888  reshape(source=recvbuffer_cc(ibuf_cc:ibufnext_cc-1),&
889  shape=(/ ixtemax^d-ixtemin^d+1 ,^nd-1 /)),&
890  idims,iside,.true.,psuse(igrid))
891  ibuf_cc=ibufnext_cc
892  end if
893  {^ifthreed end do}
894  ! Set CoCorner to false again for next step
895  pcorner(idir)=.false.
896  end if
897 
898  if (neighbor_type(mi^d,igrid)==neighbor_fine&
899  .and.neighbor_type(mh^d,igrid)==neighbor_sibling&
900  .and.neighbor_pole(mi^d,igrid)==0) then
901  mcorner(idir)=.true.
902  ! Remove coarse part
903  ! Set indices
904  call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.false.,.true.,0^d&,pcorner,mcorner)
905  ! Remove
906  call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iside,idims,igrid)%edge,idims,iside,.false.,psuse(igrid))
907  ! Add fine part
908  ! Find relative position of finer grid
909  {^ifthreed do ix=1,2}
910  inc^d=kr(idims,^d)*3*(iside-1){^ifthreed+kr(6-idir-idims,^d)*ix};
911  ineighbor=neighbor_child(1,inc^d,igrid)
912  ipe_neighbor=neighbor_child(2,inc^d,igrid)
913  iotherside=3-iside
914  ! Set indices
915  call set_ix_circ(ixf^l,ixte^l,ixe^l,ixfe^l,igrid,idims,iside,.true.,.true.,inc^d,pcorner,mcorner)
916  ! add
917  if (ipe_neighbor==mype) then
918  call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,pflux(iotherside,idims,ineighbor)%edge,idims,iside,.true.,psuse(igrid))
919  else
920  ibufnext_cc=ibuf_cc+isize_stg(idims)
921  call add_sub_circ(ixf^l,ixte^l,ixe^l,ixfe^l,&
922  reshape(source=recvbuffer_cc(ibuf_cc:ibufnext_cc-1),&
923  shape=(/ ixtemax^d-ixtemin^d+1 ,^nd-1 /)),&
924  idims,iside,.true.,psuse(igrid))
925  ibuf_cc=ibufnext_cc
926  end if
927  {^ifthreed end do}
928  ! Set CoCorner to false again for next step
929  mcorner(idir)=.false.
930  end if
931  end do
932  end select
933  end do
934  end do
935  end do
936 
937  if (nsend_ct>0) call mpi_waitall(nsend_ct,cc_sendreq,cc_sendstat,ierrmpi)
938 
939  end subroutine fix_edges
940 
941  !> This routine sets the indexes for the correction
942  !> of the circulation according to several different
943  !> cases, as grids located in different cpus,
944  !> presence of corners, and different relative locations
945  !> of the fine grid respect to the coarse one
946  subroutine set_ix_circ(ixF^L,ixtE^L,ixE^L,ixfE^L,igrid,idims,iside,add,CoCorner,inc^D,pcorner,mcorner)
948 
949  integer,intent(in) :: igrid,idims,iside,inc^D
950  logical,intent(in) :: add,CoCorner
951  logical,intent(inout) :: pcorner(1:ndim),mcorner(1:ndim)
952  integer,intent(out) :: ixF^L(1:ndim),ixtE^L,ixE^L(1:3),ixfE^L(1:3) ! Indices for faces and edges
953  integer :: icor^D,idim1,idir,nx^D,middle^D
954  integer :: ixtfE^L
955 
956  ! ixF -> Indices for the _F_aces, and
957  ! depends on the field component
958  ! ixtE -> are the _t_otal range of the 'edge' array
959  ! ixE -> are the ranges of the edge array,
960  ! depending on the component
961  ! ixfE -> are the ranges of the fE array (3D),
962  ! and also depend on the component
963 
964  ! ... General ...
965  ! Assign indices for the size of the E field array
966 
967  ixtfemin^d=ixmlo^d-1;
968  ixtfemax^d=ixmhi^d;
969 
970  if(add) then
971  nx^d=(ixmhi^d-ixmlo^d+1)/2;
972  else
973  nx^d=ixmhi^d-ixmlo^d+1;
974  end if
975 
976  do idim1=1,ndim
977  ixtemin^d=0;
978  ixtemax^d=nx^d;
979  select case(idims)
980  {case(^d)
981  ixtemin^d=1;ixtemax^d=1;
982  if (iside==1) ixtfemax^d=ixtfemin^d;
983  if (iside==2) ixtfemin^d=ixtfemax^d;
984  \}
985  end select
986  end do
987 
988  ! Assign indices, considering only the face
989  ! (idims and iside)
990  do idim1=1,ndim
991  ixfmin^d(idim1)=ixmlo^d-kr(idim1,^d);
992  ixfmax^d(idim1)=ixmhi^d;
993  select case(idims)
994  {case(^d)
995  select case(iside)
996  case(1)
997  ixfmax^d(idim1)=ixfmin^d(idim1)
998  case(2)
999  ixfmin^d(idim1)=ixfmax^d(idim1)
1000  end select
1001  \}
1002  end select
1003  end do
1004  ! ... Relative position ...
1005  ! Restrict range using relative position
1006  if(add) then
1007  middle^d=(ixmhi^d+ixmlo^d)/2;
1008  {
1009  if(inc^d==1) then
1010  ixfmax^d(:)=middle^d
1011  ixtfemax^d=middle^d
1012  end if
1013  if(inc^d==2) then
1014  ixfmin^d(:)=middle^d+1
1015  ixtfemin^d=middle^d
1016  end if
1017  \}
1018  end if
1019  ! ... Adjust ranges of edges according to direction ...
1020  do idim1=1,3
1021  ixfemax^d(idim1)=ixtfemax^d;
1022  ixemax^d(idim1)=ixtemax^d;
1023  ixfemin^d(idim1)=ixtfemin^d+kr(idim1,^d);
1024  ixemin^d(idim1)=ixtemin^d+kr(idim1,^d);
1025  end do
1026  ! ... Corners ...
1027  ! 'Coarse' corners
1028  if (cocorner) then
1029  do idim1=idims+1,ndim
1030  if (pcorner(idim1)) then
1031  do idir=1,3!Index arrays have size ndim
1032  if (idir==6-idim1-idims) then
1033  !!! Something here has to change
1034  !!! Array ixfE must have size 3, while
1035  !!! ixE must have size ndim
1036  {if (^d==idim1) then
1037  ixfemin^d(idir)=ixfemax^d(idir)
1038  if (add) then
1039  ixemax^d(idir) =ixemin^d(idir)
1040  else
1041  ixemin^d(idir) =ixemax^d(idir)
1042  end if
1043  end if\}
1044  else
1045  ixemin^d(idir)=1;
1046  ixemax^d(idir)=0;
1047  ixfemin^d(idir)=1;
1048  ixfemax^d(idir)=0;
1049  end if
1050  end do
1051  end if
1052  if (mcorner(idim1)) then
1053  do idir=1,3
1054  if (idir==6-idim1-idims) then
1055  {if (^d==idim1) then
1056  ixfemax^d(idir)=ixfemin^d(idir)
1057  if (add) then
1058  ixemin^d(idir) =ixemax^d(idir)
1059  else
1060  ixemax^d(idir) =ixemin^d(idir)
1061  end if
1062  end if\}
1063  else
1064  ixemin^d(idir)=1;
1065  ixemax^d(idir)=0;
1066  ixfemin^d(idir)=1;
1067  ixfemax^d(idir)=0;
1068  end if
1069  end do
1070  end if
1071  end do
1072  else
1073  ! Other kinds of corners
1074  ! Crop ranges to account for corners
1075  ! When the fine fluxes are added, we consider
1076  ! whether they come from the same cpu or from
1077  ! a different one, in order to minimise the
1078  ! amount of communication
1079  ! Case for different processors still not implemented!!!
1080  {if((idims.gt.^d).and.pcorner(^d)) then
1081  if((.not.add).or.(inc^d==2)) then
1082  !ixFmax^DD(:)=ixFmax^DD(:)-kr(^D,^DD);
1083  do idir=1,3
1084  if ((idir==idims).or.(idir==^d)) cycle
1085  ixfemax^d(idir)=ixfemax^d(idir)-1
1086  ixemax^d(idir)=ixemax^d(idir)-1
1087  end do
1088  end if
1089  end if\}
1090  {if((idims>^d).and.mcorner(^d)) then
1091  if((.not.add).or.(inc^d==1)) then
1092  !ixFmin^DD(:)=ixFmin^DD(:)+kr(^D,^DD);
1093  do idir=1,3
1094  if ((idir==idims).or.(idir==^d)) cycle
1095  ixfemin^d(idir)=ixfemin^d(idir)+1
1096  ixemin^d(idir)=ixemin^d(idir)+1
1097  end do
1098  end if
1099  end if\}
1100  end if
1101 
1102  end subroutine set_ix_circ
1103 
1104  subroutine add_sub_circ(ixF^L,ixtE^L,ixE^L,ixfE^L,edge,idims,iside,add,s)
1106 
1107  type(state) :: s
1108  integer,intent(in) :: idims,iside
1109  integer :: ixF^L(1:ndim),ixtE^L,ixE^L(1:3),ixfE^L(1:3)
1110  double precision :: edge(ixtE^S,1:ndim-1)
1111  logical,intent(in) :: add
1112 
1113  integer :: idim1,idim2,idir,middle^D
1114  integer :: ixfEC^L,ixEC^L
1115  double precision :: fE(ixG^T,sdim:3)
1116  double precision :: circ(ixG^T,1:ndim)
1117  integer :: ix^L,hx^L,ixC^L,hxC^L ! Indices for edges
1118 
1119  ! ixF -> Indices for the faces, depends on the field component
1120  ! ixE -> Total range for the edges
1121  ! ixfE -> Edges in fE (3D) array
1122  ! ix,hx,ixC,hxC -> Auxiliary indices
1123  ! Assign quantities stored ad edges to make it as similar as
1124  ! possible to the routine updatefaces.
1125  fe(:^d&,:)=zero
1126  do idim1=1,ndim-1
1127  {^ifthreed ! 3D: rotate indices (see routine flux_to_edge)
1128  idir=mod(idim1+idims-1,3)+1}
1129  {^iftwod ! 2D: move E back to directon 3
1130  idir=3}
1131  ixfec^l=ixfe^l(idir);
1132  ixec^l=ixe^l(idir);
1133  fe(ixfec^s,idir)=edge(ixec^s,idim1)
1134  end do
1135 
1136  ! Calculate part of circulation needed
1137  circ=zero
1138  do idim1=1,ndim
1139  do idim2=1,ndim
1140  do idir=sdim,3
1141  if (lvc(idim1,idim2,idir)==0) cycle
1142  ! Assemble indices
1143  ixc^l=ixf^l(idim1);
1144  hxc^l=ixc^l-kr(idim2,^d);
1145  if(idim1==idims) then
1146  circ(ixc^s,idim1)=circ(ixc^s,idim1)+lvc(idim1,idim2,idir)&
1147  *(fe(ixc^s,idir)-fe(hxc^s,idir))
1148  else
1149  select case(iside)
1150  case(2)
1151  circ(ixc^s,idim1)=circ(ixc^s,idim1)+lvc(idim1,idim2,idir)*fe(ixc^s,idir)
1152  case(1)
1153  circ(ixc^s,idim1)=circ(ixc^s,idim1)-lvc(idim1,idim2,idir)*fe(hxc^s,idir)
1154  end select
1155  end if
1156  end do
1157  end do
1158  end do
1159 
1160  ! Divide circulation by surface and add
1161  do idim1=1,ndim
1162  ixc^l=ixf^l(idim1);
1163  where(s%surfaceC(ixc^s,idim1)>1.0d-9*s%dvolume(ixc^s))
1164  circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
1165  elsewhere
1166  circ(ixc^s,idim1)=zero
1167  end where
1168  ! Add/subtract to field at face
1169  if (add) then
1170  s%ws(ixc^s,idim1)=s%ws(ixc^s,idim1)-circ(ixc^s,idim1)
1171  else
1172  s%ws(ixc^s,idim1)=s%ws(ixc^s,idim1)+circ(ixc^s,idim1)
1173  end if
1174  end do
1175 
1176  end subroutine add_sub_circ
1177 
1178 end module mod_fix_conserve
Module for flux conservation near refinement boundaries.
subroutine, public recvflux(idimLIM)
subroutine, public fix_edges(psuse, idimLIM)
subroutine, public sendflux(idimLIM)
subroutine set_ix_circ(ixFL, ixtEL, ixEL, ixfEL, igrid, idims, iside, add, CoCorner, incD, pcorner, mcorner)
This routine sets the indexes for the correction of the circulation according to several different ca...
subroutine, public fix_conserve(psb, idimLIM, nw0, nwfluxin)
subroutine, public deallocatebflux
subroutine add_sub_circ(ixFL, ixtEL, ixEL, ixfEL, edge, idims, iside, add, s)
type(fluxalloc), dimension(:,:,:), allocatable, public pflux
store flux to fix conservation
subroutine, public store_edge(igrid, ixIL, fE, idimLIM)
subroutine flux_to_edge(igrid, ixIL, idims, iside, restrict, fE)
subroutine, public init_comm_fix_conserve(idimLIM, nwfluxin)
subroutine, public allocatebflux
subroutine, public store_flux(igrid, fC, idimLIM, nwfluxin)
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer ixghi
Upper index of grid block arrays.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical stagger_grid
True for using stagger grid.
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
integer, dimension(:), allocatable, parameter d
integer ierrmpi
A global MPI error return code.
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer max_blocks
The maximum number of grid blocks in a processor.
integer, dimension(:,:), allocatable node