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
140 subroutine mf_write_info(fh)
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)
155 end subroutine mf_write_info
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')
211 mom(:) = var_set_momentum(
ndir)
215 mag(:) = var_set_bfield(
ndir)
219 allocate(start_indices(number_species),stop_indices(number_species))
221 start_indices(1)=mag(1)
224 psi_ = var_set_fluxvar(
'psi',
'psi', need_bc=.false.)
233 stop_indices(1)=nwflux
239 allocate(iw_vector(nvector))
240 iw_vector(1) =
mom(1) - 1
241 iw_vector(2) = mag(1) - 1
248 call mpistop(
"phys_check error: flux_type has wrong shape")
267 if(type_divb==divb_glm)
then
281 write(*,*)
'*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
300 call mf_physical_units()
304 subroutine mf_check_params
307 end subroutine mf_check_params
309 subroutine mf_physical_units()
311 double precision :: mp,kb,miu0,c_lightspeed
366 end subroutine mf_physical_units
371 integer,
intent(in) :: ixi^
l, ixo^
l
372 double precision,
intent(inout) :: w(ixi^s, nw)
373 double precision,
intent(in) :: x(ixi^s, 1:
ndim)
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)
389 subroutine mf_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
392 integer,
intent(in) :: ixi^
l, ixo^
l, idim
393 double precision,
intent(in) :: w(ixi^s, nw), x(ixi^s,1:
ndim)
394 double precision,
intent(inout) :: cmax(ixi^s)
396 cmax(ixo^s)=abs(w(ixo^s,
mom(idim)))+one
398 end subroutine mf_get_cmax
401 subroutine mf_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
405 integer,
intent(in) :: ixi^
l, ixo^
l, idim
406 double precision,
intent(in) :: wlc(ixi^s,
nw), wrc(ixi^s,
nw)
407 double precision,
intent(in) :: wlp(ixi^s,
nw), wrp(ixi^s,
nw)
408 double precision,
intent(in) :: x(ixi^s,1:
ndim)
410 double precision,
intent(inout),
optional :: cmin(ixi^s,1:
number_species)
413 double precision,
dimension(ixI^S) :: tmp1
415 tmp1(ixo^s)=0.5d0*(wlc(ixo^s,
mom(idim))+wrc(ixo^s,
mom(idim)))
416 if(
present(cmin))
then
417 cmax(ixo^s,1)=max(tmp1(ixo^s)+one,zero)
418 cmin(ixo^s,1)=min(tmp1(ixo^s)-one,zero)
420 cmax(ixo^s,1)=abs(tmp1(ixo^s))+one
423 end subroutine mf_get_cbounds
426 subroutine mf_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
429 integer,
intent(in) :: ixi^
l, ixo^
l, idim
430 double precision,
intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
431 double precision,
intent(in) :: cmax(ixi^s)
432 double precision,
intent(in),
optional :: cmin(ixi^s)
433 type(ct_velocity),
intent(inout):: vcts
435 integer :: idime,idimn
441 if(.not.
allocated(vcts%vnorm))
allocate(vcts%vnorm(ixi^s,1:
ndim))
443 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,
mom(idim))+wrp(ixo^s,
mom(idim)))
445 if(.not.
allocated(vcts%vbarC))
then
446 allocate(vcts%vbarC(ixi^s,1:
ndir,2),vcts%vbarLC(ixi^s,1:
ndir,2),vcts%vbarRC(ixi^s,1:
ndir,2))
447 allocate(vcts%cbarmin(ixi^s,1:
ndim),vcts%cbarmax(ixi^s,1:
ndim))
450 if(
present(cmin))
then
451 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
452 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
454 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
455 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
458 idimn=mod(idim,
ndir)+1
459 idime=mod(idim+1,
ndir)+1
461 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,
mom(idimn))
462 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,
mom(idimn))
463 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
464 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
465 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
467 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,
mom(idime))
468 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,
mom(idime))
469 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
470 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
471 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
473 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
476 end subroutine mf_get_ct_velocity
479 subroutine mf_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
484 integer,
intent(in) :: ixi^
l, ixo^
l, idim
486 double precision,
intent(in) :: wc(ixi^s,nw)
488 double precision,
intent(in) :: w(ixi^s,nw)
489 double precision,
intent(in) :: x(ixi^s,1:
ndim)
490 double precision,
intent(out) :: f(ixi^s,nwflux)
492 integer :: idir, ix^
d
495 {
do ix^db=ixomin^db,ixomax^db\}
498 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))
502 {
do ix^db=ixomin^db,ixomax^db\}
503 f(ix^d,mag(idim))=w(ix^d,
psi_)
505 f(ix^d,
psi_)=cmax_global**2*w(ix^d,mag(idim))
509 end subroutine mf_get_flux
512 subroutine mf_velocity_update(qt,psa)
515 double precision,
intent(in) :: qt
518 integer :: iigrid,igrid
519 logical :: stagger_flag
520 logical :: firstmf=.true.
535 do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
537 call frictional_velocity(psa(igrid)%w,psa(igrid)%x,ixg^
ll,
ixm^
ll)
561 end subroutine mf_velocity_update
564 subroutine mf_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
567 integer,
intent(in) :: ixi^
l, ixo^
l
568 double precision,
intent(in) :: qdt, dtfactor
569 double precision,
intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:
ndim)
570 double precision,
intent(inout) :: w(ixi^s,1:nw)
571 logical,
intent(in) :: qsourcesplit
572 logical,
intent(inout) :: active
574 if (.not. qsourcesplit)
then
577 if (abs(
mf_eta)>smalldouble)
then
579 call add_source_res2(qdt,ixi^
l,ixo^
l,wct,w,x)
584 call add_source_hyperres(qdt,ixi^
l,ixo^
l,wct,w,x)
592 select case (type_divb)
597 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
600 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
603 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
606 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
607 case (divb_lindejanhunen)
609 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
610 call add_source_janhunen(qdt,ixi^
l,ixo^
l,wct,w,x)
611 case (divb_lindepowel)
613 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
614 call add_source_powel(qdt,ixi^
l,ixo^
l,wct,w,x)
617 call add_source_linde(qdt,ixi^
l,ixo^
l,wct,w,x)
618 call add_source_glm(qdt,ixi^
l,ixo^
l,wct,w,x)
621 case (divb_multigrid)
624 call mpistop(
'Unknown divB fix')
629 end subroutine mf_add_source
631 subroutine frictional_velocity(w,x,ixI^L,ixO^L)
634 integer,
intent(in) :: ixi^
l, ixo^
l
635 double precision,
intent(in) :: x(ixi^s,1:
ndim)
636 double precision,
intent(inout) :: w(ixi^s,1:nw)
638 double precision :: xmin(
ndim),xmax(
ndim)
639 double precision :: decay(ixo^s)
640 double precision :: current(ixi^s,7-2*
ndir:3),tmp(ixo^s)
641 integer :: ix^
d, idirmin,idir,jdir,kdir
647 if(
block%is_physical_boundary(2*^
d-1))
then
650 if(2*^
d-1==5.and.
ndim==3)
then
652 current(ixomin^
d^
d%ixO^s,:)=current(ixomin^
d+1^
d%ixO^s,:)
657 current(ixomin^
d^
d%ixO^s,:)=current(ixomin^
d+1^
d%ixO^s,:)
661 if(
block%is_physical_boundary(2*^
d))
then
662 current(ixomax^
d^
d%ixO^s,:)=current(ixomax^
d-1^
d%ixO^s,:)
667 do idir=1,
ndir;
do jdir=idirmin,3;
do kdir=1,
ndir
668 if(
lvc(idir,jdir,kdir)/=0)
then
669 if(
lvc(idir,jdir,kdir)==1)
then
670 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))+current(ixo^s,jdir)*w(ixo^s,mag(kdir))
672 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))-current(ixo^s,jdir)*w(ixo^s,mag(kdir))
675 end do;
end do;
end do
677 tmp(ixo^s)=sum(w(ixo^s,mag(:))**2,dim=ndim+1)
679 where(tmp(ixo^s)/=0.d0)
680 tmp(ixo^s)=1.d0/(tmp(ixo^s)*
mf_nu)
684 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))*tmp(ixo^s)
688 ^d&xmin(^d)=xprobmin^d\
689 ^d&xmax(^d)=xprobmax^d\
699 tmp(ixo^s)=sqrt(sum(w(ixo^s,
mom(:))**2,dim=ndim+1))/
mf_vmax+1.d-12
700 tmp(ixo^s)=dtanh(tmp(ixo^s))/tmp(ixo^s)
702 w(ixo^s,
mom(idir))=w(ixo^s,
mom(idir))*tmp(ixo^s)*decay(ixo^s)
705 end subroutine frictional_velocity
711 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
716 integer,
intent(in) :: ixi^
l, ixo^
l
717 double precision,
intent(in) :: qdt
718 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
719 double precision,
intent(inout) :: w(ixi^s,1:nw)
722 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
723 double precision :: gradeta(ixi^s,1:
ndim), bf(ixi^s,1:
ndir)
724 double precision :: tmp(ixi^s),tmp2(ixi^s)
725 integer :: ixa^
l,idir,jdir,kdir,idirmin,idim,jxo^
l,hxo^
l,ix
726 integer :: lxo^
l, kxo^
l
735 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
736 call mpistop(
"Error in add_source_res1: Non-conforming input limits")
743 gradeta(ixo^s,1:
ndim)=zero
749 gradeta(ixo^s,idim)=tmp(ixo^s)
753 bf(ixi^s,1:
ndir)=wct(ixi^s,mag(1:
ndir))
759 tmp2(ixi^s)=bf(ixi^s,idir)
761 lxo^
l=ixo^
l+2*
kr(idim,^
d);
762 jxo^
l=ixo^
l+
kr(idim,^
d);
763 hxo^
l=ixo^
l-
kr(idim,^
d);
764 kxo^
l=ixo^
l-2*
kr(idim,^
d);
765 tmp(ixo^s)=tmp(ixo^s)+&
766 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
771 tmp2(ixi^s)=bf(ixi^s,idir)
773 jxo^
l=ixo^
l+
kr(idim,^
d);
774 hxo^
l=ixo^
l-
kr(idim,^
d);
775 tmp(ixo^s)=tmp(ixo^s)+&
776 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/
dxlevel(idim)**2
781 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
785 do jdir=1,
ndim;
do kdir=idirmin,3
786 if (
lvc(idir,jdir,kdir)/=0)
then
787 if (
lvc(idir,jdir,kdir)==1)
then
788 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
790 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
797 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
801 end subroutine add_source_res1
805 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
810 integer,
intent(in) :: ixi^
l, ixo^
l
811 double precision,
intent(in) :: qdt
812 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
813 double precision,
intent(inout) :: w(ixi^s,1:nw)
816 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
817 double precision :: tmpvec(ixi^s,1:3)
818 integer :: ixa^
l,idir,idirmin,idirmin1
825 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
826 call mpistop(
"Error in add_source_res2: Non-conforming input limits")
840 tmpvec(ixa^s,1:
ndir)=zero
842 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
849 w(ixo^s,mag(
ndir)) = w(ixo^s,mag(
ndir))-qdt*curlj(ixo^s,
ndir)
852 w(ixo^s,mag(1:
ndir)) = w(ixo^s,mag(1:
ndir))-qdt*curlj(ixo^s,1:
ndir)
855 end subroutine add_source_res2
859 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
863 integer,
intent(in) :: ixi^
l, ixo^
l
864 double precision,
intent(in) :: qdt
865 double precision,
intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
866 double precision,
intent(inout) :: w(ixi^s,1:nw)
868 double precision :: current(ixi^s,7-2*
ndir:3)
869 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
870 integer :: ixa^
l,idir,jdir,kdir,idirmin,idirmin1
873 if (iximin^
d>ixamin^
d.or.iximax^
d<ixamax^
d|.or.) &
874 call mpistop(
"Error in add_source_hyperres: Non-conforming input limits")
877 tmpvec(ixa^s,1:
ndir)=zero
879 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
883 call curlvector(tmpvec,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
886 tmpvec(ixa^s,1:
ndir)=zero
887 call curlvector(tmpvec2,ixi^
l,ixa^
l,tmpvec,idirmin1,1,3)
891 tmpvec2(ixa^s,1:
ndir)=zero
892 call curlvector(ehyper,ixi^
l,ixa^
l,tmpvec2,idirmin1,1,3)
895 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
898 end subroutine add_source_hyperres
900 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
907 integer,
intent(in) :: ixi^
l, ixo^
l
908 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
909 double precision,
intent(inout) :: w(ixi^s,1:nw)
910 double precision:: divb(ixi^s)
911 double precision :: gradpsi(ixi^s)
939 end subroutine add_source_glm
942 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
945 integer,
intent(in) :: ixi^
l, ixo^
l
946 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
947 double precision,
intent(inout) :: w(ixi^s,1:nw)
949 double precision :: divb(ixi^s)
957 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,
mom(idir))*divb(ixo^s)
960 end subroutine add_source_powel
962 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
967 integer,
intent(in) :: ixi^
l, ixo^
l
968 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
969 double precision,
intent(inout) :: w(ixi^s,1:nw)
970 double precision :: divb(ixi^s)
978 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,
mom(idir))*divb(ixo^s)
981 end subroutine add_source_janhunen
983 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
988 integer,
intent(in) :: ixi^
l, ixo^
l
989 double precision,
intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:
ndim)
990 double precision,
intent(inout) :: w(ixi^s,1:nw)
991 double precision :: divb(ixi^s),graddivb(ixi^s)
992 integer :: idim, idir, ixp^
l, i^
d, iside
993 logical,
dimension(-1:1^D&) :: leveljump
1001 if(i^
d==0|.and.) cycle
1002 if(neighbor_type(i^
d,
block%igrid)==2 .or. neighbor_type(i^
d,
block%igrid)==4)
then
1003 leveljump(i^
d)=.true.
1005 leveljump(i^
d)=.false.
1014 i^dd=kr(^dd,^d)*(2*iside-3);
1015 if (leveljump(i^dd))
then
1017 ixpmin^d=ixomin^d-i^d
1019 ixpmax^d=ixomax^d-i^d
1030 select case(typegrad)
1032 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
1034 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
1038 if (slab_uniform)
then
1039 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
1041 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
1042 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
1045 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
1048 end subroutine add_source_linde
1056 integer,
intent(in) :: ixi^
l, ixo^
l
1057 double precision,
intent(in) :: w(ixi^s,1:nw)
1058 double precision :: divb(ixi^s), dsurface(ixi^s)
1060 double precision :: invb(ixo^s)
1061 integer :: ixa^
l,idims
1063 call get_divb(w,ixi^
l,ixo^
l,divb)
1064 invb(ixo^s)=sqrt(sum(w(ixo^s, mag(:))**2, dim=
ndim+1))
1065 where(invb(ixo^s)/=0.d0)
1066 invb(ixo^s)=1.d0/invb(ixo^s)
1069 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/
dxlevel(:))
1071 ixamin^
d=ixomin^
d-1;
1072 ixamax^
d=ixomax^
d-1;
1073 dsurface(ixo^s)= sum(
block%surfaceC(ixo^s,:),dim=
ndim+1)
1075 ixa^
l=ixo^
l-
kr(idims,^
d);
1076 dsurface(ixo^s)=dsurface(ixo^s)+
block%surfaceC(ixa^s,idims)
1078 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
1079 block%dvolume(ixo^s)/dsurface(ixo^s)
1090 integer,
intent(in) :: ixo^
l, ixi^
l
1091 double precision,
intent(in) :: w(ixi^s,1:nw)
1092 integer,
intent(out) :: idirmin
1093 logical,
intent(in),
optional :: fourthorder
1096 double precision :: current(ixi^s,7-2*
ndir:3)
1097 integer :: idir, idirmin0
1101 if(
present(fourthorder))
then
1110 subroutine mf_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
1114 integer,
intent(in) :: ixi^
l, ixo^
l
1115 double precision,
intent(inout) :: dtnew
1116 double precision,
intent(in) ::
dx^
d
1117 double precision,
intent(in) :: w(ixi^s,1:nw)
1118 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1120 double precision :: dxarr(
ndim)
1121 double precision :: current(ixi^s,7-2*
ndir:3),eta(ixi^s)
1122 integer :: idirmin,idim
1129 else if (
mf_eta<zero)
then
1136 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
1139 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/
block%ds(ixo^s,idim)**2)))
1151 end subroutine mf_get_dt
1154 subroutine mf_add_source_geom(qdt,dtfactor, ixI^L,ixO^L,wCT,wprim,w,x)
1158 integer,
intent(in) :: ixi^
l, ixo^
l
1159 double precision,
intent(in) :: qdt, dtfactor, x(ixi^s,1:
ndim)
1160 double precision,
intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
1162 double precision :: tmp,invr,cs2
1163 integer :: iw,idir,ix^
d
1165 integer :: mr_,mphi_
1166 integer :: br_,bphi_
1169 br_=mag(1); bphi_=mag(1)-1+
phi_
1175 {
do ix^db=ixomin^db,ixomax^db\}
1176 w(ix^
d,bphi_)=w(ix^
d,bphi_)+qdt/x(ix^
d,1)*&
1177 (wct(ix^
d,bphi_)*wct(ix^
d,mr_) &
1178 -wct(ix^
d,br_)*wct(ix^
d,mphi_))
1182 if(
mf_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,
psi_)/x(ixo^s,1)
1184 {
do ix^db=ixomin^db,ixomax^db\}
1188 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wct(ix^d,
psi_)
1193 if(.not.stagger_grid)
then
1194 tmp=wct(ix^d,
mom(1))*wct(ix^d,mag(2))-wct(ix^d,
mom(2))*wct(ix^d,mag(1))
1195 cs2=1.d0/tan(x(ix^d,2))
1197 tmp=tmp+cs2*wct(ix^d,
psi_)
1199 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
1206 if(.not.stagger_grid)
then
1207 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
1208 (wct(ix^d,
mom(1))*wct(ix^d,mag(3)) &
1209 -wct(ix^d,
mom(3))*wct(ix^d,mag(1)))
1214 if(.not.stagger_grid)
then
1215 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
1216 (wct(ix^d,
mom(1))*wct(ix^d,mag(3)) &
1217 -wct(ix^d,
mom(3))*wct(ix^d,mag(1)) &
1218 -(wct(ix^d,
mom(3))*wct(ix^d,mag(2)) &
1219 -wct(ix^d,
mom(2))*wct(ix^d,mag(3)))*cs2)
1226 end subroutine mf_add_source_geom
1228 subroutine mf_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
1231 integer,
intent(in) :: ixi^
l, ixo^
l, idir
1232 double precision,
intent(in) :: qt
1233 double precision,
intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
1234 double precision,
intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
1236 double precision :: db(ixi^s), dpsi(ixi^s)
1244 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
1245 dpsi(ixo^s) = wrp(ixo^s,
psi_) - wlp(ixo^s,
psi_)
1247 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
1249 wlp(ixo^s,
psi_) = 0.5d0 * (wrp(ixo^s,
psi_) + wlp(ixo^s,
psi_)) &
1252 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
1255 end subroutine mf_modify_wlr
1257 subroutine mf_boundary_adjust(igrid,psb)
1259 integer,
intent(in) :: igrid
1262 integer :: ib, idims, iside, ixo^
l, i^
d
1271 i^
d=
kr(^
d,idims)*(2*iside-3);
1272 if (neighbor_type(i^
d,igrid)/=1) cycle
1273 ib=(idims-1)*2+iside
1291 call fixdivb_boundary(ixg^
ll,ixo^
l,psb(igrid)%w,psb(igrid)%x,ib)
1296 end subroutine mf_boundary_adjust
1298 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
1301 integer,
intent(in) :: ixg^
l,ixo^
l,ib
1302 double precision,
intent(inout) :: w(ixg^s,1:nw)
1303 double precision,
intent(in) :: x(ixg^s,1:
ndim)
1305 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
1306 integer :: ix^
d,ixf^
l
1318 do ix1=ixfmax1,ixfmin1,-1
1319 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
1320 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
1321 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
1324 do ix1=ixfmax1,ixfmin1,-1
1325 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
1326 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
1327 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
1328 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
1329 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
1330 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
1331 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
1345 do ix1=ixfmax1,ixfmin1,-1
1346 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1347 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
1348 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
1349 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
1350 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
1351 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
1354 do ix1=ixfmax1,ixfmin1,-1
1355 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1356 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
1357 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
1358 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
1359 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
1360 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
1361 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
1362 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
1363 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
1364 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
1365 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
1366 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
1367 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
1368 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
1369 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1370 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
1371 /
block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
1372 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
1384 do ix1=ixfmin1,ixfmax1
1385 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
1386 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
1387 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
1390 do ix1=ixfmin1,ixfmax1
1391 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
1392 w(ix1,ixfmin2:ixfmax2,mag(1)))*
block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
1393 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
1394 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
1395 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
1396 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
1397 /
block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
1411 do ix1=ixfmin1,ixfmax1
1412 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1413 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
1414 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
1415 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
1416 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
1417 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
1420 do ix1=ixfmin1,ixfmax1
1421 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1422 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
1423 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
1424 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
1425 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
1426 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
1427 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
1428 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
1429 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
1430 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
1431 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
1432 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
1433 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
1434 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
1435 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1436 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
1437 /
block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
1438 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
1450 do ix2=ixfmax2,ixfmin2,-1
1451 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
1452 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
1453 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
1456 do ix2=ixfmax2,ixfmin2,-1
1457 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
1458 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
1459 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
1460 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
1461 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
1462 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
1463 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
1477 do ix2=ixfmax2,ixfmin2,-1
1478 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
1479 ix2+1,ixfmin3:ixfmax3,mag(2)) &
1480 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
1481 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
1482 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
1483 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
1486 do ix2=ixfmax2,ixfmin2,-1
1487 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
1488 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
1489 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
1490 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
1491 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
1492 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1493 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
1494 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
1495 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1496 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
1497 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
1498 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
1499 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
1500 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
1501 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1502 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
1503 /
block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
1504 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
1516 do ix2=ixfmin2,ixfmax2
1517 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
1518 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
1519 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
1522 do ix2=ixfmin2,ixfmax2
1523 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
1524 w(ixfmin1:ixfmax1,ix2,mag(2)))*
block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
1525 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
1526 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
1527 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
1528 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
1529 /
block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
1543 do ix2=ixfmin2,ixfmax2
1544 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
1545 ix2-1,ixfmin3:ixfmax3,mag(2)) &
1546 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
1547 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
1548 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
1549 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
1552 do ix2=ixfmin2,ixfmax2
1553 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
1554 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
1555 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
1556 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
1557 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
1558 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1559 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
1560 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
1561 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1562 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
1563 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
1564 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
1565 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
1566 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
1567 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1568 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
1569 /
block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
1570 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
1585 do ix3=ixfmax3,ixfmin3,-1
1586 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
1587 ixfmin2:ixfmax2,ix3+1,mag(3)) &
1588 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
1589 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
1590 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
1591 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
1594 do ix3=ixfmax3,ixfmin3,-1
1595 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
1596 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
1597 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
1598 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
1599 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
1600 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1601 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
1602 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
1603 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1604 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
1605 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
1606 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
1607 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
1608 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
1609 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
1610 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
1611 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
1612 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
1625 do ix3=ixfmin3,ixfmax3
1626 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
1627 ixfmin2:ixfmax2,ix3-1,mag(3)) &
1628 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
1629 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
1630 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
1631 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
1634 do ix3=ixfmin3,ixfmax3
1635 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
1636 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
1637 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
1638 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
1639 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
1640 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1641 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
1642 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
1643 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1644 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
1645 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
1646 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
1647 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
1648 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
1649 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
1650 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
1651 /
block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
1652 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
1657 call mpistop(
"Special boundary is not defined for this region")
1660 end subroutine fixdivb_boundary
1669 double precision,
intent(in) :: qdt
1670 double precision,
intent(in) :: qt
1671 logical,
intent(inout) :: active
1673 double precision :: tmp(ixg^t), grad(ixg^t,
ndim)
1674 double precision :: res
1675 double precision,
parameter :: max_residual = 1
d-3
1676 double precision,
parameter :: residual_reduction = 1
d-10
1678 integer,
parameter :: max_its = 50
1679 double precision :: residual_it(max_its), max_divb
1680 integer :: iigrid, igrid
1681 integer :: n, nc, lvl, ix^
l, ixc^
l, idim
1684 mg%operator_type = mg_laplacian
1692 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1693 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1696 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
1697 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1699 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1700 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1703 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
1704 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1708 write(*,*)
"mf_clean_divb_multigrid warning: unknown boundary type"
1709 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1710 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1718 do iigrid = 1, igridstail
1719 igrid = igrids(iigrid);
1722 lvl =
mg%boxes(id)%lvl
1723 nc =
mg%box_size_lvl(lvl)
1729 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^
ll,
ixm^
ll, tmp, &
1731 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(
ixm^t)
1732 max_divb = max(max_divb, maxval(abs(tmp(
ixm^t))))
1737 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
1740 if (
mype == 0) print *,
"Performing multigrid divB cleaning"
1741 if (
mype == 0) print *,
"iteration vs residual"
1744 call mg_fas_fmg(
mg, n>1, max_res=residual_it(n))
1745 if (
mype == 0)
write(*,
"(I4,E11.3)") n, residual_it(n)
1746 if (residual_it(n) < residual_reduction * max_divb)
exit
1748 if (
mype == 0 .and. n > max_its)
then
1749 print *,
"divb_multigrid warning: not fully converged"
1750 print *,
"current amplitude of divb: ", residual_it(max_its)
1751 print *,
"multigrid smallest grid: ", &
1752 mg%domain_size_lvl(:,
mg%lowest_lvl)
1753 print *,
"note: smallest grid ideally has <= 8 cells"
1754 print *,
"multigrid dx/dy/dz ratio: ",
mg%dr(:, 1)/
mg%dr(1, 1)
1755 print *,
"note: dx/dy/dz should be similar"
1759 call mg_fas_vcycle(
mg, max_res=res)
1760 if (res < max_residual)
exit
1762 if (res > max_residual)
call mpistop(
"divb_multigrid: no convergence")
1767 do iigrid = 1, igridstail
1768 igrid = igrids(iigrid);
1777 tmp(ix^s) =
mg%boxes(id)%cc({:,}, mg_iphi)
1781 ixcmin^
d=ixmlo^
d-
kr(idim,^
d);
1783 call gradientf(tmp,ps(igrid)%x,ixg^
ll,ixc^
l,idim,grad(ixg^t,idim))
1784 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
1792 ps(igrid)%w(
ixm^t, mag(1:
ndim)) = &
1793 ps(igrid)%w(
ixm^t, mag(1:
ndim)) - grad(
ixm^t, :)
1803 subroutine mf_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
1806 integer,
intent(in) :: ixi^
l, ixo^
l
1807 double precision,
intent(in) :: qt, qdt
1809 double precision,
intent(in) :: wprim(ixi^s,1:nw)
1810 type(state) :: sct, s
1811 type(ct_velocity) :: vcts
1812 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
1813 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
1817 call update_faces_average(ixi^
l,ixo^
l,qt,qdt,fc,fe,sct,s)
1819 call update_faces_contact(ixi^
l,ixo^
l,qt,qdt,wprim,fc,fe,sct,s,vcts)
1821 call update_faces_hll(ixi^
l,ixo^
l,qt,qdt,fe,sct,s,vcts)
1823 call mpistop(
'choose average, uct_contact,or uct_hll for type_ct!')
1826 end subroutine mf_update_faces
1829 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
1833 integer,
intent(in) :: ixi^
l, ixo^
l
1834 double precision,
intent(in) :: qt, qdt
1835 type(state) :: sct, s
1836 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
1837 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
1839 double precision :: circ(ixi^s,1:
ndim)
1840 double precision :: e_resi(ixi^s,
sdim:3)
1841 integer :: ix^
d,ixc^
l,ixa^
l,i1kr^
d,i2kr^
d
1842 integer :: idim1,idim2,idir,iwdim1,iwdim2
1844 associate(bfaces=>s%ws,x=>s%x)
1851 if(
mf_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,e_resi)
1855 i1kr^
d=
kr(idim1,^
d);
1858 i2kr^
d=
kr(idim2,^
d);
1861 if (
lvc(idim1,idim2,idir)==1)
then
1863 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
1865 {
do ix^db=ixcmin^db,ixcmax^db\}
1866 fe(ix^
d,idir)=quarter*&
1867 (fc(ix^
d,iwdim1,idim2)+fc({ix^
d+i1kr^
d},iwdim1,idim2)&
1868 -fc(ix^
d,iwdim2,idim1)-fc({ix^
d+i2kr^
d},iwdim2,idim1))
1870 if(
mf_eta/=zero) fe(ix^
d,idir)=fe(ix^
d,idir)+e_resi(ix^
d,idir)
1873 fe(ix^
d,idir)=fe(ix^
d,idir)*qdt*s%dsC(ix^
d,idir)
1881 if(
associated(usr_set_electric_field)) &
1882 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
1884 circ(ixi^s,1:ndim)=zero
1890 ixcmin^d=ixomin^d-kr(idim1,^d);
1892 ixa^l=ixc^l-kr(idim2,^d);
1895 if(lvc(idim1,idim2,idir)==1)
then
1897 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
1900 else if(lvc(idim1,idim2,idir)==-1)
then
1902 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
1909 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
1911 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
1916 end subroutine update_faces_average
1919 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
1923 integer,
intent(in) :: ixi^
l, ixo^
l
1924 double precision,
intent(in) :: qt, qdt
1926 double precision,
intent(in) :: wp(ixi^s,1:nw)
1927 type(state) :: sct, s
1928 type(ct_velocity) :: vcts
1929 double precision,
intent(in) :: fc(ixi^s,1:nwflux,1:
ndim)
1930 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
1932 double precision :: circ(ixi^s,1:
ndim)
1934 double precision :: ecc(ixi^s,
sdim:3)
1936 double precision :: el(ixi^s),er(ixi^s)
1938 double precision :: elc(ixi^s),erc(ixi^s)
1940 double precision :: jce(ixi^s,
sdim:3)
1941 integer :: hxc^
l,ixc^
l,jxc^
l,ixa^
l,ixb^
l
1942 integer :: idim1,idim2,idir,iwdim1,iwdim2
1944 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
1949 if(
lvc(idim1,idim2,idir)==1)
then
1950 ecc(ixi^s,idir)=ecc(ixi^s,idir)+wp(ixi^s,mag(idim1))*wp(ixi^s,
mom(idim2))
1951 else if(
lvc(idim1,idim2,idir)==-1)
then
1952 ecc(ixi^s,idir)=ecc(ixi^s,idir)-wp(ixi^s,mag(idim1))*wp(ixi^s,
mom(idim2))
1957 if(
mf_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,jce)
1969 if (
lvc(idim1,idim2,idir)==1)
then
1971 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
1973 jxc^
l=ixc^
l+
kr(idim1,^
d);
1974 hxc^
l=ixc^
l+
kr(idim2,^
d);
1976 fe(ixc^s,idir)=quarter*&
1977 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
1978 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
1982 ixamax^
d=ixcmax^
d+
kr(idim1,^
d);
1983 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
1984 hxc^
l=ixa^
l+
kr(idim2,^
d);
1985 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
1986 where(vnorm(ixc^s,idim1)>0.d0)
1987 elc(ixc^s)=el(ixc^s)
1988 else where(vnorm(ixc^s,idim1)<0.d0)
1989 elc(ixc^s)=el(jxc^s)
1991 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
1993 hxc^
l=ixc^
l+
kr(idim2,^
d);
1994 where(vnorm(hxc^s,idim1)>0.d0)
1995 erc(ixc^s)=er(ixc^s)
1996 else where(vnorm(hxc^s,idim1)<0.d0)
1997 erc(ixc^s)=er(jxc^s)
1999 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
2001 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
2004 jxc^
l=ixc^
l+
kr(idim2,^
d);
2006 ixamax^
d=ixcmax^
d+
kr(idim2,^
d);
2007 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
2008 hxc^
l=ixa^
l+
kr(idim1,^
d);
2009 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
2010 where(vnorm(ixc^s,idim2)>0.d0)
2011 elc(ixc^s)=el(ixc^s)
2012 else where(vnorm(ixc^s,idim2)<0.d0)
2013 elc(ixc^s)=el(jxc^s)
2015 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
2017 hxc^
l=ixc^
l+
kr(idim1,^
d);
2018 where(vnorm(hxc^s,idim2)>0.d0)
2019 erc(ixc^s)=er(ixc^s)
2020 else where(vnorm(hxc^s,idim2)<0.d0)
2021 erc(ixc^s)=er(jxc^s)
2023 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
2025 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
2028 if(
mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2032 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
2047 circ(ixi^s,1:
ndim)=zero
2052 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
2056 if(
lvc(idim1,idim2,idir)/=0)
then
2057 hxc^
l=ixc^
l-
kr(idim2,^
d);
2059 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2060 +
lvc(idim1,idim2,idir)&
2067 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
2068 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2070 circ(ixc^s,idim1)=zero
2073 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2078 end subroutine update_faces_contact
2081 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
2086 integer,
intent(in) :: ixi^
l, ixo^
l
2087 double precision,
intent(in) :: qt, qdt
2088 double precision,
intent(inout) :: fe(ixi^s,
sdim:3)
2089 type(state) :: sct, s
2090 type(ct_velocity) :: vcts
2092 double precision :: vtill(ixi^s,2)
2093 double precision :: vtilr(ixi^s,2)
2094 double precision :: btill(ixi^s,
ndim)
2095 double precision :: btilr(ixi^s,
ndim)
2096 double precision :: cp(ixi^s,2)
2097 double precision :: cm(ixi^s,2)
2098 double precision :: circ(ixi^s,1:
ndim)
2100 double precision :: jce(ixi^s,
sdim:3)
2101 integer :: hxc^
l,ixc^
l,ixcp^
l,jxc^
l,ixcm^
l
2102 integer :: idim1,idim2,idir
2104 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
2105 cbarmax=>vcts%cbarmax)
2118 if(
mf_eta/=zero)
call get_resistive_electric_field(ixi^
l,ixo^
l,sct,s,jce)
2131 ixcmin^
d=ixomin^
d-1+
kr(idir,^
d);
2135 idim2=mod(idir+1,3)+1
2137 jxc^
l=ixc^
l+
kr(idim1,^
d);
2138 ixcp^
l=ixc^
l+
kr(idim2,^
d);
2142 vtill(ixi^s,2),vtilr(ixi^s,2))
2145 vtill(ixi^s,1),vtilr(ixi^s,1))
2151 btill(ixi^s,idim1),btilr(ixi^s,idim1))
2154 btill(ixi^s,idim2),btilr(ixi^s,idim2))
2158 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
2159 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
2161 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
2162 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
2166 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
2167 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
2168 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
2169 /(cp(ixc^s,1)+cm(ixc^s,1)) &
2170 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
2171 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
2172 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
2173 /(cp(ixc^s,2)+cm(ixc^s,2))
2176 if(
mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2180 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
2194 circ(ixi^s,1:
ndim)=zero
2200 ixcmin^
d=ixomin^
d-
kr(idim1,^
d);
2204 if(
lvc(idim1,idim2,idir)/=0)
then
2205 hxc^
l=ixc^
l-
kr(idim2,^
d);
2207 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2208 +
lvc(idim1,idim2,idir)&
2215 where(s%surfaceC(ixc^s,idim1) > 1.0
d-9*s%dvolume(ixc^s))
2216 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2218 circ(ixc^s,idim1)=zero
2221 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2225 end subroutine update_faces_hll
2228 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
2233 integer,
intent(in) :: ixi^
l, ixo^
l
2234 type(state),
intent(in) :: sct, s
2236 double precision :: jce(ixi^s,
sdim:3)
2239 double precision :: jcc(ixi^s,7-2*
ndir:3)
2241 double precision :: xs(ixgs^t,1:
ndim)
2243 double precision :: eta(ixi^s)
2244 double precision :: gradi(ixgs^t)
2245 integer :: ix^
d,ixc^
l,ixa^
l,ixb^
l,idir,idirmin,idim1,idim2
2247 associate(x=>s%x,
dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
2253 if (
lvc(idim1,idim2,idir)==0) cycle
2254 ixcmax^
d=ixomax^
d+1;
2255 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-2;
2256 ixbmax^
d=ixcmax^
d-
kr(idir,^
d)+1;
2259 xs(ixb^s,:)=x(ixb^s,:)
2260 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*
dx(ixb^s,idim2)
2261 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^
ll,ixc^
l,idim1,gradi)
2262 if (
lvc(idim1,idim2,idir)==1)
then
2263 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
2265 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
2272 jce(ixi^s,:)=jce(ixi^s,:)*
mf_eta
2280 ixcmin^
d=ixomin^
d+
kr(idir,^
d)-1;
2281 jcc(ixc^s,idir)=0.d0
2283 if({ ix^
d==1 .and. ^
d==idir | .or.}) cycle
2284 ixamin^
d=ixcmin^
d+ix^
d;
2285 ixamax^
d=ixcmax^
d+ix^
d;
2286 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
2288 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
2289 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
2294 end subroutine get_resistive_electric_field
2300 integer,
intent(in) :: ixo^
l
2307 do ix^db=ixomin^db,ixomax^db\}
2309 s%w(ix^
d,mag(1))=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
2310 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
2311 s%w(ix^
d,mag(2))=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
2312 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
2313 s%w(ix^
d,mag(3))=half/s%surface(ix^
d,3)*(s%ws(ix^
d,3)*s%surfaceC(ix^
d,3)&
2314 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
2317 s%w(ix^
d,mag(1))=half/s%surface(ix^
d,1)*(s%ws(ix^
d,1)*s%surfaceC(ix^
d,1)&
2318 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
2319 s%w(ix^
d,mag(2))=half/s%surface(ix^
d,2)*(s%ws(ix^
d,2)*s%surfaceC(ix^
d,2)&
2320 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
2363 integer,
intent(in) :: ixis^
l, ixi^
l, ixo^
l
2364 double precision,
intent(inout) :: ws(ixis^s,1:nws)
2365 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2367 double precision :: adummy(ixis^s,1:3)
2376 double precision :: sum_jbb_ipe, sum_j_ipe, sum_l_ipe, f_i_ipe, volume_pe
2377 double precision :: sum_jbb, sum_j, sum_l, f_i, volume, cw_sin_theta
2378 integer :: iigrid, igrid, ix^
d
2379 integer :: amode, istatus(mpi_status_size)
2380 integer,
save :: fhmf
2381 logical :: patchwi(ixg^t)
2382 logical,
save :: logmfopened=.false.
2383 character(len=800) :: filename,filehead
2384 character(len=800) :: line,datastr
2391 do iigrid=1,igridstail; igrid=igrids(iigrid);
2394 call mask_inner(ixg^
ll,
ixm^
ll,ps(igrid)%w,ps(igrid)%x,patchwi,volume_pe)
2395 sum_jbb_ipe = sum_jbb_ipe+integral_grid_mf(ixg^
ll,
ixm^
ll,ps(igrid)%w,&
2396 ps(igrid)%x,1,patchwi)
2397 sum_j_ipe = sum_j_ipe+integral_grid_mf(ixg^
ll,
ixm^
ll,ps(igrid)%w,&
2398 ps(igrid)%x,2,patchwi)
2399 f_i_ipe=f_i_ipe+integral_grid_mf(ixg^
ll,
ixm^
ll,ps(igrid)%w,&
2400 ps(igrid)%x,3,patchwi)
2401 sum_l_ipe = sum_l_ipe+integral_grid_mf(ixg^
ll,
ixm^
ll,ps(igrid)%w,&
2402 ps(igrid)%x,4,patchwi)
2404 call mpi_reduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
2406 call mpi_reduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,0,&
2408 call mpi_reduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,0,&
2410 call mpi_reduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,0,&
2412 call mpi_reduce(volume_pe,volume,1,mpi_double_precision,mpi_sum,0,&
2418 cw_sin_theta = sum_jbb/sum_j
2424 if(.not.logmfopened)
then
2426 write(filename,
"(a,a)") trim(
base_filename),
"_mf_metrics.csv"
2430 open(unit=20,file=trim(filename),status=
'replace')
2431 close(20, status=
'delete')
2434 amode=ior(mpi_mode_create,mpi_mode_wronly)
2435 amode=ior(amode,mpi_mode_append)
2436 call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,
ierrmpi)
2438 filehead=
" it,time,CW_sin_theta,current,Lorenz_force,f_i"
2440 call mpi_file_write(fhmf,filehead,len_trim(filehead), &
2441 mpi_character,istatus,
ierrmpi)
2442 call mpi_file_write(fhmf,achar(10),1,mpi_character,istatus,
ierrmpi)
2446 write(datastr,
'(i6,a)')
it,
','
2447 line=trim(line)//trim(datastr)
2449 line=trim(line)//trim(datastr)
2450 write(datastr,
'(es13.6,a)') cw_sin_theta,
','
2451 line=trim(line)//trim(datastr)
2452 write(datastr,
'(es13.6,a)') sum_j,
','
2453 line=trim(line)//trim(datastr)
2454 write(datastr,
'(es13.6,a)') sum_l,
','
2455 line=trim(line)//trim(datastr)
2456 write(datastr,
'(es13.6)') f_i
2457 line=trim(line)//trim(datastr)//new_line(
'A')
2458 call mpi_file_write(fhmf,line,len_trim(line),mpi_character,istatus,
ierrmpi)
2460 call mpi_file_close(fhmf,
ierrmpi)
2467 subroutine mask_inner(ixI^L,ixO^L,w,x,patchwi,volume_pe)
2470 integer,
intent(in) :: ixi^
l,ixo^
l
2471 double precision,
intent(in):: w(ixi^s,nw),x(ixi^s,1:
ndim)
2472 double precision,
intent(inout) :: volume_pe
2473 logical,
intent(out) :: patchwi(ixi^s)
2475 double precision :: xo^
l
2478 {xomin^
d = xprobmin^
d + 0.05d0*(xprobmax^
d-xprobmin^
d)\}
2479 {xomax^
d = xprobmax^
d - 0.05d0*(xprobmax^
d-xprobmin^
d)\}
2481 xomin^nd = xprobmin^nd
2486 {
do ix^db=ixomin^db,ixomax^db\}
2487 if({ x(ix^dd,^
d) > xomin^
d .and. x(ix^dd,^
d) < xomax^
d | .and. })
then
2488 patchwi(ix^
d)=.true.
2489 volume_pe=volume_pe+
block%dvolume(ix^
d)
2491 patchwi(ix^
d)=.false.
2495 end subroutine mask_inner
2497 function integral_grid_mf(ixI^L,ixO^L,w,x,iw,patchwi)
2500 integer,
intent(in) :: ixi^
l,ixo^
l,iw
2501 double precision,
intent(in) :: x(ixi^s,1:
ndim)
2502 double precision :: w(ixi^s,nw+
nwauxio)
2503 logical,
intent(in) :: patchwi(ixi^s)
2505 double precision,
dimension(ixI^S,7-2*ndir:3) :: current
2506 double precision,
dimension(ixI^S,1:ndir) :: bvec,qvec
2507 double precision :: integral_grid_mf,tmp(ixi^s),bm2
2508 integer :: ix^
d,idirmin0,idirmin,idir,jdir,kdir
2510 integral_grid_mf=0.d0
2514 bvec(ixo^s,:)=w(ixo^s,mag(:))
2517 qvec(ixo^s,1:
ndir)=zero
2518 do idir=1,
ndir;
do jdir=idirmin,3;
do kdir=1,
ndir
2519 if(
lvc(idir,jdir,kdir)/=0)
then
2520 if(
lvc(idir,jdir,kdir)==1)
then
2521 qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
2523 qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
2526 end do;
end do;
end do
2528 {
do ix^db=ixomin^db,ixomax^db\}
2529 if(patchwi(ix^d))
then
2530 bm2=sum(bvec(ix^d,:)**2)
2531 if(bm2/=0.d0) bm2=1.d0/bm2
2532 integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2)*&
2533 bm2)*block%dvolume(ix^d)
2539 {
do ix^db=ixomin^db,ixomax^db\}
2540 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(current(ix^d,:)**2))*&
2547 {
do ix^db=ixomin^db,ixomax^db\}
2548 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+tmp(ix^d)*block%dvolume(ix^d)
2552 bvec(ixo^s,:)=w(ixo^s,mag(:))
2555 qvec(ixo^s,1:ndir)=zero
2556 do idir=1,ndir;
do jdir=idirmin,3;
do kdir=1,ndir
2557 if(lvc(idir,jdir,kdir)/=0)
then
2558 if(lvc(idir,jdir,kdir)==1)
then
2559 qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
2561 qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
2564 end do;
end do;
end do
2566 {
do ix^db=ixomin^db,ixomax^db\}
2567 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+&
2568 sqrt(sum(qvec(ix^d,:)**2))*block%dvolume(ix^d)
2572 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)
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...
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.
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.