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 :: ^
c&m^C_
46 integer,
public,
protected :: ^
c&b^C_
49 integer,
public,
protected ::
psi_
55 integer,
parameter :: divb_none = 0
56 integer,
parameter :: divb_multigrid = -1
57 integer,
parameter :: divb_glm = 1
58 integer,
parameter :: divb_powel = 2
59 integer,
parameter :: divb_janhunen = 3
60 integer,
parameter :: divb_linde = 4
61 integer,
parameter :: divb_lindejanhunen = 5
62 integer,
parameter :: divb_lindepowel = 6
63 integer,
parameter :: divb_lindeglm = 7
64 integer,
parameter :: divb_ct = 8
70 logical,
public,
protected ::
mf_glm = .false.
88 logical :: compactres = .false.
97 character(len=std_len),
public,
protected ::
typedivbfix =
'ct'
100 character(len=std_len),
public,
protected ::
type_ct =
'average'
103 character(len=std_len) :: typedivbdiff =
'all'
123 subroutine mf_read_params(files)
126 character(len=*),
intent(in) :: files(:)
137 do n = 1,
size(files)
138 open(
unitpar, file=trim(files(n)), status=
"old")
139 read(
unitpar, mf_list,
end=111)
143 end subroutine mf_read_params
146 subroutine mf_write_info(fh)
148 integer,
intent(in) :: fh
149 integer,
parameter :: n_par = 1
150 double precision :: values(n_par)
151 character(len=name_len) :: names(n_par)
152 integer,
dimension(MPI_STATUS_SIZE) :: st
155 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
159 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
160 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
161 end subroutine mf_write_info
182 type_divb = divb_none
185 type_divb = divb_multigrid
187 mg%operator_type = mg_laplacian
194 case (
'powel',
'powell')
195 type_divb = divb_powel
197 type_divb = divb_janhunen
199 type_divb = divb_linde
200 case (
'lindejanhunen')
201 type_divb = divb_lindejanhunen
203 type_divb = divb_lindepowel
207 type_divb = divb_lindeglm
212 call mpistop(
'Unknown divB fix')
217 mom(:) = var_set_momentum(
ndir)
222 mag(:) = var_set_bfield(
ndir)
227 allocate(start_indices(number_species),stop_indices(number_species))
229 start_indices(1)=mag(1)
232 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
241 stop_indices(1)=nwflux
247 allocate(iw_vector(nvector))
248 iw_vector(1) =
mom(1) - 1
249 iw_vector(2) = mag(1) - 1
256 call mpistop(
"phys_check error: flux_type has wrong shape")
275 if(type_divb==divb_glm)
then
289 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
308 call mf_physical_units()
312 subroutine mf_check_params
315 end subroutine mf_check_params
317 subroutine mf_physical_units()
319 double precision :: mp,kb,miu0,c_lightspeed
320 double precision :: a,
b
456 end subroutine mf_physical_units
461 integer,
intent(in) :: ixi^
l, ixo^
l
462 double precision,
intent(inout) :: w(ixi^s, nw)
463 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
471 integer,
intent(in) :: ixi^
l, ixo^
l
472 double precision,
intent(inout) :: w(ixi^s, nw)
473 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
479 subroutine mf_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
482 integer,
intent(in) :: ixi^
l, ixo^
l, idim
483 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
484 double precision,
intent(inout) :: cmax(ixi^s)
486 cmax(ixo^s)=abs(w(ixo^s,
mom(idim)))+one
488 end subroutine mf_get_cmax
491 subroutine mf_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
495 integer,
intent(in) :: ixi^
l, ixo^
l, idim
496 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
497 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
498 double precision,
intent(in) :: x(ixi^s,1:
ndim)
500 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
503 double precision,
dimension(ixI^S) :: tmp1
505 tmp1(ixo^s)=0.5d0*(wlc(ixo^s,
mom(idim))+wrc(ixo^s,
mom(idim)))
506 if(
present(cmin))
then
507 cmax(ixo^s,1)=max(tmp1(ixo^s)+one,zero)
508 cmin(ixo^s,1)=min(tmp1(ixo^s)-one,zero)
510 cmax(ixo^s,1)=abs(tmp1(ixo^s))+one
513 end subroutine mf_get_cbounds
516 subroutine mf_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
519 integer,
intent(in) :: ixi^
l, ixo^
l, idim
520 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
521 double precision,
intent(in) :: cmax(ixi^s)
522 double precision,
intent(in),
optional :: cmin(ixi^s)
523 type(ct_velocity),
intent(inout):: vcts
525 integer :: idime,idimn
531 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
533 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
535 if(.not.
allocated(vcts%vbarC))
then
536 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
537 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
540 if(
present(cmin))
then
541 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
542 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
544 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
545 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
548 idimn=mod(idim,
ndir)+1
549 idime=mod(idim+1,
ndir)+1
551 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
552 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
553 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
554 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
555 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
557 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
558 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
559 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
560 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
561 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
563 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
566 end subroutine mf_get_ct_velocity
569 subroutine mf_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
574 integer,
intent(in) :: ixi^
l, ixo^
l, idim
576 double precision,
intent(in) :: wc(ixi^s,nw)
578 double precision,
intent(in) :: w(ixi^s,nw)
579 double precision,
intent(in) :: x(ixi^s,1:
ndim)
580 double precision,
intent(out) :: f(ixi^s,nwflux)
584 {
do ix^db=ixomin^db,ixomax^db\}
590 {
do ix^db=ixomin^db,ixomax^db\}
591 f(ix^d,mag(idim))=w(ix^d,
psi_)
593 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
597 end subroutine mf_get_flux
600 subroutine mf_velocity_update(qt,psa)
603 double precision,
intent(in) :: qt
606 integer :: iigrid,igrid
607 logical :: stagger_flag
608 logical :: firstmf=.true.
623 do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
625 call frictional_velocity(psa(igrid)%w,psa(igrid)%x,ixg^
ll,
ixm^
ll)
649 end subroutine mf_velocity_update
652 subroutine mf_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
655 integer,
intent(in) :: ixi^
l, ixo^
l
656 double precision,
intent(in) :: qdt, dtfactor
657 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
658 double precision,
intent(inout) :: w(ixi^s,1:nw)
659 logical,
intent(in) :: qsourcesplit
660 logical,
intent(inout) :: active
662 if (.not. qsourcesplit)
then
665 if (abs(
mf_eta)>smalldouble)
then
667 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
672 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
680 select case (type_divb)
685 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
688 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
691 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
694 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
695 case (divb_lindejanhunen)
697 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
698 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
699 case (divb_lindepowel)
701 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
702 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
705 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
706 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
709 case (divb_multigrid)
712 call mpistop(
'Unknown divB fix')
717 end subroutine mf_add_source
719 subroutine frictional_velocity(w,x,ixI^L,ixO^L)
722 integer,
intent(in) :: ixi^
l, ixo^
l
723 double precision,
intent(in) :: x(ixi^s,1:
ndim)
724 double precision,
intent(inout) :: w(ixi^s,1:nw)
726 double precision :: xmin(
ndim),xmax(
ndim)
727 double precision :: decay(ixo^s)
728 double precision :: current(ixi^s,7-2*
ndir:3),tmp(ixo^s)
729 integer :: ix^
d, idirmin,idir,jdir,kdir
735 if(
block%is_physical_boundary(2*^
d-1))
then
738 if(2*^
d-1==5.and.
ndim==3)
then
740 current(ixomin^
d^
d%ixO^s,:)=current(ixomin^
d+1^
d%ixO^s,:)
745 current(ixomin^
d^
d%ixO^s,:)=current(ixomin^
d+1^
d%ixO^s,:)
749 if(
block%is_physical_boundary(2*^
d))
then
750 current(ixomax^
d^
d%ixO^s,:)=current(ixomax^
d-1^
d%ixO^s,:)
755 do idir=1,
ndir;
do jdir=idirmin,3;
do kdir=1,
ndir
756 if(
lvc(idir,jdir,kdir)/=0)
then
757 if(
lvc(idir,jdir,kdir)==1)
then
758 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+current(ixo^s,jdir)*w(ixo^s,mag(kdir))
760 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-current(ixo^s,jdir)*w(ixo^s,mag(kdir))
763 end do;
end do;
end do
765 tmp(ixo^s)=sum(w(ixo^s,mag(:))**2,dim=ndim+1)
767 where(tmp(ixo^s)/=0.d0)
768 tmp(ixo^s)=1.d0/(tmp(ixo^s)*
mf_nu)
772 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))*tmp(ixo^s)
776 ^d&xmin(^d)=xprobmin^d\
777 ^d&xmax(^d)=xprobmax^d\
787 tmp(ixo^s)=sqrt(sum(w(ixo^s,
mom(:))**2,dim=ndim+1))/
mf_vmax+1.d-12
788 tmp(ixo^s)=dtanh(tmp(ixo^s))/tmp(ixo^s)
790 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))*tmp(ixo^s)*decay(ixo^s)
793 end subroutine frictional_velocity
799 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
804 integer,
intent(in) :: ixi^
l, ixo^
l
805 double precision,
intent(in) :: qdt
806 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
807 double precision,
intent(inout) :: w(ixi^s,1:nw)
810 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
811 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
812 double precision :: tmp(ixi^s),tmp2(ixi^s)
813 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
814 integer :: lxo^
l, kxo^
l
823 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
824 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
831 gradeta(ixo^s,1:
ndim)=zero
837 gradeta(ixo^s,idim)=tmp(ixo^s)
841 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
847 tmp2(ixi^s)=bf(ixi^s,idir)
849 lxo^
l=ixo^
l+2*
kr(idim,^
d);
850 jxo^
l=ixo^
l+
kr(idim,^
d);
851 hxo^
l=ixo^
l-
kr(idim,^
d);
852 kxo^
l=ixo^
l-2*
kr(idim,^
d);
853 tmp(ixo^s)=tmp(ixo^s)+&
854 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
859 tmp2(ixi^s)=bf(ixi^s,idir)
861 jxo^
l=ixo^
l+
kr(idim,^
d);
862 hxo^
l=ixo^
l-
kr(idim,^
d);
863 tmp(ixo^s)=tmp(ixo^s)+&
864 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
869 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
873 do jdir=1,
ndim;
do kdir=idirmin,3
874 if (
lvc(idir,jdir,kdir)/=0)
then
875 if (
lvc(idir,jdir,kdir)==1)
then
876 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
878 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
885 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
889 end subroutine add_source_res1
893 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
898 integer,
intent(in) :: ixi^
l, ixo^
l
899 double precision,
intent(in) :: qdt
900 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
901 double precision,
intent(inout) :: w(ixi^s,1:nw)
904 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
905 double precision :: tmpvec(ixi^s,1:3)
906 integer :: ixa^
l,idir,idirmin,idirmin1
913 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
914 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
928 tmpvec(ixa^s,1:
ndir)=zero
930 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
937 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
940 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
943 end subroutine add_source_res2
947 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
951 integer,
intent(in) :: ixi^
l, ixo^
l
952 double precision,
intent(in) :: qdt
953 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
954 double precision,
intent(inout) :: w(ixi^s,1:nw)
956 double precision :: current(ixi^s,7-2*
ndir:3)
957 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
958 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
961 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
962 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
965 tmpvec(ixa^s,1:
ndir)=zero
967 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
971 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
974 tmpvec(ixa^s,1:
ndir)=zero
975 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
979 tmpvec2(ixa^s,1:
ndir)=zero
980 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
983 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
986 end subroutine add_source_hyperres
988 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
995 integer,
intent(in) :: ixi^
l, ixo^
l
996 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
997 double precision,
intent(inout) :: w(ixi^s,1:nw)
998 double precision:: divb(ixi^s)
999 double precision :: gradpsi(ixi^s)
1000 integer :: idim,idir
1027 end subroutine add_source_glm
1030 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
1033 integer,
intent(in) :: ixi^
l, ixo^
l
1034 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
1035 double precision,
intent(inout) :: w(ixi^s,1:nw)
1037 double precision :: divb(ixi^s)
1045 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,
mom(idir))*divb(ixo^s)
1048 end subroutine add_source_powel
1050 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
1055 integer,
intent(in) :: ixi^
l, ixo^
l
1056 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
1057 double precision,
intent(inout) :: w(ixi^s,1:nw)
1058 double precision :: divb(ixi^s)
1066 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,
mom(idir))*divb(ixo^s)
1069 end subroutine add_source_janhunen
1071 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
1076 integer,
intent(in) :: ixi^
l, ixo^
l
1077 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
1078 double precision,
intent(inout) :: w(ixi^s,1:nw)
1079 double precision :: divb(ixi^s),graddivb(ixi^s)
1080 integer :: idim, idir, ixp^
l, i^
d, iside
1081 logical,
dimension(-1:1^D&) :: leveljump
1089 if(i^
d==0|.and.) cycle
1090 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
1091 leveljump(i^
d)=.true.
1093 leveljump(i^
d)=.false.
1102 i^dd=kr(^dd,^d)*(2*iside-3);
1103 if (leveljump(i^dd))
then
1105 ixpmin^d=ixomin^d-i^d
1107 ixpmax^d=ixomax^d-i^d
1118 select case(typegrad)
1120 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
1122 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
1126 if (slab_uniform)
then
1127 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
1129 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
1130 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
1133 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
1136 end subroutine add_source_linde
1144 integer,
intent(in) :: ixi^
l, ixo^
l
1145 double precision,
intent(in) :: w(ixi^s,1:nw)
1146 double precision :: divb(ixi^s), dsurface(ixi^s)
1148 double precision :: invb(ixo^s)
1149 integer :: ixa^
l,idims
1151 call get_divb(w,ixi^
l,ixo^
l,divb)
1152 invb(ixo^s)=sqrt(sum(w(ixo^s, mag(:))**2, dim=
ndim+1))
1153 where(invb(ixo^s)/=0.d0)
1154 invb(ixo^s)=1.d0/invb(ixo^s)
1157 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
1159 ixamin^
d=ixomin^
d-1;
1160 ixamax^
d=ixomax^
d-1;
1161 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
1163 ixa^
l=ixo^
l-
kr(idims,^
d);
1164 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
1166 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
1167 block%dvolume(ixo^s)/dsurface(ixo^s)
1178 integer,
intent(in) :: ixo^
l, ixi^
l
1179 double precision,
intent(in) :: w(ixi^s,1:nw)
1180 integer,
intent(out) :: idirmin
1181 logical,
intent(in),
optional :: fourthorder
1184 double precision :: current(ixi^s,7-2*
ndir:3)
1185 integer :: idir, idirmin0
1189 if(
present(fourthorder))
then
1198 subroutine mf_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
1202 integer,
intent(in) :: ixi^
l, ixo^
l
1203 double precision,
intent(inout) :: dtnew
1204 double precision,
intent(in) ::
dx^
d
1205 double precision,
intent(in) :: w(ixi^s,1:nw)
1206 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1208 double precision :: dxarr(
ndim)
1209 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
1210 integer :: idirmin,idim
1217 else if (
mf_eta<zero)
then
1224 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
1227 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
1239 end subroutine mf_get_dt
1242 subroutine mf_add_source_geom(qdt,dtfactor, ixI^L,ixO^L,wCT,wprim,w,x)
1246 integer,
intent(in) :: ixi^
l, ixo^
l
1247 double precision,
intent(in) :: qdt, dtfactor, x(ixi^s,1:
ndim)
1248 double precision,
intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
1250 double precision :: tmp,invr,cs2
1251 integer :: iw,idir,ix^
d
1253 integer :: mr_,mphi_
1254 integer :: br_,bphi_
1257 br_=mag(1); bphi_=mag(1)-1+
phi_
1263 {
do ix^db=ixomin^db,ixomax^db\}
1264 w(ix^
d,bphi_)=w(ix^
d,bphi_)+qdt/x(ix^
d,1)*&
1265 (wct(ix^
d,bphi_)*wct(ix^
d,mr_) &
1266 -wct(ix^
d,br_)*wct(ix^
d,mphi_))
1270 if(
mf_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
1272 {
do ix^db=ixomin^db,ixomax^db\}
1276 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wct(ix^d,
psi_)
1281 if(.not.stagger_grid)
then
1282 tmp=wct(ix^d,
mom(1))*wct(ix^d,mag(2))-wct(ix^d,
mom(2))*wct(ix^d,mag(1))
1283 cs2=1.d0/tan(x(ix^d,2))
1285 tmp=tmp+cs2*wct(ix^d,
psi_)
1287 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
1294 if(.not.stagger_grid)
then
1295 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
1296 (wct(ix^d,
mom(1))*wct(ix^d,mag(3)) &
1297 -wct(ix^d,
mom(3))*wct(ix^d,mag(1)))
1302 if(.not.stagger_grid)
then
1303 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
1304 (wct(ix^d,
mom(1))*wct(ix^d,mag(3)) &
1305 -wct(ix^d,
mom(3))*wct(ix^d,mag(1)) &
1306 -(wct(ix^d,
mom(3))*wct(ix^d,mag(2)) &
1307 -wct(ix^d,
mom(2))*wct(ix^d,mag(3)))*cs2)
1314 end subroutine mf_add_source_geom
1316 subroutine mf_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
1319 integer,
intent(in) :: ixi^
l, ixo^
l, idir
1320 double precision,
intent(in) :: qt
1321 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
1322 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
1324 double precision :: db(ixi^s), dpsi(ixi^s)
1332 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
1333 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
1335 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
1337 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
1340 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
1343 end subroutine mf_modify_wlr
1345 subroutine mf_boundary_adjust(igrid,psb)
1347 integer,
intent(in) :: igrid
1350 integer :: ib, idims, iside, ixo^
l, i^
d
1359 i^
d=
kr(^
d,idims)*(2*iside-3);
1360 if (neighbor_type(i^
d,igrid)/=1) cycle
1361 ib=(idims-1)*2+iside
1379 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
1384 end subroutine mf_boundary_adjust
1386 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
1389 integer,
intent(in) :: ixg^
l,ixo^
l,ib
1390 double precision,
intent(inout) :: w(ixg^s,1:nw)
1391 double precision,
intent(in) :: x(ixg^s,1:
ndim)
1393 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
1394 integer :: ix^
d,ixf^
l
1406 do ix1=ixfmax1,ixfmin1,-1
1407 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
1408 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
1409 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
1412 do ix1=ixfmax1,ixfmin1,-1
1413 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
1414 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
1415 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
1416 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
1417 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
1418 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
1419 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
1433 do ix1=ixfmax1,ixfmin1,-1
1434 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1435 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
1436 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
1437 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
1438 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
1439 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
1442 do ix1=ixfmax1,ixfmin1,-1
1443 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1444 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
1445 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
1446 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
1447 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
1448 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
1449 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
1450 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
1451 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
1452 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
1453 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
1454 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
1455 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
1456 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
1457 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1458 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
1459 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
1460 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
1472 do ix1=ixfmin1,ixfmax1
1473 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
1474 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
1475 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
1478 do ix1=ixfmin1,ixfmax1
1479 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
1480 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
1481 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
1482 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
1483 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
1484 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
1485 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
1499 do ix1=ixfmin1,ixfmax1
1500 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1501 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
1502 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
1503 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
1504 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
1505 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
1508 do ix1=ixfmin1,ixfmax1
1509 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1510 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
1511 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
1512 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
1513 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
1514 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
1515 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
1516 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
1517 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
1518 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
1519 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
1520 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
1521 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
1522 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
1523 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1524 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
1525 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
1526 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
1538 do ix2=ixfmax2,ixfmin2,-1
1539 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
1540 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
1541 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
1544 do ix2=ixfmax2,ixfmin2,-1
1545 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
1546 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
1547 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
1548 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
1549 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
1550 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
1551 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
1565 do ix2=ixfmax2,ixfmin2,-1
1566 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
1567 ix2+1,ixfmin3:ixfmax3,mag(2)) &
1568 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
1569 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
1570 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
1571 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
1574 do ix2=ixfmax2,ixfmin2,-1
1575 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
1576 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
1577 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
1578 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
1579 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
1580 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1581 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
1582 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
1583 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1584 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
1585 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
1586 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
1587 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
1588 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
1589 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1590 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
1591 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
1592 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
1604 do ix2=ixfmin2,ixfmax2
1605 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
1606 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
1607 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
1610 do ix2=ixfmin2,ixfmax2
1611 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
1612 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
1613 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
1614 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
1615 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
1616 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
1617 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
1631 do ix2=ixfmin2,ixfmax2
1632 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
1633 ix2-1,ixfmin3:ixfmax3,mag(2)) &
1634 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
1635 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
1636 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
1637 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
1640 do ix2=ixfmin2,ixfmax2
1641 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
1642 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
1643 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
1644 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
1645 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
1646 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1647 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
1648 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
1649 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1650 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
1651 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
1652 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
1653 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
1654 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
1655 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1656 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
1657 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
1658 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
1673 do ix3=ixfmax3,ixfmin3,-1
1674 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
1675 ixfmin2:ixfmax2,ix3+1,mag(3)) &
1676 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
1677 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
1678 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
1679 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
1682 do ix3=ixfmax3,ixfmin3,-1
1683 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
1684 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
1685 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
1686 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
1687 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
1688 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1689 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
1690 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
1691 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1692 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
1693 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
1694 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
1695 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
1696 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
1697 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
1698 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
1699 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
1700 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
1713 do ix3=ixfmin3,ixfmax3
1714 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
1715 ixfmin2:ixfmax2,ix3-1,mag(3)) &
1716 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
1717 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
1718 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
1719 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
1722 do ix3=ixfmin3,ixfmax3
1723 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
1724 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
1725 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
1726 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
1727 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
1728 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1729 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
1730 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
1731 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1732 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
1733 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
1734 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
1735 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
1736 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
1737 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
1738 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
1739 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
1740 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
1745 call mpistop(
"Special boundary is not defined for this region")
1748 end subroutine fixdivb_boundary
1757 double precision,
intent(in) :: qdt
1758 double precision,
intent(in) :: qt
1759 logical,
intent(inout) :: active
1761 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
1762 double precision :: res
1763 double precision,
parameter :: max_residual = 1
d-3
1764 double precision,
parameter :: residual_reduction = 1
d-10
1766 integer,
parameter :: max_its = 50
1767 double precision :: residual_it(max_its), max_divb
1768 integer :: iigrid, igrid
1769 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
1772 mg%operator_type = mg_laplacian
1780 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1781 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1784 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
1785 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1787 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1788 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1791 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
1792 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1796 write(*,*)
"mf_clean_divb_multigrid warning: unknown boundary type"
1797 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1798 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1806 do iigrid = 1, igridstail
1807 igrid = igrids(iigrid);
1810 lvl =
mg%boxes(id)%lvl
1811 nc =
mg%box_size_lvl(lvl)
1817 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
1819 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
1820 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
1825 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
1828 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
1829 if (
mype == 0) print *,
"iteration vs residual"
1832 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
1833 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
1834 if (residual_it(n) < residual_reduction * max_divb)
exit
1836 if (
mype == 0 .and. n > max_its)
then
1837 print *,
"divb_multigrid warning: not fully converged"
1838 print *,
"current amplitude of divb: ", residual_it(max_its)
1839 print *,
"multigrid smallest grid: ", &
1840 mg%domain_size_lvl(:,
mg%lowest_lvl)
1841 print *,
"note: smallest grid ideally has <= 8 cells"
1842 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
1843 print *,
"note: dx/dy/dz should be similar"
1847 call mg_fas_vcycle(
mg, max_res=res)
1848 if (res < max_residual)
exit
1850 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
1855 do iigrid = 1, igridstail
1856 igrid = igrids(iigrid);
1865 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
1869 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
1871 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
1872 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
1880 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
1881 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
1891 subroutine mf_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
1894 integer,
intent(in) :: ixi^
l, ixo^
l
1895 double precision,
intent(in) :: qt, qdt
1897 double precision,
intent(in) :: wprim(ixi^s,1:nw)
1898 type(state) :: sct, s
1899 type(ct_velocity) :: vcts
1900 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
1901 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
1905 call update_faces_average(ixi^
l,ixo^
l,qt,qdt,fc,fe,sct,s)
1907 call update_faces_contact(ixi^
l,ixo^
l,qt,qdt,wprim,fc,fe,sct,s,vcts)
1909 call update_faces_hll(ixi^
l,ixo^
l,qt,qdt,fe,sct,s,vcts)
1911 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
1914 end subroutine mf_update_faces
1917 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
1921 integer,
intent(in) :: ixi^
l, ixo^
l
1922 double precision,
intent(in) :: qt, qdt
1923 type(state) :: sct, s
1924 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
1925 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
1927 double precision :: circ(ixi^s,1:
ndim)
1928 double precision :: e_resi(ixi^s,
sdim:3)
1929 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
1930 integer :: idim1,idim2,idir,iwdim1,iwdim2
1932 associate(bfaces=>s%ws,x=>s%x)
1939 if(
mf_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
1943 i1kr^
d=
kr(idim1,^
d);
1946 i2kr^
d=
kr(idim2,^
d);
1949 if (
lvc(idim1,idim2,idir)==1)
then
1951 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
1953 {
do ix^db=ixcmin^db,ixcmax^db\}
1954 fe(ix^
d,idir)=quarter*&
1955 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
1956 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
1958 if(
mf_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
1961 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
1969 if(
associated(usr_set_electric_field)) &
1970 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
1972 circ(ixi^s,1:ndim)=zero
1978 ixcmin^d=ixomin^d-kr(idim1,^d);
1980 ixa^l=ixc^l-kr(idim2,^d);
1983 if(lvc(idim1,idim2,idir)==1)
then
1985 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
1988 else if(lvc(idim1,idim2,idir)==-1)
then
1990 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
1997 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
1999 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2004 end subroutine update_faces_average
2007 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
2011 integer,
intent(in) :: ixi^
l, ixo^
l
2012 double precision,
intent(in) :: qt, qdt
2014 double precision,
intent(in) :: wp(ixi^s,1:nw)
2015 type(state) :: sct, s
2016 type(ct_velocity) :: vcts
2017 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
2018 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
2020 double precision :: circ(ixi^s,1:
ndim)
2022 double precision :: ecc(ixi^s,
sdim:3)
2024 double precision :: el(ixi^s),er(ixi^s)
2026 double precision :: elc(ixi^s),erc(ixi^s)
2028 double precision :: jce(ixi^s,
sdim:3)
2029 integer :: hxc^
l,ixc^
l,jxc^
l,ixa^
l,ixb^
l
2030 integer :: idim1,idim2,idir,iwdim1,iwdim2
2032 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
2037 if(
lvc(idim1,idim2,idir)==1)
then
2038 ecc(ixi^s,idir)=ecc(ixi^s,idir)+wp(ixi^s,mag(idim1))*wp(ixi^s,
mom(idim2))
2039 else if(
lvc(idim1,idim2,idir)==-1)
then
2040 ecc(ixi^s,idir)=ecc(ixi^s,idir)-wp(ixi^s,mag(idim1))*wp(ixi^s,
mom(idim2))
2045 if(
mf_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,jce)
2057 if (
lvc(idim1,idim2,idir)==1)
then
2059 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
2061 jxc^
l=ixc^
l+
kr(idim1,^
d);
2062 hxc^
l=ixc^
l+
kr(idim2,^
d);
2064 fe(ixc^s,idir)=quarter*&
2065 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
2066 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
2070 ixamax^
d=ixcmax^
d+
kr(idim1,^
d);
2071 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
2072 hxc^
l=ixa^
l+
kr(idim2,^
d);
2073 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
2074 where(vnorm(ixc^s,idim1)>0.d0)
2075 elc(ixc^s)=el(ixc^s)
2076 else where(vnorm(ixc^s,idim1)<0.d0)
2077 elc(ixc^s)=el(jxc^s)
2079 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
2081 hxc^
l=ixc^
l+
kr(idim2,^
d);
2082 where(vnorm(hxc^s,idim1)>0.d0)
2083 erc(ixc^s)=er(ixc^s)
2084 else where(vnorm(hxc^s,idim1)<0.d0)
2085 erc(ixc^s)=er(jxc^s)
2087 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
2089 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
2092 jxc^
l=ixc^
l+
kr(idim2,^
d);
2094 ixamax^
d=ixcmax^
d+
kr(idim2,^
d);
2095 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
2096 hxc^
l=ixa^
l+
kr(idim1,^
d);
2097 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
2098 where(vnorm(ixc^s,idim2)>0.d0)
2099 elc(ixc^s)=el(ixc^s)
2100 else where(vnorm(ixc^s,idim2)<0.d0)
2101 elc(ixc^s)=el(jxc^s)
2103 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
2105 hxc^
l=ixc^
l+
kr(idim1,^
d);
2106 where(vnorm(hxc^s,idim2)>0.d0)
2107 erc(ixc^s)=er(ixc^s)
2108 else where(vnorm(hxc^s,idim2)<0.d0)
2109 erc(ixc^s)=er(jxc^s)
2111 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
2113 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
2116 if(
mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2120 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
2135 circ(ixi^s,1:
ndim)=zero
2140 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
2144 if(
lvc(idim1,idim2,idir)/=0)
then
2145 hxc^
l=ixc^
l-
kr(idim2,^
d);
2147 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2148 +
lvc(idim1,idim2,idir)&
2155 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
2156 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2158 circ(ixc^s,idim1)=zero
2161 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2166 end subroutine update_faces_contact
2169 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
2174 integer,
intent(in) :: ixi^
l, ixo^
l
2175 double precision,
intent(in) :: qt, qdt
2176 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
2177 type(state) :: sct, s
2178 type(ct_velocity) :: vcts
2180 double precision :: vtill(ixi^s,2)
2181 double precision :: vtilr(ixi^s,2)
2182 double precision :: btill(ixi^s,
ndim)
2183 double precision :: btilr(ixi^s,
ndim)
2184 double precision :: cp(ixi^s,2)
2185 double precision :: cm(ixi^s,2)
2186 double precision :: circ(ixi^s,1:
ndim)
2188 double precision :: jce(ixi^s,
sdim:3)
2189 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
2190 integer :: idim1,idim2,idir
2192 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
2193 cbarmax=>vcts%cbarmax)
2206 if(
mf_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,jce)
2219 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
2223 idim2=mod(idir+1,3)+1
2225 jxc^
l=ixc^
l+
kr(idim1,^
d);
2226 ixcp^
l=ixc^
l+
kr(idim2,^
d);
2230 vtill(ixi^s,2),vtilr(ixi^s,2))
2233 vtill(ixi^s,1),vtilr(ixi^s,1))
2239 btill(ixi^s,idim1),btilr(ixi^s,idim1))
2242 btill(ixi^s,idim2),btilr(ixi^s,idim2))
2246 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
2247 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
2249 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
2250 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
2254 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
2255 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
2256 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
2257 /(cp(ixc^s,1)+cm(ixc^s,1)) &
2258 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
2259 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
2260 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
2261 /(cp(ixc^s,2)+cm(ixc^s,2))
2264 if(
mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2268 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
2282 circ(ixi^s,1:
ndim)=zero
2288 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
2292 if(
lvc(idim1,idim2,idir)/=0)
then
2293 hxc^
l=ixc^
l-
kr(idim2,^
d);
2295 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2296 +
lvc(idim1,idim2,idir)&
2303 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
2304 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2306 circ(ixc^s,idim1)=zero
2309 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2313 end subroutine update_faces_hll
2316 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
2321 integer,
intent(in) :: ixi^
l, ixo^
l
2322 type(state),
intent(in) :: sct, s
2324 double precision :: jce(ixi^s,
sdim:3)
2327 double precision :: jcc(ixi^s,7-2*
ndir:3)
2329 double precision :: xs(ixgs^t,1:
ndim)
2331 double precision :: eta(ixi^s)
2332 double precision :: gradi(ixgs^t)
2333 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
2335 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
2341 if (
lvc(idim1,idim2,idir)==0) cycle
2342 ixcmax^
d=ixomax^
d+1;
2343 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-2;
2344 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
2347 xs(ixb^s,:)=x(ixb^s,:)
2348 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
2349 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi)
2350 if (
lvc(idim1,idim2,idir)==1)
then
2351 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
2353 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
2360 jce(ixi^s,:)=jce(ixi^s,:)*
mf_eta
2368 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
2369 jcc(ixc^s,idir)=0.d0
2371 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
2372 ixamin^
d=ixcmin^
d+ix^
d;
2373 ixamax^
d=ixcmax^
d+ix^
d;
2374 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
2376 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
2377 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
2382 end subroutine get_resistive_electric_field
2388 integer,
intent(in) :: ixo^
l
2395 do ix^db=ixomin^db,ixomax^db\}
2397 s%w(ix^
d,mag(1))=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
2398 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
2399 s%w(ix^
d,mag(2))=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
2400 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
2401 s%w(ix^
d,mag(3))=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
2402 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
2405 s%w(ix^
d,mag(1))=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
2406 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
2407 s%w(ix^
d,mag(2))=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
2408 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
2451 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
2452 double precision,
intent(inout) :: ws(ixis^s,1:nws)
2453 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2455 double precision :: adummy(ixis^s,1:3)
2464 double precision :: sum_jbb_ipe, sum_j_ipe, sum_l_ipe, f_i_ipe, volume_pe
2465 double precision :: sum_jbb, sum_j, sum_l, f_i, volume, cw_sin_theta
2466 integer :: iigrid, igrid, ix^
d
2467 integer :: amode, istatus(mpi_status_size)
2468 integer,
save :: fhmf
2469 logical :: patchwi(ixg^t)
2470 logical,
save :: logmfopened=.false.
2471 character(len=800) :: filename,filehead
2472 character(len=800) :: line,datastr
2479 do iigrid=1,igridstail; igrid=igrids(iigrid);
2482 call mask_inner(ixg^
ll,
ixm^
ll,ps(igrid)%w,ps(igrid)%x,patchwi,volume_pe)
2483 sum_jbb_ipe = sum_jbb_ipe+integral_grid_mf(ixg^
ll,
ixm^
ll,ps(igrid)%w,&
2484 ps(igrid)%x,1,patchwi)
2485 sum_j_ipe = sum_j_ipe+integral_grid_mf(ixg^
ll,
ixm^
ll,ps(igrid)%w,&
2486 ps(igrid)%x,2,patchwi)
2487 f_i_ipe=f_i_ipe+integral_grid_mf(ixg^
ll,
ixm^
ll,ps(igrid)%w,&
2488 ps(igrid)%x,3,patchwi)
2489 sum_l_ipe = sum_l_ipe+integral_grid_mf(ixg^
ll,
ixm^
ll,ps(igrid)%w,&
2490 ps(igrid)%x,4,patchwi)
2492 call mpi_reduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
2494 call mpi_reduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,0,&
2496 call mpi_reduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,0,&
2498 call mpi_reduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,0,&
2500 call mpi_reduce(volume_pe,volume,1,mpi_double_precision,mpi_sum,0,&
2506 cw_sin_theta = sum_jbb/sum_j
2512 if(.not.logmfopened)
then
2514 write(filename,
"(a,a)") trim(
base_filename),
"_mf_metrics.csv"
2518 open(unit=20,file=trim(filename),status=
'replace')
2519 close(20, status=
'delete')
2522 amode=ior(mpi_mode_create,mpi_mode_wronly)
2523 amode=ior(amode,mpi_mode_append)
2524 call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,
ierrmpi)
2526 filehead=
" it,time,CW_sin_theta,current,Lorenz_force,f_i"
2528 call mpi_file_write(fhmf,filehead,len_trim(filehead), &
2529 mpi_character,istatus,
ierrmpi)
2530 call mpi_file_write(fhmf,achar(10),1,mpi_character,istatus,
ierrmpi)
2534 write(datastr,
'(i6,a)')
it,
','
2535 line=trim(line)//trim(datastr)
2537 line=trim(line)//trim(datastr)
2538 write(datastr,
'(es13.6,a)') cw_sin_theta,
','
2539 line=trim(line)//trim(datastr)
2540 write(datastr,
'(es13.6,a)') sum_j,
','
2541 line=trim(line)//trim(datastr)
2542 write(datastr,
'(es13.6,a)') sum_l,
','
2543 line=trim(line)//trim(datastr)
2544 write(datastr,
'(es13.6)') f_i
2545 line=trim(line)//trim(datastr)//new_line(
'A')
2546 call mpi_file_write(fhmf,line,len_trim(line),mpi_character,istatus,
ierrmpi)
2548 call mpi_file_close(fhmf,
ierrmpi)
2555 subroutine mask_inner(ixI^L,ixO^L,w,x,patchwi,volume_pe)
2558 integer,
intent(in) :: ixi^
l,ixo^
l
2559 double precision,
intent(in):: w(ixi^s,nw),x(ixi^s,1:
ndim)
2560 double precision,
intent(inout) :: volume_pe
2561 logical,
intent(out) :: patchwi(ixi^s)
2563 double precision :: xo^
l
2566 {xomin^
d = xprobmin^
d + 0.05d0*(xprobmax^
d-xprobmin^
d)\}
2567 {xomax^
d = xprobmax^
d - 0.05d0*(xprobmax^
d-xprobmin^
d)\}
2569 xomin^nd = xprobmin^nd
2574 {
do ix^db=ixomin^db,ixomax^db\}
2575 if({ x(ix^dd,^
d) > xomin^
d .and. x(ix^dd,^
d) < xomax^
d | .and. })
then
2576 patchwi(ix^
d)=.true.
2577 volume_pe=volume_pe+
block%dvolume(ix^
d)
2579 patchwi(ix^
d)=.false.
2583 end subroutine mask_inner
2585 function integral_grid_mf(ixI^L,ixO^L,w,x,iw,patchwi)
2588 integer,
intent(in) :: ixi^
l,ixo^
l,iw
2589 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2590 double precision :: w(ixi^s,nw+
nwauxio)
2591 logical,
intent(in) :: patchwi(ixi^s)
2593 double precision,
dimension(ixI^S,7-2*ndir:3) :: current
2594 double precision,
dimension(ixI^S,1:ndir) :: bvec,qvec
2595 double precision :: integral_grid_mf,tmp(ixi^s),bm2
2596 integer :: ix^
d,idirmin0,idirmin,idir,jdir,kdir
2598 integral_grid_mf=0.d0
2602 bvec(ixo^s,:)=w(ixo^s,mag(:))
2605 qvec(ixo^s,1:
ndir)=zero
2606 do idir=1,
ndir;
do jdir=idirmin,3;
do kdir=1,
ndir
2607 if(
lvc(idir,jdir,kdir)/=0)
then
2608 if(
lvc(idir,jdir,kdir)==1)
then
2609 qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
2611 qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
2614 end do;
end do;
end do
2616 {
do ix^db=ixomin^db,ixomax^db\}
2617 if(patchwi(ix^d))
then
2618 bm2=sum(bvec(ix^d,:)**2)
2619 if(bm2/=0.d0) bm2=1.d0/bm2
2620 integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2)*&
2621 bm2)*block%dvolume(ix^d)
2627 {
do ix^db=ixomin^db,ixomax^db\}
2628 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(current(ix^d,:)**2))*&
2635 {
do ix^db=ixomin^db,ixomax^db\}
2636 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+tmp(ix^d)*block%dvolume(ix^d)
2640 bvec(ixo^s,:)=w(ixo^s,mag(:))
2643 qvec(ixo^s,1:ndir)=zero
2644 do idir=1,ndir;
do jdir=idirmin,3;
do kdir=1,ndir
2645 if(lvc(idir,jdir,kdir)/=0)
then
2646 if(lvc(idir,jdir,kdir)==1)
then
2647 qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
2649 qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
2652 end do;
end do;
end do
2654 {
do ix^db=ixomin^db,ixomax^db\}
2655 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+&
2656 sqrt(sum(qvec(ix^d,:)**2))*block%dvolume(ix^d)
2660 end function integral_grid_mf
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
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...
subroutine b_from_vector_potentiala(ixisl, ixil, ixol, ws, x, a)
calculate magnetic field from vector potential A at cell edges
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.
subroutine, public get_divb(w, ixil, ixol, divb, nth_in)
Calculate div B within ixO.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
Module with geometry-related routines (e.g., divergence, curl)
integer, parameter cylindrical
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 gradient(q, ixil, ixol, idir, gradq, nth_in)
subroutine gradientf(q, x, ixil, ixol, idir, gradq, nth_in, pm_in)
subroutine gradientl(q, ixil, ixol, idir, gradq)
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
subroutine getbc(time, qdt, psb, nwstart, nwbc)
do update ghost cells of all blocks including physical boundaries
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
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 sdim
starting dimension for electric field
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 max_blocks
The maximum number of grid blocks in a processor.
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)
integer, public, protected b
double precision, public mf_nu
viscosity coefficient in s cm^-2 for solar corona (Cheung 2012 ApJ)
subroutine, public mf_phys_init()
integer, public, protected mf_divb_nth
Whether divB is computed with a fourth order approximation.
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
logical, public, protected mf_record_electric_field
set to true if need to record electric field on cell edges
logical, public, protected clean_initial_divb
clean divb in the initial condition
logical, public divbwave
Add divB wave in Roe solver.
subroutine, public mf_face_to_center(ixol, s)
calculate cell-center values from face-center values
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.
double precision, public mf_eta_hyper
The hyper-resistivity.
double precision, public mf_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
integer, public, protected c
Indices of the momentum density for the form of better vectorization.
integer, public, protected c_
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
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 get_normalized_divb(w, ixil, ixol, divb)
get dimensionless div B = |divB| * volume / area / |B|
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
double precision, dimension(2 *^nd), public mf_decay_scale
decay scale of frictional velocity near boundaries
subroutine, public record_force_free_metrics()
subroutine, public mf_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
logical, public, protected mf_4th_order
MHD fourth order.
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)...
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.
integer, public, protected m
subroutine, public mf_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
double precision, public mf_eta
The resistivity.
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.
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
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
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
The data structure that contains information about a tree node/grid block.