12 double precision,
public ::
mf_nu = 1.d-15
15 double precision,
public ::
mf_vmax = 3.d6
25 double precision,
public ::
mf_eta = 0.0d0
31 double precision :: divbdiff = 0.8d0
40 integer,
allocatable,
public,
protected ::
mom(:)
43 integer,
public,
protected ::
psi_
49 integer,
parameter :: divb_none = 0
50 integer,
parameter :: divb_multigrid = -1
51 integer,
parameter :: divb_glm = 1
52 integer,
parameter :: divb_powel = 2
53 integer,
parameter :: divb_janhunen = 3
54 integer,
parameter :: divb_linde = 4
55 integer,
parameter :: divb_lindejanhunen = 5
56 integer,
parameter :: divb_lindepowel = 6
57 integer,
parameter :: divb_lindeglm = 7
58 integer,
parameter :: divb_ct = 8
64 logical,
public,
protected ::
mf_glm = .false.
82 logical :: compactres = .false.
91 character(len=std_len),
public,
protected ::
typedivbfix =
'ct'
94 character(len=std_len),
public,
protected ::
type_ct =
'average'
97 character(len=std_len) :: typedivbdiff =
'all'
117 subroutine mf_read_params(files)
120 character(len=*),
intent(in) :: files(:)
131 do n = 1,
size(files)
132 open(
unitpar, file=trim(files(n)), status=
"old")
133 read(
unitpar, mf_list,
end=111)
137 end subroutine mf_read_params
142 integer,
intent(in) :: fh
143 integer,
parameter :: n_par = 1
144 double precision :: values(n_par)
145 character(len=name_len) :: names(n_par)
146 integer,
dimension(MPI_STATUS_SIZE) :: st
149 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
153 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
154 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
176 type_divb = divb_none
179 type_divb = divb_multigrid
181 mg%operator_type = mg_laplacian
188 case (
'powel',
'powell')
189 type_divb = divb_powel
191 type_divb = divb_janhunen
193 type_divb = divb_linde
194 case (
'lindejanhunen')
195 type_divb = divb_lindejanhunen
197 type_divb = divb_lindepowel
201 type_divb = divb_lindeglm
206 call mpistop(
'Unknown divB fix')
214 mag(:) = var_set_bfield(
ndir)
218 allocate(start_indices(number_species),stop_indices(number_species))
220 start_indices(1)=mag(1)
223 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
230 mom(:) = var_set_momentum(
ndir)
242 stop_indices(1)=nwflux
248 allocate(iw_vector(nvector))
249 iw_vector(1) =
mom(1) - 1
250 iw_vector(2) = mag(1) - 1
257 call mpistop(
"phys_check error: flux_type has wrong shape")
276 if(type_divb==divb_glm)
then
291 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
321 double precision :: mp,kB,miu0,c_lightspeed
381 integer,
intent(in) :: ixi^
l, ixo^
l
382 double precision,
intent(inout) :: w(ixi^s, nw)
383 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
391 integer,
intent(in) :: ixi^
l, ixo^
l
392 double precision,
intent(inout) :: w(ixi^s, nw)
393 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
402 integer,
intent(in) :: ixI^L, ixO^L, idim
403 double precision,
intent(in) :: w(ixI^S, nw), x(ixI^S,1:ndim)
404 double precision,
intent(inout) :: cmax(ixI^S)
406 cmax(ixo^s)=abs(w(ixo^s,
mom(idim)))+one
411 subroutine mf_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
415 integer,
intent(in) :: ixI^L, ixO^L, idim
416 double precision,
intent(in) :: wLC(ixI^S, nw), wRC(ixI^S, nw)
417 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
418 double precision,
intent(in) :: x(ixI^S,1:ndim)
419 double precision,
intent(inout) :: cmax(ixI^S,1:number_species)
420 double precision,
intent(inout),
optional :: cmin(ixI^S,1:number_species)
421 double precision,
intent(in) :: Hspeed(ixI^S,1:number_species)
423 double precision,
dimension(ixI^S) :: tmp1
425 tmp1(ixo^s)=0.5d0*(wlc(ixo^s,
mom(idim))+wrc(ixo^s,
mom(idim)))
426 if(
present(cmin))
then
427 cmax(ixo^s,1)=max(tmp1(ixo^s)+one,zero)
428 cmin(ixo^s,1)=min(tmp1(ixo^s)-one,zero)
430 cmax(ixo^s,1)=abs(tmp1(ixo^s))+one
439 integer,
intent(in) :: ixI^L, ixO^L, idim
440 double precision,
intent(in) :: wLp(ixI^S, nw), wRp(ixI^S, nw)
441 double precision,
intent(in) :: cmax(ixI^S)
442 double precision,
intent(in),
optional :: cmin(ixI^S)
443 type(ct_velocity),
intent(inout):: vcts
445 integer :: idimE,idimN
451 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
453 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
455 if(.not.
allocated(vcts%vbarC))
then
456 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
457 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
460 if(
present(cmin))
then
461 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
462 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
464 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
465 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
468 idimn=mod(idim,
ndir)+1
469 idime=mod(idim+1,
ndir)+1
471 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
472 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
473 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
474 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
475 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
477 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
478 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
479 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
480 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
481 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
483 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
494 integer,
intent(in) :: ixI^L, ixO^L, idim
496 double precision,
intent(in) :: wC(ixI^S,nw)
498 double precision,
intent(in) :: w(ixI^S,nw)
499 double precision,
intent(in) :: x(ixI^S,1:ndim)
500 double precision,
intent(out) :: f(ixI^S,nwflux)
502 integer :: idir, ix^D
505 {
do ix^db=ixomin^db,ixomax^db\}
508 f(ix^d,mag(idir))=w(ix^d,
mom(idim))*w(ix^d,mag(idir))-w(ix^d,mag(idim))*w(ix^d,
mom(idir))
512 {
do ix^db=ixomin^db,ixomax^db\}
513 f(ix^d,mag(idim))=w(ix^d,
psi_)
515 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
525 double precision,
intent(in) :: qt
526 type(state),
target :: psa(max_blocks)
528 integer :: iigrid,igrid
529 logical :: stagger_flag
530 logical :: firstmf=.true.
545 do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
574 subroutine mf_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
577 integer,
intent(in) :: ixI^L, ixO^L
578 double precision,
intent(in) :: qdt, dtfactor
579 double precision,
intent(in) :: wCT(ixI^S,1:nw),wCTprim(ixI^S,1:nw), x(ixI^S,1:ndim)
580 double precision,
intent(inout) :: w(ixI^S,1:nw)
581 logical,
intent(in) :: qsourcesplit
582 logical,
intent(inout) :: active
584 if (.not. qsourcesplit)
then
587 if (abs(
mf_eta)>smalldouble)
then
602 select case (type_divb)
617 case (divb_lindejanhunen)
621 case (divb_lindepowel)
631 case (divb_multigrid)
634 call mpistop(
'Unknown divB fix')
644 integer,
intent(in) :: ixI^L, ixO^L
645 double precision,
intent(in) :: x(ixI^S,1:ndim)
646 double precision,
intent(inout) :: w(ixI^S,1:nw)
648 double precision :: xmin(ndim),xmax(ndim)
649 double precision :: decay(ixO^S)
650 double precision :: current(ixI^S,7-2*ndir:3),tmp(ixO^S)
651 integer :: ix^D, idirmin,idir,jdir,kdir
657 if(
block%is_physical_boundary(2*^d-1))
then
660 if(2*^d-1==5.and.ndim==3)
then
662 current(ixomin^d^d%ixO^s,:)=current(ixomin^d+1^d%ixO^s,:)
667 current(ixomin^d^d%ixO^s,:)=current(ixomin^d+1^d%ixO^s,:)
671 if(
block%is_physical_boundary(2*^d))
then
672 current(ixomax^d^d%ixO^s,:)=current(ixomax^d-1^d%ixO^s,:)
677 do idir=1,ndir;
do jdir=idirmin,3;
do kdir=1,ndir
678 if(
lvc(idir,jdir,kdir)/=0)
then
679 if(
lvc(idir,jdir,kdir)==1)
then
680 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+current(ixo^s,jdir)*w(ixo^s,mag(kdir))
682 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-current(ixo^s,jdir)*w(ixo^s,mag(kdir))
685 end do;
end do;
end do
687 tmp(ixo^s)=sum(w(ixo^s,mag(:))**2,dim=ndim+1)
689 where(tmp(ixo^s)/=0.d0)
690 tmp(ixo^s)=1.d0/(tmp(ixo^s)*
mf_nu)
694 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))*tmp(ixo^s)
698 ^d&xmin(^d)=xprobmin^d\
699 ^d&xmax(^d)=xprobmax^d\
709 tmp(ixo^s)=sqrt(sum(w(ixo^s,
mom(:))**2,dim=ndim+1))/
mf_vmax+1.d-12
710 tmp(ixo^s)=dtanh(tmp(ixo^s))/tmp(ixo^s)
712 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))*tmp(ixo^s)*decay(ixo^s)
726 integer,
intent(in) :: ixI^L, ixO^L
727 double precision,
intent(in) :: qdt
728 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
729 double precision,
intent(inout) :: w(ixI^S,1:nw)
732 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
733 double precision :: gradeta(ixI^S,1:ndim), Bf(ixI^S,1:ndir)
734 double precision :: tmp(ixI^S),tmp2(ixI^S)
735 integer :: ixA^L,idir,jdir,kdir,idirmin,idim,jxO^L,hxO^L,ix
736 integer :: lxO^L, kxO^L
745 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
746 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
753 gradeta(ixo^s,1:ndim)=zero
758 call gradient(eta,ixi^l,ixo^l,idim,tmp)
759 gradeta(ixo^s,idim)=tmp(ixo^s)
763 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
769 tmp2(ixi^s)=bf(ixi^s,idir)
771 lxo^l=ixo^l+2*
kr(idim,^
d);
772 jxo^l=ixo^l+
kr(idim,^
d);
773 hxo^l=ixo^l-
kr(idim,^
d);
774 kxo^l=ixo^l-2*
kr(idim,^
d);
775 tmp(ixo^s)=tmp(ixo^s)+&
776 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
781 tmp2(ixi^s)=bf(ixi^s,idir)
783 jxo^l=ixo^l+
kr(idim,^
d);
784 hxo^l=ixo^l-
kr(idim,^
d);
785 tmp(ixo^s)=tmp(ixo^s)+&
786 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
791 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
795 do jdir=1,ndim;
do kdir=idirmin,3
796 if (
lvc(idir,jdir,kdir)/=0)
then
797 if (
lvc(idir,jdir,kdir)==1)
then
798 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
800 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
807 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
820 integer,
intent(in) :: ixI^L, ixO^L
821 double precision,
intent(in) :: qdt
822 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
823 double precision,
intent(inout) :: w(ixI^S,1:nw)
826 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S),curlj(ixI^S,1:3)
827 double precision :: tmpvec(ixI^S,1:3)
828 integer :: ixA^L,idir,idirmin,idirmin1
835 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
836 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
850 tmpvec(ixa^s,1:ndir)=zero
852 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
855 call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
857 if(ndim==2.and.ndir==3)
then
859 w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
862 w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
873 integer,
intent(in) :: ixI^L, ixO^L
874 double precision,
intent(in) :: qdt
875 double precision,
intent(in) :: wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
876 double precision,
intent(inout) :: w(ixI^S,1:nw)
878 double precision :: current(ixI^S,7-2*ndir:3)
879 double precision :: tmpvec(ixI^S,1:3),tmpvec2(ixI^S,1:3),tmp(ixI^S),ehyper(ixI^S,1:3)
880 integer :: ixA^L,idir,jdir,kdir,idirmin,idirmin1
883 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
884 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
887 tmpvec(ixa^s,1:ndir)=zero
889 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
893 call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
896 tmpvec(ixa^s,1:ndir)=zero
897 call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
898 ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*
mf_eta_hyper
901 tmpvec2(ixa^s,1:ndir)=zero
902 call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
905 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
917 integer,
intent(in) :: ixI^L, ixO^L
918 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
919 double precision,
intent(inout) :: w(ixI^S,1:nw)
920 double precision:: divb(ixI^S)
921 double precision :: gradPsi(ixI^S)
956 integer,
intent(in) :: ixI^L, ixO^L
957 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
958 double precision,
intent(inout) :: w(ixI^S,1:nw)
960 double precision :: divb(ixI^S)
968 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,
mom(idir))*divb(ixo^s)
978 integer,
intent(in) :: ixI^L, ixO^L
979 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
980 double precision,
intent(inout) :: w(ixI^S,1:nw)
981 double precision :: divb(ixI^S)
989 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,
mom(idir))*divb(ixo^s)
999 integer,
intent(in) :: ixI^L, ixO^L
1000 double precision,
intent(in) :: qdt, wCT(ixI^S,1:nw), x(ixI^S,1:ndim)
1001 double precision,
intent(inout) :: w(ixI^S,1:nw)
1002 double precision :: divb(ixI^S),graddivb(ixI^S)
1003 integer :: idim, idir, ixp^L, i^D, iside
1004 logical,
dimension(-1:1^D&) :: leveljump
1012 if(i^d==0|.and.) cycle
1013 if(neighbor_type(i^d,
block%igrid)==2 .or. neighbor_type(i^d,
block%igrid)==4)
then
1014 leveljump(i^d)=.true.
1016 leveljump(i^d)=.false.
1025 i^dd=kr(^dd,^d)*(2*iside-3);
1026 if (leveljump(i^dd))
then
1028 ixpmin^d=ixomin^d-i^d
1030 ixpmax^d=ixomax^d-i^d
1041 select case(typegrad)
1043 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
1045 call gradients(divb,ixi^l,ixp^l,idim,graddivb)
1049 if (slab_uniform)
then
1050 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
1052 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
1053 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
1056 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
1067 integer,
intent(in) :: ixi^
l, ixo^
l
1068 double precision,
intent(in) :: w(ixi^s,1:nw)
1069 double precision :: divb(ixi^s), dsurface(ixi^s)
1071 double precision :: invb(ixo^s)
1072 integer :: ixa^
l,idims
1074 call get_divb(w,ixi^
l,ixo^
l,divb)
1075 invb(ixo^s)=sqrt(sum(w(ixo^s, mag(:))**2, dim=
ndim+1))
1076 where(invb(ixo^s)/=0.d0)
1077 invb(ixo^s)=1.d0/invb(ixo^s)
1080 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
1082 ixamin^
d=ixomin^
d-1;
1083 ixamax^
d=ixomax^
d-1;
1084 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
1086 ixa^
l=ixo^
l-
kr(idims,^
d);
1087 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
1089 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
1090 block%dvolume(ixo^s)/dsurface(ixo^s)
1101 integer,
intent(in) :: ixo^
l, ixi^
l
1102 double precision,
intent(in) :: w(ixi^s,1:nw)
1103 integer,
intent(out) :: idirmin
1104 logical,
intent(in),
optional :: fourthorder
1107 double precision :: current(ixi^s,7-2*
ndir:3)
1108 integer :: idir, idirmin0
1112 if(
present(fourthorder))
then
1125 integer,
intent(in) :: ixI^L, ixO^L
1126 double precision,
intent(inout) :: dtnew
1127 double precision,
intent(in) :: dx^D
1128 double precision,
intent(in) :: w(ixI^S,1:nw)
1129 double precision,
intent(in) :: x(ixI^S,1:ndim)
1131 double precision :: dxarr(ndim)
1132 double precision :: current(ixI^S,7-2*ndir:3),eta(ixI^S)
1133 integer :: idirmin,idim
1140 else if (
mf_eta<zero)
then
1147 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
1150 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
1169 integer,
intent(in) :: ixI^L, ixO^L
1170 double precision,
intent(in) :: qdt, dtfactor, x(ixI^S,1:ndim)
1171 double precision,
intent(inout) :: wCT(ixI^S,1:nw), wprim(ixI^S,1:nw), w(ixI^S,1:nw)
1173 double precision :: tmp,invr,cs2
1174 integer :: iw,idir,ix^D
1176 integer :: mr_,mphi_
1177 integer :: br_,bphi_
1180 br_=mag(1); bphi_=mag(1)-1+
phi_
1186 {
do ix^db=ixomin^db,ixomax^db\}
1187 w(ix^d,bphi_)=w(ix^d,bphi_)+qdt/x(ix^d,1)*&
1188 (wct(ix^d,bphi_)*wct(ix^d,mr_) &
1189 -wct(ix^d,br_)*wct(ix^d,mphi_))
1193 if(
mf_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
1195 {
do ix^db=ixomin^db,ixomax^db\}
1199 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wct(ix^d,
psi_)
1204 if(.not.stagger_grid)
then
1205 tmp=wct(ix^d,
mom(1))*wct(ix^d,mag(2))-wct(ix^d,
mom(2))*wct(ix^d,mag(1))
1206 cs2=1.d0/tan(x(ix^d,2))
1208 tmp=tmp+cs2*wct(ix^d,
psi_)
1210 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
1217 if(.not.stagger_grid)
then
1218 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
1219 (wct(ix^d,
mom(1))*wct(ix^d,mag(3)) &
1220 -wct(ix^d,
mom(3))*wct(ix^d,mag(1)))
1225 if(.not.stagger_grid)
then
1226 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
1227 (wct(ix^d,
mom(1))*wct(ix^d,mag(3)) &
1228 -wct(ix^d,
mom(3))*wct(ix^d,mag(1)) &
1229 -(wct(ix^d,
mom(3))*wct(ix^d,mag(2)) &
1230 -wct(ix^d,
mom(2))*wct(ix^d,mag(3)))*cs2)
1242 integer,
intent(in) :: ixI^L, ixO^L, idir
1243 double precision,
intent(in) :: qt
1244 double precision,
intent(inout) :: wLC(ixI^S,1:nw), wRC(ixI^S,1:nw)
1245 double precision,
intent(inout) :: wLp(ixI^S,1:nw), wRp(ixI^S,1:nw)
1247 double precision :: dB(ixI^S), dPsi(ixI^S)
1255 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
1256 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
1258 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
1260 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
1263 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
1270 integer,
intent(in) :: igrid
1271 type(state),
target :: psb(max_blocks)
1273 integer :: iB, idims, iside, ixO^L, i^D
1282 i^d=
kr(^d,idims)*(2*iside-3);
1283 if (neighbor_type(i^d,igrid)/=1) cycle
1284 ib=(idims-1)*2+iside
1312 integer,
intent(in) :: ixG^L,ixO^L,iB
1313 double precision,
intent(inout) :: w(ixG^S,1:nw)
1314 double precision,
intent(in) :: x(ixG^S,1:ndim)
1316 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
1317 integer :: ix^D,ixF^L
1329 do ix1=ixfmax1,ixfmin1,-1
1330 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
1331 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
1332 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
1335 do ix1=ixfmax1,ixfmin1,-1
1336 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
1337 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
1338 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
1339 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
1340 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
1341 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
1342 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
1356 do ix1=ixfmax1,ixfmin1,-1
1357 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1358 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
1359 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
1360 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
1361 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
1362 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
1365 do ix1=ixfmax1,ixfmin1,-1
1366 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1367 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
1368 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
1369 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
1370 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
1371 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
1372 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
1373 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
1374 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
1375 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
1376 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
1377 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
1378 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
1379 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
1380 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1381 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
1382 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
1383 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
1395 do ix1=ixfmin1,ixfmax1
1396 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
1397 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
1398 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
1401 do ix1=ixfmin1,ixfmax1
1402 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
1403 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
1404 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
1405 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
1406 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
1407 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
1408 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
1422 do ix1=ixfmin1,ixfmax1
1423 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1424 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
1425 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
1426 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
1427 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
1428 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
1431 do ix1=ixfmin1,ixfmax1
1432 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1433 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
1434 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
1435 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
1436 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
1437 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
1438 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
1439 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
1440 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
1441 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
1442 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
1443 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
1444 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
1445 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
1446 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1447 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
1448 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
1449 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
1461 do ix2=ixfmax2,ixfmin2,-1
1462 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
1463 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
1464 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
1467 do ix2=ixfmax2,ixfmin2,-1
1468 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
1469 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
1470 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
1471 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
1472 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
1473 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
1474 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
1488 do ix2=ixfmax2,ixfmin2,-1
1489 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
1490 ix2+1,ixfmin3:ixfmax3,mag(2)) &
1491 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
1492 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
1493 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
1494 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
1497 do ix2=ixfmax2,ixfmin2,-1
1498 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
1499 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
1500 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
1501 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
1502 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
1503 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1504 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
1505 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
1506 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1507 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
1508 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
1509 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
1510 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
1511 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
1512 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1513 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
1514 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
1515 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
1527 do ix2=ixfmin2,ixfmax2
1528 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
1529 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
1530 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
1533 do ix2=ixfmin2,ixfmax2
1534 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
1535 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
1536 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
1537 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
1538 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
1539 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
1540 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
1554 do ix2=ixfmin2,ixfmax2
1555 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
1556 ix2-1,ixfmin3:ixfmax3,mag(2)) &
1557 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
1558 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
1559 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
1560 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
1563 do ix2=ixfmin2,ixfmax2
1564 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
1565 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
1566 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
1567 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
1568 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
1569 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1570 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
1571 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
1572 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1573 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
1574 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
1575 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
1576 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
1577 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
1578 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1579 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
1580 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
1581 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
1596 do ix3=ixfmax3,ixfmin3,-1
1597 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
1598 ixfmin2:ixfmax2,ix3+1,mag(3)) &
1599 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
1600 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
1601 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
1602 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
1605 do ix3=ixfmax3,ixfmin3,-1
1606 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
1607 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
1608 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
1609 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
1610 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
1611 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1612 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
1613 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
1614 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1615 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
1616 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
1617 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
1618 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
1619 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
1620 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
1621 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
1622 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
1623 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
1636 do ix3=ixfmin3,ixfmax3
1637 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
1638 ixfmin2:ixfmax2,ix3-1,mag(3)) &
1639 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
1640 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
1641 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
1642 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
1645 do ix3=ixfmin3,ixfmax3
1646 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
1647 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
1648 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
1649 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
1650 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
1651 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1652 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
1653 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
1654 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1655 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
1656 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
1657 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
1658 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
1659 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
1660 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
1661 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
1662 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
1663 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
1668 call mpistop(
"Special boundary is not defined for this region")
1680 double precision,
intent(in) :: qdt
1681 double precision,
intent(in) :: qt
1682 logical,
intent(inout) :: active
1684 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
1685 double precision :: res
1686 double precision,
parameter :: max_residual = 1
d-3
1687 double precision,
parameter :: residual_reduction = 1
d-10
1689 integer,
parameter :: max_its = 50
1690 double precision :: residual_it(max_its), max_divb
1691 integer :: iigrid, igrid
1692 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
1695 mg%operator_type = mg_laplacian
1703 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1704 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1707 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
1708 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1710 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1711 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1714 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
1715 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1719 write(*,*)
"mf_clean_divb_multigrid warning: unknown boundary type"
1720 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1721 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1729 do iigrid = 1, igridstail
1730 igrid = igrids(iigrid);
1733 lvl =
mg%boxes(id)%lvl
1734 nc =
mg%box_size_lvl(lvl)
1740 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
1742 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
1743 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
1748 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
1751 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
1752 if (
mype == 0) print *,
"iteration vs residual"
1755 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
1756 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
1757 if (residual_it(n) < residual_reduction * max_divb)
exit
1759 if (
mype == 0 .and. n > max_its)
then
1760 print *,
"divb_multigrid warning: not fully converged"
1761 print *,
"current amplitude of divb: ", residual_it(max_its)
1762 print *,
"multigrid smallest grid: ", &
1763 mg%domain_size_lvl(:,
mg%lowest_lvl)
1764 print *,
"note: smallest grid ideally has <= 8 cells"
1765 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
1766 print *,
"note: dx/dy/dz should be similar"
1770 call mg_fas_vcycle(
mg, max_res=res)
1771 if (res < max_residual)
exit
1773 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
1778 do iigrid = 1, igridstail
1779 igrid = igrids(iigrid);
1788 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
1792 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
1794 call gradientx(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim),.false.)
1795 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
1803 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
1804 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
1817 integer,
intent(in) :: ixI^L, ixO^L
1818 double precision,
intent(in) :: qt, qdt
1820 double precision,
intent(in) :: wprim(ixI^S,1:nw)
1821 type(state) :: sCT, s
1822 type(ct_velocity) :: vcts
1823 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
1824 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
1834 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
1844 integer,
intent(in) :: ixI^L, ixO^L
1845 double precision,
intent(in) :: qt, qdt
1846 type(state) :: sCT, s
1847 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
1848 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
1850 double precision :: circ(ixI^S,1:ndim)
1851 double precision :: E_resi(ixI^S,sdim:3)
1852 integer :: ix^D,ixC^L,ixA^L,i1kr^D,i2kr^D
1853 integer :: idim1,idim2,idir,iwdim1,iwdim2
1855 associate(bfaces=>s%ws,x=>s%x)
1866 i1kr^d=
kr(idim1,^d);
1869 i2kr^d=
kr(idim2,^d);
1872 if (
lvc(idim1,idim2,idir)==1)
then
1874 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
1876 {
do ix^db=ixcmin^db,ixcmax^db\}
1877 fe(ix^d,idir)=quarter*&
1878 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
1879 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
1881 if(
mf_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
1884 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
1892 if(
associated(usr_set_electric_field)) &
1893 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
1895 circ(ixi^s,1:ndim)=zero
1901 ixcmin^d=ixomin^d-kr(idim1,^d);
1903 ixa^l=ixc^l-kr(idim2,^d);
1906 if(lvc(idim1,idim2,idir)==1)
then
1908 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
1911 else if(lvc(idim1,idim2,idir)==-1)
then
1913 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
1920 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
1922 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
1934 integer,
intent(in) :: ixI^L, ixO^L
1935 double precision,
intent(in) :: qt, qdt
1937 double precision,
intent(in) :: wp(ixI^S,1:nw)
1938 type(state) :: sCT, s
1939 type(ct_velocity) :: vcts
1940 double precision,
intent(in) :: fC(ixI^S,1:nwflux,1:ndim)
1941 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
1943 double precision :: circ(ixI^S,1:ndim)
1945 double precision :: ECC(ixI^S,sdim:3)
1947 double precision :: EL(ixI^S),ER(ixI^S)
1949 double precision :: ELC(ixI^S),ERC(ixI^S)
1951 double precision :: jce(ixI^S,sdim:3)
1952 integer :: hxC^L,ixC^L,jxC^L,ixA^L,ixB^L
1953 integer :: idim1,idim2,idir,iwdim1,iwdim2
1955 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
1959 do idim1=1,ndim;
do idim2=1,ndim;
do idir=sdim,3
1960 if(
lvc(idim1,idim2,idir)==1)
then
1961 ecc(ixi^s,idir)=ecc(ixi^s,idir)+wp(ixi^s,mag(idim1))*wp(ixi^s,
mom(idim2))
1962 else if(
lvc(idim1,idim2,idir)==-1)
then
1963 ecc(ixi^s,idir)=ecc(ixi^s,idir)-wp(ixi^s,mag(idim1))*wp(ixi^s,
mom(idim2))
1980 if (
lvc(idim1,idim2,idir)==1)
then
1982 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
1984 jxc^l=ixc^l+
kr(idim1,^
d);
1985 hxc^l=ixc^l+
kr(idim2,^
d);
1987 fe(ixc^s,idir)=quarter*&
1988 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
1989 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
1993 ixamax^
d=ixcmax^
d+
kr(idim1,^
d);
1994 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
1995 hxc^l=ixa^l+
kr(idim2,^
d);
1996 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
1997 where(vnorm(ixc^s,idim1)>0.d0)
1998 elc(ixc^s)=el(ixc^s)
1999 else where(vnorm(ixc^s,idim1)<0.d0)
2000 elc(ixc^s)=el(jxc^s)
2002 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
2004 hxc^l=ixc^l+
kr(idim2,^
d);
2005 where(vnorm(hxc^s,idim1)>0.d0)
2006 erc(ixc^s)=er(ixc^s)
2007 else where(vnorm(hxc^s,idim1)<0.d0)
2008 erc(ixc^s)=er(jxc^s)
2010 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
2012 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
2015 jxc^l=ixc^l+
kr(idim2,^
d);
2017 ixamax^
d=ixcmax^
d+
kr(idim2,^
d);
2018 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
2019 hxc^l=ixa^l+
kr(idim1,^
d);
2020 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
2021 where(vnorm(ixc^s,idim2)>0.d0)
2022 elc(ixc^s)=el(ixc^s)
2023 else where(vnorm(ixc^s,idim2)<0.d0)
2024 elc(ixc^s)=el(jxc^s)
2026 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
2028 hxc^l=ixc^l+
kr(idim1,^
d);
2029 where(vnorm(hxc^s,idim2)>0.d0)
2030 erc(ixc^s)=er(ixc^s)
2031 else where(vnorm(hxc^s,idim2)<0.d0)
2032 erc(ixc^s)=er(jxc^s)
2034 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
2036 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
2039 if(
mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2043 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
2058 circ(ixi^s,1:ndim)=zero
2063 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
2067 if(
lvc(idim1,idim2,idir)/=0)
then
2068 hxc^l=ixc^l-
kr(idim2,^
d);
2070 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2071 +
lvc(idim1,idim2,idir)&
2078 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
2079 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2081 circ(ixc^s,idim1)=zero
2084 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2097 integer,
intent(in) :: ixI^L, ixO^L
2098 double precision,
intent(in) :: qt, qdt
2099 double precision,
intent(inout) :: fE(ixI^S,sdim:3)
2100 type(state) :: sCT, s
2101 type(ct_velocity) :: vcts
2103 double precision :: vtilL(ixI^S,2)
2104 double precision :: vtilR(ixI^S,2)
2105 double precision :: btilL(ixI^S,ndim)
2106 double precision :: btilR(ixI^S,ndim)
2107 double precision :: cp(ixI^S,2)
2108 double precision :: cm(ixI^S,2)
2109 double precision :: circ(ixI^S,1:ndim)
2111 double precision :: jce(ixI^S,sdim:3)
2112 integer :: hxC^L,ixC^L,ixCp^L,jxC^L,ixCm^L
2113 integer :: idim1,idim2,idir
2115 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
2116 cbarmax=>vcts%cbarmax)
2142 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
2146 idim2=mod(idir+1,3)+1
2148 jxc^l=ixc^l+
kr(idim1,^
d);
2149 ixcp^l=ixc^l+
kr(idim2,^
d);
2152 call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
2153 vtill(ixi^s,2),vtilr(ixi^s,2))
2155 call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
2156 vtill(ixi^s,1),vtilr(ixi^s,1))
2161 call reconstruct(ixi^l,ixc^l,idim2,bfacesct(ixi^s,idim1),&
2162 btill(ixi^s,idim1),btilr(ixi^s,idim1))
2164 call reconstruct(ixi^l,ixc^l,idim1,bfacesct(ixi^s,idim2),&
2165 btill(ixi^s,idim2),btilr(ixi^s,idim2))
2169 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
2170 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
2172 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
2173 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
2177 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
2178 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
2179 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
2180 /(cp(ixc^s,1)+cm(ixc^s,1)) &
2181 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
2182 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
2183 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
2184 /(cp(ixc^s,2)+cm(ixc^s,2))
2187 if(
mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2191 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
2205 circ(ixi^s,1:ndim)=zero
2211 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
2215 if(
lvc(idim1,idim2,idir)/=0)
then
2216 hxc^l=ixc^l-
kr(idim2,^
d);
2218 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2219 +
lvc(idim1,idim2,idir)&
2226 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
2227 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2229 circ(ixc^s,idim1)=zero
2232 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2244 integer,
intent(in) :: ixI^L, ixO^L
2245 type(state),
intent(in) :: sCT, s
2247 double precision :: jce(ixI^S,sdim:3)
2250 double precision :: jcc(ixI^S,7-2*ndir:3)
2252 double precision :: xs(ixGs^T,1:ndim)
2254 double precision :: eta(ixI^S)
2255 double precision :: gradi(ixGs^T)
2256 integer :: ix^D,ixC^L,ixA^L,ixB^L,idir,idirmin,idim1,idim2
2258 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
2264 if (
lvc(idim1,idim2,idir)==0) cycle
2265 ixcmax^d=ixomax^d+1;
2266 ixcmin^d=ixomin^d+
kr(idir,^d)-2;
2267 ixbmax^d=ixcmax^d-
kr(idir,^d)+1;
2270 xs(ixb^s,:)=x(ixb^s,:)
2271 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
2272 call gradientx(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^l,idim1,gradi,.false.)
2273 if (
lvc(idim1,idim2,idir)==1)
then
2274 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
2276 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
2283 jce(ixi^s,:)=jce(ixi^s,:)*
mf_eta
2291 ixcmin^d=ixomin^d+
kr(idir,^d)-1;
2292 jcc(ixc^s,idir)=0.d0
2294 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
2295 ixamin^d=ixcmin^d+ix^d;
2296 ixamax^d=ixcmax^d+ix^d;
2297 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
2299 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
2300 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
2311 integer,
intent(in) :: ixo^
l
2318 do ix^db=ixomin^db,ixomax^db\}
2320 s%w(ix^
d,mag(1))=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
2321 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
2322 s%w(ix^
d,mag(2))=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
2323 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
2324 s%w(ix^
d,mag(3))=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
2325 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
2328 s%w(ix^
d,mag(1))=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
2329 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
2330 s%w(ix^
d,mag(2))=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
2331 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
2374 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
2375 double precision,
intent(inout) :: ws(ixis^s,1:nws)
2376 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2378 double precision :: adummy(ixis^s,1:3)
2387 double precision :: sum_jbb_ipe, sum_j_ipe, sum_l_ipe, f_i_ipe, volume_pe
2388 double precision :: sum_jbb, sum_j, sum_l, f_i, volume, cw_sin_theta
2389 integer :: iigrid, igrid, ix^
d
2390 integer :: amode, istatus(mpi_status_size)
2391 integer,
save :: fhmf
2392 logical :: patchwi(ixg^t)
2393 logical,
save :: logmfopened=.false.
2394 character(len=800) :: filename,filehead
2395 character(len=800) :: line,datastr
2402 do iigrid=1,igridstail; igrid=igrids(iigrid);
2407 ps(igrid)%x,1,patchwi)
2409 ps(igrid)%x,2,patchwi)
2411 ps(igrid)%x,3,patchwi)
2413 ps(igrid)%x,4,patchwi)
2415 call mpi_reduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
2417 call mpi_reduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,0,&
2419 call mpi_reduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,0,&
2421 call mpi_reduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,0,&
2423 call mpi_reduce(volume_pe,volume,1,mpi_double_precision,mpi_sum,0,&
2429 cw_sin_theta = sum_jbb/sum_j
2435 if(.not.logmfopened)
then
2437 write(filename,
"(a,a)") trim(
base_filename),
"_mf_metrics.csv"
2441 open(unit=20,file=trim(filename),status=
'replace')
2442 close(20, status=
'delete')
2445 amode=ior(mpi_mode_create,mpi_mode_wronly)
2446 amode=ior(amode,mpi_mode_append)
2447 call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,
ierrmpi)
2449 filehead=
" it,time,CW_sin_theta,current,Lorenz_force,f_i"
2451 call mpi_file_write(fhmf,filehead,len_trim(filehead), &
2452 mpi_character,istatus,
ierrmpi)
2453 call mpi_file_write(fhmf,achar(10),1,mpi_character,istatus,
ierrmpi)
2457 write(datastr,
'(i6,a)')
it,
','
2458 line=trim(line)//trim(datastr)
2460 line=trim(line)//trim(datastr)
2461 write(datastr,
'(es13.6,a)') cw_sin_theta,
','
2462 line=trim(line)//trim(datastr)
2463 write(datastr,
'(es13.6,a)') sum_j,
','
2464 line=trim(line)//trim(datastr)
2465 write(datastr,
'(es13.6,a)') sum_l,
','
2466 line=trim(line)//trim(datastr)
2467 write(datastr,
'(es13.6)') f_i
2468 line=trim(line)//trim(datastr)//new_line(
'A')
2469 call mpi_file_write(fhmf,line,len_trim(line),mpi_character,istatus,
ierrmpi)
2471 call mpi_file_close(fhmf,
ierrmpi)
2481 integer,
intent(in) :: ixI^L,ixO^L
2482 double precision,
intent(in):: w(ixI^S,nw),x(ixI^S,1:ndim)
2483 double precision,
intent(inout) :: volume_pe
2484 logical,
intent(out) :: patchwi(ixI^S)
2486 double precision :: xO^L
2489 {xomin^d = xprobmin^d + 0.05d0*(xprobmax^d-xprobmin^d)\}
2490 {xomax^d = xprobmax^d - 0.05d0*(xprobmax^d-xprobmin^d)\}
2492 xomin^nd = xprobmin^nd
2497 {
do ix^db=ixomin^db,ixomax^db\}
2498 if({ x(ix^dd,^d) > xomin^d .and. x(ix^dd,^d) < xomax^d | .and. })
then
2499 patchwi(ix^d)=.true.
2500 volume_pe=volume_pe+
block%dvolume(ix^d)
2502 patchwi(ix^d)=.false.
2511 integer,
intent(in) :: ixi^
l,ixo^
l,iw
2512 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2513 double precision :: w(ixi^s,nw+
nwauxio)
2514 logical,
intent(in) :: patchwi(ixi^s)
2516 double precision,
dimension(ixI^S,7-2*ndir:3) :: current
2517 double precision,
dimension(ixI^S,1:ndir) :: bvec,qvec
2519 integer :: ix^
d,idirmin0,idirmin,idir,jdir,kdir
2525 bvec(ixo^s,:)=w(ixo^s,mag(:))
2528 qvec(ixo^s,1:
ndir)=zero
2529 do idir=1,
ndir;
do jdir=idirmin,3;
do kdir=1,
ndir
2530 if(
lvc(idir,jdir,kdir)/=0)
then
2531 if(
lvc(idir,jdir,kdir)==1)
then
2532 qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
2534 qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
2537 end do;
end do;
end do
2539 {
do ix^db=ixomin^db,ixomax^db\}
2540 if(patchwi(ix^d))
then
2541 bm2=sum(bvec(ix^d,:)**2)
2542 if(bm2/=0.d0) bm2=1.d0/bm2
2544 bm2)*block%dvolume(ix^d)
2550 {
do ix^db=ixomin^db,ixomax^db\}
2558 {
do ix^db=ixomin^db,ixomax^db\}
2563 bvec(ixo^s,:)=w(ixo^s,mag(:))
2566 qvec(ixo^s,1:ndir)=zero
2567 do idir=1,ndir;
do jdir=idirmin,3;
do kdir=1,ndir
2568 if(lvc(idir,jdir,kdir)/=0)
then
2569 if(lvc(idir,jdir,kdir)==1)
then
2570 qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
2572 qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
2575 end do;
end do;
end do
2577 {
do ix^db=ixomin^db,ixomax^db\}
2579 sqrt(sum(qvec(ix^d,:)**2))*block%dvolume(ix^d)
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
subroutine b_from_vector_potentiala(ixIsL, ixIL, ixOL, ws, x, A)
calculate magnetic field from vector potential A at cell edges
subroutine reconstruct(ixIL, ixCL, idir, q, qL, qR)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
Module with basic grid data structures.
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
subroutine, public get_divb(w, ixIL, ixOL, divb, fourthorder)
Calculate div B within ixO.
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter cylindrical
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine gradients(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir first use limiter to go from cell cente...
subroutine gradientx(q, x, ixIL, ixOL, idir, gradq, fourth_order)
Calculate gradient of a scalar q in direction idir at cell interfaces.
update ghost cells of all blocks including physical boundaries
integer, dimension(0:3^d &), target type_recv_p_p1
integer, dimension( :^d &), pointer type_recv_r
integer, dimension(-1:1^d &), target type_send_srl_p1
integer, dimension(0:3^d &), target type_send_p_f
integer, dimension( :^d &), pointer type_send_p
integer, dimension( :^d &), pointer type_send_srl
integer, dimension(-1:1^d &), target type_send_r_p1
subroutine create_bc_mpi_datatype(nwstart, nwbc)
integer, dimension(0:3^d &), target type_recv_r_f
integer, dimension(-1:1^d &), target type_recv_srl_f
integer, dimension(-1:1^d &), target type_send_r_f
integer, dimension(-1:1^d &), target type_send_srl_f
integer, dimension(0:3^d &), target type_send_p_p1
integer, dimension( :^d &), pointer type_send_r
integer, dimension(0:3^d &), target type_recv_r_p1
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
integer, dimension( :^d &), pointer type_recv_p
integer, dimension(-1:1^d &), target type_recv_srl_p1
integer, dimension(0:3^d &), target type_recv_p_f
integer, dimension( :^d &), pointer type_recv_srl
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
integer ixghi
Upper index of grid block arrays.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
integer, parameter bc_asymm
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
integer nwauxio
Number of auxiliary variables that are only included in output.
double precision unit_velocity
Physical scaling factor for velocity.
logical time_advance
do time evolving
character(len=std_len) restart_from_file
If not 'unavailable', resume from snapshot with this base file name.
double precision c_norm
Normalised speed of light.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
integer, parameter bc_cont
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter bc_symm
character(len= *), parameter undefined
logical need_global_cmax
need global maximal wave speed
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
character(len=std_len) base_filename
Base file name for simulation output, which will be followed by a number.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
logical record_electric_field
True for record electric field.
integer, parameter ixglo
Lower index of grid block arrays (always 1)
double precision, public mf_nu
viscosity coefficient in s cm^-2 for solar corona (Cheung 2012 ApJ)
logical, public, protected mf_divb_4thorder
Whether divB is computed with a fourth order approximation.
subroutine frictional_velocity(w, x, ixIL, ixOL)
subroutine mf_get_ct_velocity(vcts, wLp, wRp, ixIL, ixOL, idim, cmax, cmin)
prepare velocities for ct methods
subroutine add_source_hyperres(qdt, ixIL, ixOL, wCT, w, x)
Add Hyper-resistive source to w within ixO Uses 9 point stencil (4 neighbours) in each direction.
subroutine, public mf_to_primitive(ixIL, ixOL, w, x)
Transform conservative variables into primitive ones.
subroutine, public mf_phys_init()
subroutine mf_get_flux(wC, w, x, ixIL, ixOL, idim, f)
Calculate fluxes within ixO^L.
subroutine add_source_glm(qdt, ixIL, ixOL, wCT, w, x)
subroutine mf_check_params
subroutine mf_add_source_geom(qdt, dtfactor, ixIL, ixOL, wCT, wprim, w, x)
subroutine add_source_janhunen(qdt, ixIL, ixOL, wCT, w, x)
subroutine, public mf_face_to_center(ixOL, s)
calculate cell-center values from face-center values
subroutine mf_modify_wlr(ixIL, ixOL, qt, wLC, wRC, wLp, wRp, s, idir)
logical, public, protected mf_record_electric_field
set to true if need to record electric field on cell edges
subroutine add_source_linde(qdt, ixIL, ixOL, wCT, w, x)
subroutine update_faces_hll(ixIL, ixOL, qt, qdt, fE, sCT, s, vcts)
update faces
subroutine, public b_from_vector_potential(ixIsL, ixIL, ixOL, ws, x)
calculate magnetic field from vector potential
logical, public, protected clean_initial_divb
clean divb in the initial condition
subroutine fixdivb_boundary(ixGL, ixOL, w, x, iB)
logical, public divbwave
Add divB wave in Roe solver.
subroutine mf_get_cmax(w, x, ixIL, ixOL, idim, cmax)
Calculate cmax_idim=csound+abs(v_idim) within ixO^L.
double precision, public mf_vmax
maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
logical, public, protected mf_glm
Whether GLM-MHD is used.
subroutine get_resistive_electric_field(ixIL, ixOL, sCT, s, jce)
calculate eta J at cell edges
double precision, public mf_eta_hyper
The hyper-resistivity.
subroutine add_source_res2(qdt, ixIL, ixOL, wCT, w, x)
Add resistive source to w within ixO Uses 5 point stencil (2 neighbours) in each direction,...
subroutine, public get_normalized_divb(w, ixIL, ixOL, divb)
get dimensionless div B = |divB| * volume / area / |B|
subroutine add_source_res1(qdt, ixIL, ixOL, wCT, w, x)
Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in each direction,...
subroutine, public get_current(w, ixIL, ixOL, idirmin, current, fourthorder)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
subroutine mf_get_cbounds(wLC, wRC, wLp, wRp, x, ixIL, ixOL, idim, Hspeed, cmax, cmin)
Estimating bounds for the minimum and maximum signal velocities.
double precision function integral_grid_mf(ixIL, ixOL, w, x, iw, patchwi)
subroutine mf_velocity_update(qt, psa)
Add global source terms to update frictional velocity on complete domain.
double precision, public mf_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
subroutine mf_add_source(qdt, dtfactor, ixIL, ixOL, wCT, wCTprim, w, x, qsourcesplit, active)
w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
subroutine mf_get_dt(w, ixIL, ixOL, dtnew, dxD, x)
If resistivity is not zero, check diffusion time limit for dt.
subroutine mask_inner(ixIL, ixOL, w, x, patchwi, volume_pe)
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
subroutine update_faces_average(ixIL, ixOL, qt, qdt, fC, fE, sCT, s)
get electric field though averaging neighors to update faces in CT
character(len=std_len), public, protected type_ct
Method type of constrained transport.
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
subroutine, public mf_to_conserved(ixIL, ixOL, w, x)
Transform primitive variables into conservative ones.
subroutine mf_boundary_adjust(igrid, psb)
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
subroutine add_source_powel(qdt, ixIL, ixOL, wCT, w, x)
Add divB related sources to w within ixO corresponding to Powel.
double precision, dimension(2 *^nd), public mf_decay_scale
decay scale of frictional velocity near boundaries
subroutine, public record_force_free_metrics()
subroutine mf_write_info(fh)
Write this module's parameters to a snapsoht.
logical, public, protected mf_4th_order
MHD fourth order.
integer, public, protected psi_
Indices of the GLM psi.
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
subroutine, public mf_clean_divb_multigrid(qdt, qt, active)
logical, public, protected mf_particles
Whether particles module is added.
double precision, public mf_eta
The resistivity.
subroutine update_faces_contact(ixIL, ixOL, qt, qdt, wp, fC, fE, sCT, s, vcts)
update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
subroutine mf_update_faces(ixIL, ixOL, qt, qdt, wprim, fC, fE, sCT, s, vcts)
subroutine mf_physical_units()
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
Module containing all the particle routines.
subroutine particles_init()
Initialize particle data and parameters.
This module defines the procedures of a physics module. It contains function pointers for the various...
procedure(sub_get_ct_velocity), pointer phys_get_ct_velocity
procedure(sub_convert), pointer phys_to_primitive
integer, parameter flux_tvdlf
Indicates the flux should be treated with tvdlf.
procedure(sub_write_info), pointer phys_write_info
procedure(sub_get_flux), pointer phys_get_flux
logical phys_req_diagonal
Whether the physics routines require diagonal ghost cells, for example for computing a curl.
procedure(sub_get_cbounds), pointer phys_get_cbounds
procedure(sub_get_dt), pointer phys_get_dt
procedure(sub_add_source_geom), pointer phys_add_source_geom
procedure(sub_check_params), pointer phys_check_params
integer, parameter flux_default
Indicates a normal flux.
integer, dimension(:, :), allocatable flux_type
Array per direction per variable, which can be used to specify that certain fluxes have to be treated...
procedure(sub_update_faces), pointer phys_update_faces
procedure(sub_convert), pointer phys_to_conserved
character(len=name_len) physics_type
String describing the physics type of the simulation.
procedure(sub_clean_divb), pointer phys_clean_divb
procedure(sub_boundary_adjust), pointer phys_boundary_adjust
procedure(sub_global_source), pointer phys_global_source_after
procedure(sub_add_source), pointer phys_add_source
procedure(sub_face_to_center), pointer phys_face_to_center
procedure(sub_modify_wlr), pointer phys_modify_wlr
procedure(sub_get_cmax), pointer phys_get_cmax
logical phys_energy
Solve energy equation or not.
procedure(sub_special_advance), pointer phys_special_advance
Module with all the methods that users can customize in AMRVAC.
procedure(special_resistivity), pointer usr_special_resistivity
procedure(set_electric_field), pointer usr_set_electric_field
The data structure that contains information about a tree node/grid block.