38 double precision ::
tmf 45 integer,
private,
protected :: rho_
48 integer,
allocatable,
private,
protected :: mom(:)
51 integer,
allocatable,
private,
protected :: mag(:)
62 character(len=*),
intent(in) :: files(:)
68 open(
unitpar, file=trim(files(n)), status=
"old")
107 double precision :: dvolume(ixg^t),dsurface(ixg^t),dvone
108 double precision :: dtfff,dtfff_pe,dtnew,dx^D
109 double precision :: cwsin_theta_new,cwsin_theta_old
110 double precision :: sum_jbb,sum_jbb_ipe,sum_j,sum_j_ipe,sum_l_ipe,sum_l
111 double precision :: f_i_ipe,f_i,volumepe,volume,tmpt,time_in
112 double precision,
external :: integral_grid
113 integer :: i,iigrid, igrid, idims,ix^D,hxM^LL,fhmf,tmpit,i^D
114 logical :: patchwi(ixg^t), stagger_flag
121 if(
mype==0)
write(*,*)
'Evolving to force-free field using magnetofricitonal method...' 153 do iigrid=1,igridstail; igrid=igrids(iigrid);
170 do iigrid=1,igridstail; igrid=igrids(iigrid);
172 pso(igrid)%w(ixg^t,mag(:))=ps(igrid)%w(ixg^t,mag(:))
175 dtfff_pe=min(dtfff_pe,dtnew)
177 call mpi_allreduce(dtfff_pe,dtfff,1,mpi_double_precision,mpi_min, &
198 if(mod(i,10)==0)
then 206 do iigrid=1,igridstail; igrid=igrids(iigrid);
212 do iigrid=1,igridstail; igrid=igrids(iigrid);
218 write(*,*)
'<CW sin theta>:',cwsin_theta_new
219 write(*,*)
'<f_i>:',f_i
220 write(*,*)
'----------------------------------------------------------' 226 if(mod(i,10)/=0)
then 232 write (*,*)
'Reach maximum iteration step!' 233 write (*,*)
'The total iteration step is:', i
247 do iigrid=1,igridstail; igrid=igrids(iigrid);
248 ps(igrid)%w(ixg^t,mom(:))=zero
257 if(
mype==0)
write(*,*)
'Magnetofriction phase took : ',mpi_wtime()-time_in,
' sec' 267 do iigrid=1,igridstail; igrid=igrids(iigrid);
269 dvolume(
ixm^t)=block%dvolume(
ixm^t)
271 hxm^ll=
ixm^ll-
kr(idims,^d);
272 dsurface(
ixm^t)=dsurface(
ixm^t)+block%surfaceC(hxm^t,idims)
277 ps(igrid)%x,1,patchwi)
279 ps(igrid)%x,2,patchwi)
281 ps(igrid)%x,3,patchwi)
283 ps(igrid)%x,4,patchwi)
285 call mpi_allreduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
287 call mpi_allreduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,&
289 call mpi_allreduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,&
291 call mpi_allreduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,&
293 call mpi_allreduce(volumepe,volume,1,mpi_double_precision,mpi_sum,&
297 cwsin_theta_new = sum_jbb/sum_j
308 integer,
intent(in) :: ixI^L,ixO^L
309 double precision,
intent(in):: w(ixi^s,nw),x(ixi^s,1:
ndim)
310 double precision :: xO^L
313 {xomin^d = xprobmin^d + 0.05d0*(xprobmax^d-xprobmin^d)\}
314 {xomax^d = xprobmax^d - 0.05d0*(xprobmax^d-xprobmin^d)\}
316 xomin^nd = xprobmin^nd
321 {
do ix^db=ixomin^db,ixomax^db\}
322 if({ x(ix^dd,^d) > xomin^d .and. x(ix^dd,^d) < xomax^d | .and. })
then 324 volumepe=volumepe+dvolume(ix^d)
326 patchwi(ix^d)=.false.
333 integer :: amode, status(mpi_status_size)
334 character(len=800) :: filename,filehead
335 character(len=2048) :: line,datastr
336 logical,
save :: logmfopened=.false.
339 if(.not.logmfopened)
then 343 amode=ior(mpi_mode_create,mpi_mode_wronly)
344 amode=ior(amode,mpi_mode_append)
345 call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,
ierrmpi)
347 filehead=
" itmf, dt, <f_i>, <CW sin theta>, <Current>, <Lorenz force>" 348 call mpi_file_write(fhmf,filehead,len_trim(filehead), &
350 call mpi_file_write(fhmf,achar(10),1,mpi_character,status,
ierrmpi)
353 write(datastr,
'(i6,a)') i,
',' 354 line=trim(line)//trim(datastr)
355 write(datastr,
'(es13.6,a)') dtfff,
',' 356 line=trim(line)//trim(datastr)
357 write(datastr,
'(es13.6,a)') f_i,
',' 358 line=trim(line)//trim(datastr)
359 write(datastr,
'(es13.6,a)') cwsin_theta_new,
',' 360 line=trim(line)//trim(datastr)
361 write(datastr,
'(es13.6,a)') sum_j,
',' 362 line=trim(line)//trim(datastr)
363 write(datastr,
'(es13.6)') sum_l
364 line=trim(line)//trim(datastr)//new_line(
'A')
365 call mpi_file_write(fhmf,line,len_trim(line),mpi_character,status,
ierrmpi)
373 integer,
intent(in) :: ixI^L,ixO^L,iw
374 double precision,
intent(in) :: x(ixi^s,1:
ndim)
375 double precision :: w(ixi^s,nw+
nwauxio)
376 logical,
intent(in) :: patchwi(ixi^s)
378 double precision,
dimension(ixI^S,1:ndir) :: bvec,qvec,current
379 double precision :: integral_grid_mf,tmp(ixi^s),b_mag(ixi^s)
380 integer :: ix^D,i,idirmin,idir,jdir,kdir
382 integral_grid_mf=0.d0
387 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
389 bvec(ixi^s,:)=w(ixi^s,mag(:))
393 qvec(ixo^s,1:
ndir)=zero
394 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
395 if(
lvc(idir,jdir,kdir)/=0)
then 396 tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
397 if(
lvc(idir,jdir,kdir)==1)
then 398 qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
400 qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
405 {
do ix^db=ixomin^db,ixomax^db\}
407 sum(bvec(ix^d,:)**2))*dvolume(ix^d)
412 {
do ix^db=ixomin^db,ixomax^db\}
420 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
422 bvec(ixi^s,:)=w(ixi^s,mag(:))
424 call divvector(bvec,ixi^
l,ixo^
l,tmp)
425 {
do ix^db=ixomin^db,ixomax^db\}
427 dvolume(ix^d)**2/sqrt(sum(bvec(ix^d,:)**2))/dsurface(ix^d)
432 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
434 bvec(ixi^s,:)=w(ixi^s,mag(:))
436 call curlvector(bvec,ixi^
l,ixo^
l,current,idirmin,1,
ndir)
438 qvec(ixo^s,1:
ndir)=zero
439 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
440 if(
lvc(idir,jdir,kdir)/=0)
then 441 tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
442 if(
lvc(idir,jdir,kdir)==1)
then 443 qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
445 qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
450 {
do ix^db=ixomin^db,ixomax^db\}
462 double precision,
intent(in) :: dtfff
463 integer :: i,iigrid, igrid
464 double precision :: vhatmax,vhatmax_pe,vhatmaxgrid
466 vhatmax_pe=smalldouble
467 do iigrid=1,igridstail; igrid=igrids(iigrid);
470 call vhat(ps(igrid)%w,ps(igrid)%x,ixg^
ll,
ixm^
ll,vhatmaxgrid)
471 vhatmax_pe=max(vhatmax_pe,vhatmaxgrid)
473 call mpi_allreduce(vhatmax_pe,vhatmax,1,mpi_double_precision,mpi_max, &
475 do iigrid=1,igridstail; igrid=igrids(iigrid);
484 subroutine vhat(w,x,ixI^L,ixO^L,vhatmaxgrid)
488 integer,
intent(in) :: ixI^L, ixO^L
489 double precision,
intent(inout) :: w(ixi^s,nw)
490 double precision,
intent(in) :: x(ixi^s,1:
ndim)
491 double precision,
intent(out) :: vhatmaxgrid
493 double precision :: current(ixi^s,7-2*
ndir:3),tmp(ixi^s),dxhm
494 double precision :: dxhms(ixo^s)
495 integer :: idirmin,idir,jdir,kdir
500 do idir=1,
ndir;
do jdir=1,
ndir;
do kdir=idirmin,3
501 if(
lvc(idir,jdir,kdir)/=0)
then 503 tmp(ixo^s)=current(ixo^s,jdir)*(w(ixo^s,mag(kdir))+block%b0(ixo^s,kdir,0))
505 tmp(ixo^s)=current(ixo^s,jdir)*w(ixo^s,mag(kdir))
507 if(
lvc(idir,jdir,kdir)==1)
then 508 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+tmp(ixo^s)
510 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-tmp(ixo^s)
516 tmp(ixo^s)=sum((w(ixo^s,mag(:))+block%b0(ixo^s,:,0))**2,dim=ndim+1)
518 tmp(ixo^s)=sum(w(ixo^s,mag(:))**2,dim=ndim+1)
521 if(slab_uniform)
then 522 dxhm=dble(ndim)/(^d&1.0d0/dxlevel(^d)+)
524 w(ixo^s,mom(idir))=dxhm*w(ixo^s,mom(idir))/tmp(ixo^s)
527 dxhms(ixo^s)=dble(ndim)/sum(1.d0/block%dx(ixo^s,:),dim=ndim+1)
529 w(ixo^s,mom(idir))=dxhms(ixo^s)*w(ixo^s,mom(idir))/tmp(ixo^s)
532 vhatmaxgrid=maxval(sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1)))
539 integer,
intent(in) :: ixI^L, ixO^L
540 double precision,
intent(in) :: x(ixi^s,1:
ndim),qdt,qvmax
541 double precision,
intent(inout) :: w(ixi^s,1:nw)
543 double precision :: dxhm,disbd(6),bfzone^D
544 double precision :: dxhms(ixo^s)
545 integer :: ix^D, idir
552 w(ixo^s,mom(:))=w(ixo^s,mom(:))*dxhm
554 dxhms(ixo^s)=dble(
ndim)/sum(1.d0/block%dx(ixo^s,:),dim=
ndim+1)
555 dxhms(ixo^s)=
mf_cc*
mf_cy/qvmax*dxhms(ixo^s)/qdt
558 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*dxhms(ixo^s)
563 bfzone1=0.05d0*(xprobmax1-xprobmin1)
564 bfzone2=0.05d0*(xprobmax2-xprobmin2)
565 bfzone3=0.05d0*(xprobmax3-xprobmin3)
566 {
do ix^db=ixomin^db,ixomax^db\}
567 disbd(1)=x(ix^d,1)-xprobmin1
568 disbd(2)=xprobmax1-x(ix^d,1)
569 disbd(3)=x(ix^d,2)-xprobmin2
570 disbd(4)=xprobmax2-x(ix^d,2)
571 disbd(5)=x(ix^d,3)-xprobmin1
572 disbd(6)=xprobmax3-x(ix^d,3)
575 if(disbd(1)<bfzone1)
then 576 w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(1))/bfzone1)**2)*w(ix^d,mom(:))
579 if(disbd(5)<bfzone3)
then 580 w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(5))/bfzone3)**2)*w(ix^d,mom(:))
583 if(disbd(2)<bfzone1)
then 584 w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(2))/bfzone1)**2)*w(ix^d,mom(:))
586 if(disbd(3)<bfzone2)
then 587 w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(3))/bfzone2)**2)*w(ix^d,mom(:))
589 if(disbd(4)<bfzone2)
then 590 w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(4))/bfzone2)**2)*w(ix^d,mom(:))
592 if(disbd(6)<bfzone3)
then 593 w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(6))/bfzone3)**2)*w(ix^d,mom(:))
599 subroutine advectmf(idim^LIM,qt,qdt)
606 integer,
intent(in) :: idim^LIM
607 double precision,
intent(in) :: qt, qdt
609 integer :: iigrid, igrid
615 do iigrid=1,igridstail; igrid=igrids(iigrid);
616 ps1(igrid)%w=ps(igrid)%w
635 do iigrid=1,igridstail; igrid=igrids(iigrid);
636 ps2(igrid)%w(ixg^t,1:nwflux)=0.75d0*ps(igrid)%w(ixg^t,1:nwflux)+0.25d0*&
637 ps1(igrid)%w(ixg^t,1:nwflux)
638 if (nw>nwflux) ps2(igrid)%w(ixg^t,nwflux+1:nw) = &
639 ps(igrid)%w(ixg^t,nwflux+1:nw)
644 do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
645 ps(igrid)%w(ixg^t,1:nwflux)=1.0d0/3.0d0*ps(igrid)%w(ixg^t,1:nwflux)+&
646 2.0d0/3.0d0*ps2(igrid)%w(ixg^t,1:nwflux)
652 write(
unitterm,*)
"Error in advectmf: Unknown time integration method" 653 call mpistop(
"Correct time_integrator")
658 subroutine advect1mf(method,dtin,dtfactor,idim^LIM,qtC,psa,qt,psb)
666 integer,
intent(in) :: idim^LIM
669 double precision,
intent(in) :: dtin,dtfactor, qtC, qt
670 character(len=*),
intent(in) :: method(
nlevelshi)
672 double precision :: qdt
673 integer :: iigrid, igrid, level, i^D
680 do iigrid=1,igridstail; igrid=igrids(iigrid);
685 psa(igrid)%w,qt,psb(igrid)%w,pso(igrid)%w)
719 subroutine process1_gridmf(method,igrid,qdt,ixG^L,idim^LIM,qtC,wCT,qt,w,wold)
724 character(len=*),
intent(in) :: method
725 integer,
intent(in) :: igrid, ixG^L, idim^LIM
726 double precision,
intent(in) :: qdt, qtC, qt
727 double precision :: wCT(ixg^s,1:nw), w(ixg^s,1:nw), wold(ixg^s,1:nw)
728 double precision :: dx^D, fC(ixg^s,1:
ndir,1:
ndim)
739 ixo^l=ixg^l^lsubnghostcells;
745 call centdiff4mf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
750 call tvdlfmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
753 call hancockmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,dx^d,ps(igrid)%x)
758 call fdmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
760 if(
mype==0)
write(
unitterm,*)
'Error in advect1_gridmf:',method,
' is not there!' 761 call mpistop(
"The scheme is not implemented.")
770 subroutine upwindlrmf(ixI^L,ixL^L,ixR^L,idim,w,wCT,wLC,wRC,x)
776 integer,
intent(in) :: ixI^L, ixL^L, ixR^L, idim
777 double precision,
dimension(ixI^S,1:nw) :: w, wCT
778 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC
779 double precision,
dimension(ixI^S,1:ndim) :: x
781 integer :: jxR^L, ixC^L, jxC^L, iw
782 double precision :: ldw(ixi^s), rdw(ixi^s), dwC(ixi^s)
785 call mp5limiter(ixi^l,ixl^l,idim,w,wlc,wrc)
787 call ppmlimiter(ixi^l,
ixm^
ll,idim,w,wct,wlc,wrc)
789 jxr^l=ixr^l+
kr(idim,^
d);
790 ixcmax^
d=jxrmax^
d; ixcmin^
d=ixlmin^
d-
kr(idim,^
d);
791 jxc^l=ixc^l+
kr(idim,^
d);
795 w(ixcmin^
d:jxcmax^
d,iw)=dlog10(w(ixcmin^
d:jxcmax^
d,iw))
796 wlc(ixl^s,iw)=dlog10(wlc(ixl^s,iw))
797 wrc(ixr^s,iw)=dlog10(wrc(ixr^s,iw))
800 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
804 wlc(ixl^s,iw)=wlc(ixl^s,iw)+half*ldw(ixl^s)
805 wrc(ixr^s,iw)=wrc(ixr^s,iw)-half*rdw(jxr^s)
808 w(ixcmin^
d:jxcmax^
d,iw)=10.0d0**w(ixcmin^
d:jxcmax^
d,iw)
809 wlc(ixl^s,iw)=10.0d0**wlc(ixl^s,iw)
810 wrc(ixr^s,iw)=10.0d0**wrc(ixr^s,iw)
818 subroutine getfluxmf(w,x,ixI^L,ixO^L,idir,idim,f)
822 integer,
intent(in) :: ixI^L, ixO^L, idir, idim
823 double precision,
intent(in) :: w(ixi^s,nw)
824 double precision,
intent(in) :: x(ixi^s,1:
ndim)
825 double precision,
intent(out) :: f(ixi^s)
832 f(ixo^s)=w(ixo^s,mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom(idir))
835 +w(ixo^s,mom(idim))*block%B0(ixo^s,idir,idim)&
836 -w(ixo^s,mom(idir))*block%B0(ixo^s,idim,idim)
842 subroutine tvdlfmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,wold,fC,dx^D,x)
847 double precision,
intent(in) :: qdt, qtC, qt, dx^D
848 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
849 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
850 double precision,
dimension(ixI^S,1:nw) :: wCT, wnew, wold
851 double precision,
dimension(ixI^S,1:ndir,1:ndim) :: fC
853 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC, wmean
854 double precision,
dimension(ixI^S) :: fLC, fRC
855 double precision,
dimension(ixI^S) :: cmaxC
856 double precision :: dxinv(1:
ndim), inv_volume(ixo^s)
857 integer :: idims, idir, ix^L, hxO^L, ixC^L, ixCR^L, jxC^L, kxC^L, kxR^L
863 ix^l=ix^l^ladd2*
kr(idims,^d);
865 if (ixi^l^ltix^l|.or.|.or.) &
866 call mpistop(
"Error in tvdlfmf: Nonconforming input limits")
868 ^d&dxinv(^d)=-qdt/dx^d;
873 hxo^l=ixo^l-
kr(idims,^d);
875 ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
877 jxc^l=ixc^l+
kr(idims,^d);
878 kxcmin^d=iximin^d; kxcmax^d=iximax^d-
kr(idims,^d);
879 kxr^l=kxc^l+
kr(idims,^d);
882 wrc(kxc^s,1:nwflux)=wct(kxr^s,1:nwflux)
883 wlc(kxc^s,1:nwflux)=wct(kxc^s,1:nwflux)
885 call upwindlrmf(ixi^l,ixcr^l,ixcr^l,idims,wct,wct,wlc,wrc,x)
890 wmean=0.5d0*(wlc+wrc)
898 flc(ixc^s)=half*(flc(ixc^s)+frc(ixc^s))
901 if (idir==idims)
then 902 flc(ixc^s)=flc(ixc^s)-
mf_tvdlfeps*
tvdlfeps*cmaxc(ixc^s)*half*(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
905 fc(ixc^s,idir,idims)=flc(ixc^s)
907 fc(ixc^s,idir,idims)=block%surfaceC(ixc^s,idims)*flc(ixc^s)
915 hxo^l=ixo^l-
kr(idims,^d);
918 fc(ixi^s,:,idims)=dxinv(idims)*fc(ixi^s,:,idims)
919 wnew(ixo^s,mag(:))=wnew(ixo^s,mag(:)) &
920 + (fc(ixo^s,:,idims)-fc(hxo^s,:,idims))
922 inv_volume = 1.0d0/block%dvolume(ixo^s)
923 fc(ixi^s,:,idims)=-qdt*fc(ixi^s,:,idims)
926 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir)) + (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims)) * &
934 call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
938 subroutine hancockmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,dx^D,x)
946 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
947 double precision,
intent(in) :: qdt, qtC, qt, dx^D, x(ixi^s,1:
ndim)
948 double precision,
intent(inout) :: wCT(ixi^s,1:nw), wnew(ixi^s,1:nw)
950 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC
951 double precision,
dimension(ixI^S) :: fLC, fRC
952 double precision :: dxinv(1:
ndim)
953 integer :: idims, idir, ix^L, hxO^L, ixtest^L
958 ix^l=ix^l^laddkr(idims,^d);
960 if (ixi^l^ltix^l|.or.|.or.) &
961 call mpistop(
"Error in Hancockmf: Nonconforming input limits")
963 ^d&dxinv(^d)=-qdt/dx^d;
969 hxo^l=ixo^l-
kr(idims,^d);
971 wrc(hxo^s,1:nwflux)=wct(ixo^s,1:nwflux)
972 wlc(ixo^s,1:nwflux)=wct(ixo^s,1:nwflux)
974 call upwindlrmf(ixi^l,ixo^l,hxo^l,idims,wct,wct,wlc,wrc,x)
979 call getfluxmf(wrc,x,ixi^l,hxo^l,idir,idims,frc)
980 call getfluxmf(wlc,x,ixi^l,ixo^l,idir,idims,flc)
983 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+dxinv(idims)* &
984 (flc(ixo^s)-frc(hxo^s))
986 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))-qdt/block%dvolume(ixo^s) &
987 *(block%surfaceC(ixo^s,idims)*flc(ixo^s) &
988 -block%surfaceC(hxo^s,idims)*frc(hxo^s))
997 subroutine fdmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,wold,fC,dx^D,x)
999 double precision,
intent(in) :: qdt, qtC, qt, dx^D
1000 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
1001 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
1003 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: wCT, wnew, wold
1004 double precision,
dimension(ixI^S,1:ndir,1:ndim),
intent(out) :: fC
1006 double precision,
dimension(ixI^S) :: fCT
1007 double precision,
dimension(ixI^S,1:nw) :: fm, fp, fmR, fpL
1008 double precision,
dimension(ixI^S) :: v
1009 double precision :: dxinv(1:
ndim)
1010 integer :: idims, idir, ixC^L, ix^L, hxO^L, ixCR^L
1012 ^d&dxinv(^d)=-qdt/dx^d;
1019 hxo^l=ixo^l-
kr(idims,^d);
1021 ixmax^d=ixomax^d; ixmin^d=hxomin^d;
1025 call getfluxmf(wct,x,ixg^
ll,ixcr^l,idir,idims,fct)
1037 fc(ix^s,idir,idims) = dxinv(idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
1038 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
1039 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1041 fc(ix^s,idir,idims)=-qdt*block%surfaceC(ix^s,idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
1042 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
1043 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/block%dvolume(ixo^s)
1050 call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
1058 integer,
intent(in) :: ixI^L, iL^L, idims
1059 double precision,
intent(in) :: w(ixi^s,1:nw)
1061 double precision,
intent(out) :: wLC(ixi^s,1:nw)
1063 double precision :: ldw(ixi^s), dwC(ixi^s)
1064 integer :: jxR^L, ixC^L, jxC^L, kxC^L, iw
1068 call mp5limiterl(ixi^l,il^l,idims,w,wlc)
1071 kxcmin^
d=iximin^
d; kxcmax^
d=iximax^
d-
kr(idims,^
d);
1073 wlc(kxc^s,1:nwflux) = w(kxc^s,1:nwflux)
1075 jxr^l=il^l+
kr(idims,^
d);
1077 ixcmax^
d=jxrmax^
d; ixcmin^
d=ilmin^
d-
kr(idims,^
d);
1078 jxc^l=ixc^l+
kr(idims,^
d);
1081 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1085 wlc(il^s,iw)=wlc(il^s,iw)+half*ldw(il^s)
1095 integer,
intent(in) :: ixI^L, iL^L, idims
1096 double precision,
intent(in) :: w(ixi^s,1:nw)
1098 double precision,
intent(out) :: wRC(ixi^s,1:nw)
1100 double precision :: rdw(ixi^s), dwC(ixi^s)
1101 integer :: jxR^L, ixC^L, jxC^L, kxC^L, kxR^L, iw
1105 call mp5limiterr(ixi^l,il^l,idims,w,wrc)
1108 kxcmin^
d=iximin^
d; kxcmax^
d=iximax^
d-
kr(idims,^
d);
1109 kxr^l=kxc^l+
kr(idims,^
d);
1111 wrc(kxc^s,1:nwflux)=w(kxr^s,1:nwflux)
1113 jxr^l=il^l+
kr(idims,^
d);
1114 ixcmax^
d=jxrmax^
d; ixcmin^
d=ilmin^
d-
kr(idims,^
d);
1115 jxc^l=ixc^l+
kr(idims,^
d);
1118 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1121 wrc(il^s,iw)=wrc(il^s,iw)-half*rdw(jxr^s)
1127 subroutine centdiff4mf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,w,wold,fC,dx^D,x)
1135 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
1136 double precision,
intent(in) :: qdt, qtC, qt, dx^D
1137 double precision :: wCT(ixi^s,1:nw), w(ixi^s,1:nw), wold(ixi^s,1:nw)
1138 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1139 double precision :: fC(ixi^s,1:
ndir,1:
ndim)
1141 double precision :: v(ixi^s,
ndim), f(ixi^s)
1142 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC
1143 double precision,
dimension(ixI^S) :: vLC, vRC,cmaxLC,cmaxRC
1144 double precision :: dxinv(1:
ndim)
1145 integer :: idims, idir, idirmin,ix^D
1146 integer :: ix^L, hxO^L, ixC^L, jxC^L, hxC^L, kxC^L, kkxC^L, kkxR^L
1151 ix^l=ix^l^ladd2*
kr(idims,^d);
1154 if (ixi^l^ltix^l|.or.|.or.)
then 1155 call mpistop(
"Error in evolve_CentDiff4: Non-conforming input limits")
1157 ^d&dxinv(^d)=-qdt/dx^d;
1162 ix^l=ixo^l^ladd2*
kr(idims,^d);
1163 hxo^l=ixo^l-
kr(idims,^d);
1165 ixcmin^d=hxomin^d; ixcmax^d=ixomax^d;
1166 hxc^l=ixc^l-
kr(idims,^d);
1167 jxc^l=ixc^l+
kr(idims,^d);
1168 kxc^l=ixc^l+2*
kr(idims,^d);
1170 kkxcmin^d=iximin^d; kkxcmax^d=iximax^d-
kr(idims,^d);
1171 kkxr^l=kkxc^l+
kr(idims,^d);
1172 wrc(kkxc^s,1:nwflux)=wct(kkxr^s,1:nwflux)
1173 wlc(kkxc^s,1:nwflux)=wct(kkxc^s,1:nwflux)
1175 call upwindlrmf(ixi^l,ixc^l,ixc^l,idims,wct,wct,wlc,wrc,x)
1181 vlc(ixc^s)=max(cmaxrc(ixc^s),cmaxlc(ixc^s))
1185 call getfluxmf(wct,x,ixi^l,ix^l,idir,idims,f)
1188 fc(ixc^s,idir,idims)=(-f(kxc^s)+7.0d0*(f(jxc^s)+f(ixc^s))-f(hxc^s))/12.0d0
1193 *(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
1196 fc(ixc^s,idir,idims)=dxinv(idims)*fc(ixc^s,idir,idims)
1198 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+(fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1200 fc(ixc^s,idir,idims)=-qdt*block%surfaceC(ixc^s,idims)*fc(ixc^s,idir,idims)
1201 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+ &
1202 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/block%dvolume(ixo^s)
1216 integer,
intent(in) :: ixI^L, ixO^L
1217 double precision,
intent(in) :: x(ixi^s,1:
ndim)
1218 double precision,
intent(inout) :: w(ixi^s,1:nw), dtnew
1220 double precision :: courantmax, dxinv(1:
ndim)
1221 double precision :: cmax(ixi^s),tmp(ixi^s),alfven(ixi^s)
1232 tmp(ixo^s)=cmax(ixo^s)/block%dx(ixo^s,idims)
1233 courantmax=max(courantmax,maxval(tmp(ixo^s)))
1235 tmp(ixo^s)=cmax(ixo^s)*dxinv(idims)
1236 courantmax=max(courantmax,maxval(tmp(ixo^s)))
1240 if (courantmax>smalldouble) dtnew=min(dtnew,
mf_cc/courantmax)
1244 subroutine getcmaxfff(w,ixI^L,ixO^L,idims,cmax)
1247 logical :: new_cmax,needcmin
1248 integer,
intent(in) :: ixI^L, ixO^L, idims
1249 double precision,
intent(in) :: w(ixi^s,1:nw)
1250 double precision,
intent(out) :: cmax(ixi^s)
1254 cmax(ixo^s)=sqrt(sum((w(ixo^s,mag(:))+block%b0(ixo^s,:,0))**2,dim=
ndim+1)/w(ixo^s,rho_))
1256 cmax(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2,dim=
ndim+1)/w(ixo^s,rho_))
1258 cmax(ixo^s)=cmax(ixo^s)+abs(w(ixo^s,mom(idims)))
1263 subroutine divbclean(qdt,ixI^L,ixO^L,wCT,w,x)
1267 integer,
intent(in) :: ixI^L, ixO^L
1268 double precision,
intent(in) :: x(ixi^s,1:
ndim),wCT(ixi^s,1:nw),qdt
1269 double precision,
intent(inout) :: w(ixi^s,1:nw)
1270 integer :: idims, ix^L, ixp^L, i^D, iside
1271 double precision :: divb(ixi^s),graddivb(ixi^s),bdivb(ixi^s,1:
ndir)
1282 call gradient(divb,ixi^l,ixp^l,idims,graddivb)
1288 graddivb(ixp^s)=graddivb(ixp^s)*
mf_cdivb &
1289 /(^d&1.0d0/block%dx(ixp^s,^d)**2+)
1296 w(ixp^s,mag(idims))=w(ixp^s,mag(idims))+&
1307 integer,
intent(in) :: ixI^L, ixO^L
1308 double precision,
intent(in) :: qdt, x(ixi^s,1:
ndim)
1309 double precision,
intent(inout) :: wCT(ixi^s,1:nw), w(ixi^s,1:nw)
1311 double precision :: tmp(ixi^s)
1313 integer :: mr_,mphi_
1314 integer :: br_,bphi_
1316 mr_=mom(1); mphi_=mom(1)-1+
phi_ 1317 br_=mag(1); bphi_=mag(1)-1+
phi_ 1323 tmp(ixo^s)=(wct(ixo^s,bphi_)*wct(ixo^s,mom(1)) &
1324 -wct(ixo^s,br_)*wct(ixo^s,mom(3)))
1325 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt*tmp(ixo^s)/x(ixo^s,1)
1331 tmp(ixo^s)= wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
1332 -wct(ixo^s,mom(2))*wct(ixo^s,mag(1))
1334 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*block%b0(ixo^s,2,0) &
1335 -wct(ixo^s,mom(2))*block%b0(ixo^s,1,0)
1338 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
1343 tmp(ixo^s)=wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
1344 -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)){^nooned &
1345 -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
1346 -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
1349 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*block%b0(ixo^s,3,0) &
1350 -wct(ixo^s,mom(3))*block%b0(ixo^s,1,0){^nooned &
1351 -(wct(ixo^s,mom(3))*block%b0(ixo^s,2,0) &
1352 -wct(ixo^s,mom(2))*block%b0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
1356 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
1364 subroutine get_current(w,ixI^L,ixO^L,idirmin,current)
1369 integer :: ixO^L, idirmin, ixI^L
1370 double precision :: w(ixi^s,1:nw)
1374 double precision :: current(ixi^s,7-2*
ndir:3),bvec(ixi^s,1:
ndir)
1378 bvec(ixi^s,1:
ndir)=w(ixi^s,mag(1:
ndir))
1382 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
1383 block%J0(ixo^s,idirmin0:3)
1388 subroutine get_divb(w,ixI^L,ixO^L,divb)
1392 integer,
intent(in) :: ixI^L, ixO^L
1393 double precision,
intent(in) :: w(ixi^s,1:nw)
1394 double precision :: divb(ixi^s)
1396 double precision :: bvec(ixi^s,1:
ndir)
1398 bvec(ixi^s,:)=w(ixi^s,mag(:))
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_p2
integer, parameter cylindrical
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_f
subroutine fdmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, wold, fC, dxD, x)
This module contains definitions of global parameters and variables and some generic functions/subrou...
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_p1
character(len=std_len) time_integrator
Which time integrator to use.
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_f
subroutine get_current(w, ixIL, ixOL, idirmin, current)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
update ghost cells of all blocks including physical boundaries
subroutine mf_velocity_update(dtfff)
procedure(p_no_args), pointer usr_before_main_loop
subroutine, public init_comm_fix_conserve(idimLIM, nwfluxin)
double precision mf_tvdlfeps_min
subroutine get_divb(w, ixIL, ixOL, divb)
Calculate div B within ixO.
integer, dimension(:^d &,:^d &), pointer type_recv_srl
character(len=std_len), dimension(:), allocatable typepred1
The spatial dicretization for the predictor step when using a two step method.
subroutine create_bc_mpi_datatype(nwstart, nwbc)
integer, dimension(:^d &,:^d &), pointer type_send_r
integer max_blocks
The maximum number of grid blocks in a processor.
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_p2
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
subroutine tvdlfmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, wold, fC, dxD, x)
integer, dimension(:), allocatable type_gradient_limiter
Type of slope limiter used for computing gradients or divergences, when typegrad or typediv are set t...
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges...
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_p2
integer, parameter limiter_mp5
Module with geometry-related routines (e.g., divergence, curl)
subroutine, public store_flux(igrid, fC, idimLIM, nwfluxin)
subroutine reconstructrmf(ixIL, iLL, idims, w, wRC)
integer, parameter plevel_
subroutine centdiff4mf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, w, wold, fC, dxD, x)
integer ndir
Number of spatial dimensions (components) for vector variables.
subroutine divbclean(qdt, ixIL, ixOL, wCT, w, x)
Clean divergence of magnetic field by Janhunen's and Linde's source terms.
integer, dimension(:), allocatable type_limiter
Type of slope limiter used for reconstructing variables on cell edges.
subroutine addgeometrymf(qdt, ixIL, ixOL, wCT, w, x)
character(len=std_len) base_filename
Base file name for simulation output, which will be followed by a number.
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_f
Module for flux conservation near refinement boundaries.
integer istep
Index of the sub-step in a multi-step time integrator.
double precision mf_cy
frictional velocity coefficient
subroutine advectmf(idimLIM, qt, qdt)
subroutine resettree
reset AMR and (de)allocate boundary flux storage at level changes
integer, dimension(:^d &,:^d &), pointer type_recv_r
integer, dimension(:^d &,:^d &), pointer type_recv_p
subroutine getfluxmf(w, x, ixIL, ixOL, idir, idim, f)
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_f
subroutine getcmaxfff(w, ixIL, ixOL, idims, cmax)
logical b0field
split magnetic field as background B0 field
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_p2
character(len=std_len), dimension(:), allocatable flux_scheme
Which spatial discretization to use (per grid level)
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision mf_cy_max
integer, parameter limiter_ppm
double precision mf_cdivb
divb cleaning coefficient controls diffusion speed of divb
subroutine, public sendflux(idimLIM)
double precision function integral_grid_mf(ixIL, ixOL, w, x, iw, patchwi)
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_p1
double precision cmax_mype
maximal speed for fd scheme
Module with slope/flux limiters.
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_p1
logical fix_conserve_at_step
integer nghostcells
Number of ghost cells surrounding a grid.
subroutine process1_gridmf(method, igrid, qdt, ixGL, idimLIM, qtC, wCT, qt, w, wold)
Module with all the methods that users can customize in AMRVAC.
logical, dimension(:), allocatable loglimit
integer ditregrid
Reconstruct the AMR grid once every ditregrid iteration(s)
double precision mf_cc
stability coefficient controls numerical stability
logical stagger_grid
True for using stagger grid.
procedure(sub_convert), pointer phys_to_conserved
integer ierrmpi
A global MPI error return code.
integer ixm
the mesh range (within a block with ghost cells)
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_p2
subroutine reconstructlmf(ixIL, iLL, idims, w, wLC)
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
integer, parameter unitterm
Unit for standard output.
integer, parameter unitpar
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, dimension(:), allocatable, parameter d
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder)
Calculate divergence of a vector qvec within ixL.
integer mype
The rank of the current MPI task.
double precision global_time
The global simulation time.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
double precision mf_tvdlfeps
TVDLF dissipation coefficient controls the dissipation term.
subroutine frictional_velocity(w, x, ixIL, ixOL, qvmax, qdt)
subroutine mf_params_read(files)
Read this module"s parameters from a file.
subroutine, public fix_conserve(psb, idimLIM, nw0, nwfluxin)
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_f
double precision, dimension(ndim) dxlevel
integer snapshotini
Resume from the snapshot with this index.
double precision cmax_global
maximal speed for fd scheme
integer, parameter ndim
Number of spatial dimensions for grid variables.
subroutine getdtfff_courant(w, x, ixIL, ixOL, dtnew)
double precision tmf
time in magnetofriction process
subroutine mask_inner(ixIL, ixOL, w, x)
integer, dimension(:^d &,:^d &), pointer type_send_p
logical slab
Cartesian geometry or not.
subroutine advect1mf(method, dtin, dtfactor, idimLIM, qtC, psa, qt, psb)
subroutine vhat(w, x, ixIL, ixOL, vhatmaxgrid)
integer nwauxio
Number of auxiliary variables that are only included in output.
integer icomm
The MPI communicator.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
character(len=std_len) typediv
subroutine, public recvflux(idimLIM)
integer, parameter nlevelshi
The maximum number of levels in the grid refinement.
integer, dimension(^nd, 0:3) l
subroutine hancockmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, dxD, x)
integer refine_max_level
Maximal number of AMR levels.
This module defines the procedures of a physics module. It contains function pointers for the various...
procedure(sub_convert), pointer phys_to_primitive
integer, dimension(:,:), allocatable node
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_f
integer, parameter spherical
integer, dimension(:^d &,:^d &), pointer type_send_srl
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_p1
subroutine dwlimiter2(dwC, ixIL, ixCL, idims, typelim, ldw, rdw)
Limit the centered dwC differences within ixC for iw in direction idim. The limiter is chosen accordi...
integer it
Number of time steps taken.
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_p2
subroutine magnetofriction_init()
Initialize the module.
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision tvdlfeps
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_p1
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_p1
subroutine magnetofriction
subroutine upwindlrmf(ixIL, ixLL, ixRL, idim, w, wCT, wLC, wRC, x)