38 double precision ::
tmf
44 double precision,
public ::
mf_vmax = 3.d6
47 integer,
private,
protected :: rho_
50 integer,
allocatable,
private,
protected :: mom(:)
53 integer,
allocatable,
private,
protected :: mag(:)
65 character(len=*),
intent(in) :: files(:)
72 open(
unitpar, file=trim(files(n)), status=
"old")
121 double precision :: dvolume(ixG^T),dsurface(ixG^T),dvone
122 double precision :: dtfff,dtfff_pe,dtnew,dx^D
123 double precision :: cwsin_theta_new,cwsin_theta_old
124 double precision :: sum_jbb,sum_jbb_ipe,sum_j,sum_j_ipe,sum_l_ipe,sum_l
125 double precision :: f_i_ipe,f_i,volumepe,volume,tmpt,time_in
126 double precision,
external :: integral_grid
127 integer :: i,iigrid, igrid, idims,ix^D,hxM^LL,fhmf,tmpit,i^D
128 logical :: patchwi(ixG^T), stagger_flag
135 if(
mype==0)
write(*,*)
'Evolving to force-free field using magnetofricitonal method...'
167 do iigrid=1,igridstail; igrid=igrids(iigrid);
184 do iigrid=1,igridstail; igrid=igrids(iigrid);
188 dtfff_pe=min(dtfff_pe,dtnew)
190 call mpi_allreduce(dtfff_pe,dtfff,1,mpi_double_precision,mpi_min, &
215 if(mod(i,10)==0)
then
223 do iigrid=1,igridstail; igrid=igrids(iigrid);
229 do iigrid=1,igridstail; igrid=igrids(iigrid);
235 write(*,*)
'<CW sin theta>:',cwsin_theta_new
236 write(*,*)
'<f_i>:',f_i
237 write(*,*)
'----------------------------------------------------------'
243 if(mod(i,10)/=0)
then
249 write (*,*)
'Reach maximum iteration step!'
250 write (*,*)
'The total iteration step is:', i
264 do iigrid=1,igridstail; igrid=igrids(iigrid);
265 ps(igrid)%w(ixg^t,mom(:))=zero
274 if(
mype==0)
write(*,*)
'Magnetofriction phase took : ',mpi_wtime()-time_in,
' sec'
285 do iigrid=1,igridstail; igrid=igrids(iigrid);
289 hxm^ll=
ixm^ll-
kr(idims,^d);
290 dsurface(
ixm^t)=dsurface(
ixm^t)+
block%surfaceC(hxm^t,idims)
295 ps(igrid)%x,1,patchwi)
297 ps(igrid)%x,2,patchwi)
299 ps(igrid)%x,3,patchwi)
301 ps(igrid)%x,4,patchwi)
303 call mpi_allreduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
305 call mpi_allreduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,&
307 call mpi_allreduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,&
309 call mpi_allreduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,&
311 call mpi_allreduce(volumepe,volume,1,mpi_double_precision,mpi_sum,&
315 cwsin_theta_new = sum_jbb/sum_j
325 integer,
intent(in) :: ixI^L,ixO^L
326 double precision,
intent(in):: w(ixI^S,nw),x(ixI^S,1:ndim)
327 double precision :: xO^L
330 {xomin^d = xprobmin^d + 0.05d0*(xprobmax^d-xprobmin^d)\}
331 {xomax^d = xprobmax^d - 0.05d0*(xprobmax^d-xprobmin^d)\}
333 xomin^nd = xprobmin^nd
338 {
do ix^db=ixomin^db,ixomax^db\}
339 if({ x(ix^dd,^d) > xomin^d .and. x(ix^dd,^d) < xomax^d | .and. })
then
341 volumepe=volumepe+dvolume(ix^d)
343 patchwi(ix^d)=.false.
350 integer :: amode, status(MPI_STATUS_SIZE)
351 character(len=800) :: filename,filehead
352 character(len=2048) :: line,datastr
353 logical,
save :: logmfopened=.false.
356 if(.not.logmfopened)
then
358 write(filename,
"(a,a)") trim(base_filename),
"_mflog.csv"
360 amode=ior(mpi_mode_create,mpi_mode_wronly)
361 amode=ior(amode,mpi_mode_append)
362 call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,ierrmpi)
364 filehead=
" itmf, dt, <f_i>, <CW sin theta>, <Current>, <Lorenz force>"
365 call mpi_file_write(fhmf,filehead,len_trim(filehead), &
366 mpi_character,status,ierrmpi)
367 call mpi_file_write(fhmf,achar(10),1,mpi_character,status,ierrmpi)
370 write(datastr,
'(i6,a)') i,
','
371 line=trim(line)//trim(datastr)
372 write(datastr,
'(es13.6,a)') dtfff,
','
373 line=trim(line)//trim(datastr)
374 write(datastr,
'(es13.6,a)') f_i,
','
375 line=trim(line)//trim(datastr)
376 write(datastr,
'(es13.6,a)') cwsin_theta_new,
','
377 line=trim(line)//trim(datastr)
378 write(datastr,
'(es13.6,a)') sum_j,
','
379 line=trim(line)//trim(datastr)
380 write(datastr,
'(es13.6)') sum_l
381 line=trim(line)//trim(datastr)//new_line(
'A')
382 call mpi_file_write(fhmf,line,len_trim(line),mpi_character,status,ierrmpi)
390 integer,
intent(in) :: ixi^l,ixo^l,iw
391 double precision,
intent(in) :: x(ixi^s,1:ndim)
392 double precision,
intent(in) :: w(ixi^s,nw+nwauxio)
393 logical,
intent(in) :: patchwi(ixi^s)
395 double precision,
dimension(ixI^S,1:ndir) :: bvec,qvec,current
397 integer :: ix^d,i,idirmin,idir,jdir,kdir
404 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
406 bvec(ixi^s,:)=w(ixi^s,mag(:))
410 qvec(ixo^s,1:ndir)=zero
411 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
412 if(lvc(idir,jdir,kdir)/=0)
then
413 tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
414 if(lvc(idir,jdir,kdir)==1)
then
415 qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
417 qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
422 {
do ix^db=ixomin^db,ixomax^db\}
424 sum(bvec(ix^d,:)**2))*dvolume(ix^d)
429 {
do ix^db=ixomin^db,ixomax^db\}
437 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
439 bvec(ixi^s,:)=w(ixi^s,mag(:))
441 call divvector(bvec,ixi^l,ixo^l,tmp)
442 {
do ix^db=ixomin^db,ixomax^db\}
444 dvolume(ix^d)**2/sqrt(sum(bvec(ix^d,:)**2))/dsurface(ix^d)
449 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
451 bvec(ixi^s,:)=w(ixi^s,mag(:))
453 call curlvector(bvec,ixi^l,ixo^l,current,idirmin,1,ndir)
455 qvec(ixo^s,1:ndir)=zero
456 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
457 if(lvc(idir,jdir,kdir)/=0)
then
458 tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
459 if(lvc(idir,jdir,kdir)==1)
then
460 qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
462 qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
467 {
do ix^db=ixomin^db,ixomax^db\}
479 double precision,
intent(in) :: dtfff
480 double precision :: vhatmax,vhatmax_pe,vhatmaxgrid
481 integer :: i,iigrid, igrid
483 vhatmax_pe=smalldouble
484 do iigrid=1,igridstail; igrid=igrids(iigrid);
487 call vhat(ps(igrid)%w,ps(igrid)%x,ixg^
ll,
ixm^
ll,vhatmaxgrid)
488 vhatmax_pe=max(vhatmax_pe,vhatmaxgrid)
490 call mpi_allreduce(vhatmax_pe,vhatmax,1,mpi_double_precision,mpi_max, &
492 do iigrid=1,igridstail; igrid=igrids(iigrid);
501 subroutine vhat(w,x,ixI^L,ixO^L,vhatmaxgrid)
505 integer,
intent(in) :: ixI^L, ixO^L
506 double precision,
intent(inout) :: w(ixI^S,nw)
507 double precision,
intent(in) :: x(ixI^S,1:ndim)
508 double precision,
intent(out) :: vhatmaxgrid
510 double precision :: current(ixI^S,7-2*ndir:3),tmp(ixI^S),dxhm
511 double precision :: dxhms(ixO^S)
512 integer :: idirmin,idir,jdir,kdir
517 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
518 if(
lvc(idir,jdir,kdir)/=0)
then
520 tmp(ixo^s)=current(ixo^s,jdir)*(w(ixo^s,mag(kdir))+
block%b0(ixo^s,kdir,0))
522 tmp(ixo^s)=current(ixo^s,jdir)*w(ixo^s,mag(kdir))
524 if(
lvc(idir,jdir,kdir)==1)
then
525 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+tmp(ixo^s)
527 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-tmp(ixo^s)
534 tmp(ixo^s)=1.d0/(sum((w(ixo^s,mag(:))+
block%b0(ixo^s,:,0))**2,dim=ndim+1)+smalldouble)
536 tmp(ixo^s)=1.d0/(sum(w(ixo^s,mag(:))**2,dim=ndim+1)+smalldouble)
542 w(ixo^s,mom(idir))=dxhm*w(ixo^s,mom(idir))*tmp(ixo^s)
545 dxhms(ixo^s)=dble(ndim)/sum(1.d0/
block%dx(ixo^s,:),dim=ndim+1)
547 w(ixo^s,mom(idir))=dxhms(ixo^s)*w(ixo^s,mom(idir))*tmp(ixo^s)
550 vhatmaxgrid=maxval(sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1)))
557 integer,
intent(in) :: ixI^L, ixO^L
558 double precision,
intent(in) :: x(ixI^S,1:ndim),qdt,qvmax
559 double precision,
intent(inout) :: w(ixI^S,1:nw)
561 double precision :: dxhm,disbd(6),bfzone^D
562 double precision :: dxhms(ixO^S)
563 integer :: ix^D, idir
567 dxhm=dble(ndim)/(^d&1.0d0/
dxlevel(^d)+)
570 w(ixo^s,mom(:))=w(ixo^s,mom(:))*dxhm
572 dxhms(ixo^s)=dble(ndim)/sum(1.d0/
block%dx(ixo^s,:),dim=ndim+1)
573 dxhms(ixo^s)=
mf_cc*
mf_cy/qvmax*dxhms(ixo^s)/qdt
576 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*dxhms(ixo^s)
581 bfzone1=0.05d0*(xprobmax1-xprobmin1)
582 bfzone2=0.05d0*(xprobmax2-xprobmin2)
583 bfzone3=0.05d0*(xprobmax3-xprobmin3)
584 {
do ix^db=ixomin^db,ixomax^db\}
585 disbd(1)=x(ix^d,1)-xprobmin1
586 disbd(2)=xprobmax1-x(ix^d,1)
587 disbd(3)=x(ix^d,2)-xprobmin2
588 disbd(4)=xprobmax2-x(ix^d,2)
589 disbd(5)=x(ix^d,3)-xprobmin1
590 disbd(6)=xprobmax3-x(ix^d,3)
593 if(disbd(1)<bfzone1)
then
594 w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(1))/bfzone1)**2)*w(ix^d,mom(:))
597 if(disbd(5)<bfzone3)
then
598 w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(5))/bfzone3)**2)*w(ix^d,mom(:))
601 if(disbd(2)<bfzone1)
then
602 w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(2))/bfzone1)**2)*w(ix^d,mom(:))
604 if(disbd(3)<bfzone2)
then
605 w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(3))/bfzone2)**2)*w(ix^d,mom(:))
607 if(disbd(4)<bfzone2)
then
608 w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(4))/bfzone2)**2)*w(ix^d,mom(:))
610 if(disbd(6)<bfzone3)
then
611 w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(6))/bfzone3)**2)*w(ix^d,mom(:))
617 dxhms(ixo^s)=sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1))/
mf_vmax+1.d-12
618 dxhms(ixo^s)=dtanh(dxhms(ixo^s))/dxhms(ixo^s)
620 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*dxhms(ixo^s)
632 integer,
intent(in) :: idim^LIM
633 double precision,
intent(in) :: qt, qdt
635 integer :: iigrid, igrid
640 do iigrid=1,igridstail; igrid=igrids(iigrid);
641 ps1(igrid)%w=ps(igrid)%w
660 do iigrid=1,igridstail; igrid=igrids(iigrid);
661 ps2(igrid)%w(ixg^t,1:nwflux)=0.75d0*ps(igrid)%w(ixg^t,1:nwflux)+0.25d0*&
662 ps1(igrid)%w(ixg^t,1:nwflux)
663 if (nw>nwflux) ps2(igrid)%w(ixg^t,nwflux+1:nw) = &
664 ps(igrid)%w(ixg^t,nwflux+1:nw)
669 do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
670 ps(igrid)%w(ixg^t,1:nwflux)=1.0d0/3.0d0*ps(igrid)%w(ixg^t,1:nwflux)+&
671 2.0d0/3.0d0*ps2(igrid)%w(ixg^t,1:nwflux)
676 call mpistop(
"unkown time_stepper in advectmf")
681 subroutine advect1mf(method,dtin,dtfactor,idim^LIM,qtC,psa,qt,psb)
689 integer,
intent(in) :: idim^LIM
690 type(state) :: psa(max_blocks)
691 type(state) :: psb(max_blocks)
692 double precision,
intent(in) :: dtin,dtfactor, qtC, qt
693 integer,
intent(in) :: method(nlevelshi)
695 double precision :: qdt
696 integer :: iigrid, igrid, level, i^D
703 do iigrid=1,igridstail; igrid=igrids(iigrid);
708 psa(igrid)%w,qt,psb(igrid)%w)
748 integer,
intent(in) :: method
749 integer,
intent(in) :: igrid, ixG^L, idim^LIM
750 double precision,
intent(in) :: qdt, qtC, qt
751 double precision :: wCT(ixG^S,1:nw), w(ixG^S,1:nw)
752 double precision :: dx^D, fC(ixG^S,1:ndir,1:ndim)
759 ixo^l=ixg^l^lsubnghostcells;
765 call centdiff4mf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,fc,dx^d,ps(igrid)%x)
770 call tvdlfmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,fc,dx^d,ps(igrid)%x)
773 call hancockmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,dx^d,ps(igrid)%x)
778 call fdmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,fc,dx^d,ps(igrid)%x)
780 call mpistop(
"unknown flux scheme in advect1_gridmf")
789 subroutine upwindlrmf(ixI^L,ixL^L,ixR^L,idim,w,wCT,wLC,wRC,x)
795 integer,
intent(in) :: ixI^L, ixL^L, ixR^L, idim
796 double precision,
dimension(ixI^S,1:nw) :: w, wCT
797 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC
798 double precision,
dimension(ixI^S,1:ndim) :: x
800 double precision :: ldw(ixI^S), rdw(ixI^S), dwC(ixI^S)
801 integer :: jxR^L, ixC^L, jxC^L, iw
804 call mp5limiter(ixi^l,ixl^l,idim,w,wlc,wrc)
806 call ppmlimiter(ixi^l,
ixm^
ll,idim,w,wct,wlc,wrc)
808 jxr^l=ixr^l+
kr(idim,^
d);
809 ixcmax^
d=jxrmax^
d; ixcmin^
d=ixlmin^
d-
kr(idim,^
d);
810 jxc^l=ixc^l+
kr(idim,^
d);
814 w(ixcmin^
d:jxcmax^
d,iw)=dlog10(w(ixcmin^
d:jxcmax^
d,iw))
815 wlc(ixl^s,iw)=dlog10(wlc(ixl^s,iw))
816 wrc(ixr^s,iw)=dlog10(wrc(ixr^s,iw))
819 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
823 wlc(ixl^s,iw)=wlc(ixl^s,iw)+half*ldw(ixl^s)
824 wrc(ixr^s,iw)=wrc(ixr^s,iw)-half*rdw(jxr^s)
827 w(ixcmin^
d:jxcmax^
d,iw)=10.0d0**w(ixcmin^
d:jxcmax^
d,iw)
828 wlc(ixl^s,iw)=10.0d0**wlc(ixl^s,iw)
829 wrc(ixr^s,iw)=10.0d0**wrc(ixr^s,iw)
841 integer,
intent(in) :: ixI^L, ixO^L, idir, idim
842 double precision,
intent(in) :: w(ixI^S,nw)
843 double precision,
intent(in) :: x(ixI^S,1:ndim)
844 double precision,
intent(out) :: f(ixI^S)
851 f(ixo^s)=w(ixo^s,mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom(idir))
854 +w(ixo^s,mom(idim))*
block%B0(ixo^s,idir,idim)&
855 -w(ixo^s,mom(idir))*
block%B0(ixo^s,idim,idim)
861 subroutine tvdlfmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,fC,dx^D,x)
867 double precision,
intent(in) :: qdt, qtC, qt, dx^D
868 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
869 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
870 double precision,
dimension(ixI^S,1:nw) :: wCT, wnew
871 double precision,
dimension(ixI^S,1:ndir,1:ndim) :: fC
873 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC, wmean
874 double precision,
dimension(ixI^S) :: fLC, fRC
875 double precision,
dimension(ixI^S) :: cmaxC
876 double precision :: dxinv(1:ndim), inv_volume(ixO^S)
877 integer :: idims, idir, ix^L, hxO^L, ixC^L, ixCR^L, jxC^L, kxC^L, kxR^L
883 ix^l=ix^l^ladd2*
kr(idims,^d);
885 if (ixi^l^ltix^l|.or.|.or.) &
886 call mpistop(
"Error in tvdlfmf: Nonconforming input limits")
888 ^d&dxinv(^d)=-qdt/dx^d;
893 hxo^l=ixo^l-
kr(idims,^d);
895 ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
897 jxc^l=ixc^l+
kr(idims,^d);
898 kxcmin^d=iximin^d; kxcmax^d=iximax^d-
kr(idims,^d);
899 kxr^l=kxc^l+
kr(idims,^d);
902 wrc(kxc^s,1:nwflux)=wct(kxr^s,1:nwflux)
903 wlc(kxc^s,1:nwflux)=wct(kxc^s,1:nwflux)
905 call upwindlrmf(ixi^l,ixcr^l,ixcr^l,idims,wct,wct,wlc,wrc,x)
910 wmean=0.5d0*(wlc+wrc)
918 flc(ixc^s)=half*(flc(ixc^s)+frc(ixc^s))
921 if (idir==idims)
then
922 flc(ixc^s)=flc(ixc^s)-
mf_tvdlfeps*
tvdlfeps*cmaxc(ixc^s)*half*(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
925 fc(ixc^s,idir,idims)=flc(ixc^s)
927 fc(ixc^s,idir,idims)=
block%surfaceC(ixc^s,idims)*flc(ixc^s)
936 hxo^l=ixo^l-
kr(idims,^d);
939 fc(ixi^s,:,idims)=dxinv(idims)*fc(ixi^s,:,idims)
940 wnew(ixo^s,mag(:))=wnew(ixo^s,mag(:)) &
941 + (fc(ixo^s,:,idims)-fc(hxo^s,:,idims))
943 inv_volume = 1.0d0/
block%dvolume(ixo^s)
944 fc(ixi^s,:,idims)=-qdt*fc(ixi^s,:,idims)
947 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir)) + (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims)) * &
955 call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
959 subroutine hancockmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,dx^D,x)
968 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
969 double precision,
intent(in) :: qdt, qtC, qt, dx^D, x(ixI^S,1:ndim)
970 double precision,
intent(inout) :: wCT(ixI^S,1:nw), wnew(ixI^S,1:nw)
972 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC
973 double precision,
dimension(ixI^S) :: fLC, fRC
974 double precision :: dxinv(1:ndim)
975 integer :: idims, idir, ix^L, hxO^L, ixtest^L
980 ix^l=ix^l^laddkr(idims,^d);
982 if (ixi^l^ltix^l|.or.|.or.) &
983 call mpistop(
"Error in Hancockmf: Nonconforming input limits")
985 ^d&dxinv(^d)=-qdt/dx^d;
991 hxo^l=ixo^l-
kr(idims,^d);
993 wrc(hxo^s,1:nwflux)=wct(ixo^s,1:nwflux)
994 wlc(ixo^s,1:nwflux)=wct(ixo^s,1:nwflux)
996 call upwindlrmf(ixi^l,ixo^l,hxo^l,idims,wct,wct,wlc,wrc,x)
1001 call getfluxmf(wrc,x,ixi^l,hxo^l,idir,idims,frc)
1002 call getfluxmf(wlc,x,ixi^l,ixo^l,idir,idims,flc)
1005 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+dxinv(idims)* &
1006 (flc(ixo^s)-frc(hxo^s))
1008 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))-qdt/
block%dvolume(ixo^s) &
1009 *(
block%surfaceC(ixo^s,idims)*flc(ixo^s) &
1010 -
block%surfaceC(hxo^s,idims)*frc(hxo^s))
1020 subroutine fdmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,fC,dx^D,x)
1022 double precision,
intent(in) :: qdt, qtC, qt, dx^D
1023 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
1024 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
1026 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: wCT, wnew
1027 double precision,
dimension(ixI^S,1:ndir,1:ndim),
intent(out) :: fC
1029 double precision,
dimension(ixI^S) :: fCT
1030 double precision,
dimension(ixI^S,1:nw) :: fm, fp, fmR, fpL
1031 double precision,
dimension(ixI^S) :: v
1032 double precision :: dxinv(1:ndim)
1033 integer :: idims, idir, ixC^L, ix^L, hxO^L, ixCR^L
1035 ^d&dxinv(^d)=-qdt/dx^d;
1042 hxo^l=ixo^l-
kr(idims,^d);
1044 ixmax^d=ixomax^d; ixmin^d=hxomin^d;
1048 call getfluxmf(wct,x,ixg^
ll,ixcr^l,idir,idims,fct)
1060 fc(ix^s,idir,idims) = dxinv(idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
1061 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
1062 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1064 fc(ix^s,idir,idims)=-qdt*
block%surfaceC(ix^s,idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
1065 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
1066 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/
block%dvolume(ixo^s)
1073 call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
1081 integer,
intent(in) :: ixI^L, iL^L, idims
1082 double precision,
intent(in) :: w(ixI^S,1:nw)
1084 double precision,
intent(out) :: wLC(ixI^S,1:nw)
1086 double precision :: ldw(ixI^S), dwC(ixI^S)
1087 integer :: jxR^L, ixC^L, jxC^L, kxC^L, iw
1091 call mp5limiterl(ixi^l,il^l,idims,w,wlc)
1093 call weno5limiterl(ixi^l,il^l,idims,w,wlc,1)
1095 call weno5limiterl(ixi^l,il^l,idims,w,wlc,2)
1098 kxcmin^
d=iximin^
d; kxcmax^
d=iximax^
d-
kr(idims,^
d);
1100 wlc(kxc^s,1:nwflux) = w(kxc^s,1:nwflux)
1102 jxr^l=il^l+
kr(idims,^
d);
1104 ixcmax^
d=jxrmax^
d; ixcmin^
d=ilmin^
d-
kr(idims,^
d);
1105 jxc^l=ixc^l+
kr(idims,^
d);
1108 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1112 wlc(il^s,iw)=wlc(il^s,iw)+half*ldw(il^s)
1122 integer,
intent(in) :: ixI^L, iL^L, idims
1123 double precision,
intent(in) :: w(ixI^S,1:nw)
1125 double precision,
intent(out) :: wRC(ixI^S,1:nw)
1127 double precision :: rdw(ixI^S), dwC(ixI^S)
1128 integer :: jxR^L, ixC^L, jxC^L, kxC^L, kxR^L, iw
1132 call mp5limiterr(ixi^l,il^l,idims,w,wrc)
1134 call weno5limiterr(ixi^l,il^l,idims,w,wrc,1)
1136 call weno5limiterr(ixi^l,il^l,idims,w,wrc,2)
1139 kxcmin^
d=iximin^
d; kxcmax^
d=iximax^
d-
kr(idims,^
d);
1140 kxr^l=kxc^l+
kr(idims,^
d);
1142 wrc(kxc^s,1:nwflux)=w(kxr^s,1:nwflux)
1144 jxr^l=il^l+
kr(idims,^
d);
1145 ixcmax^
d=jxrmax^
d; ixcmin^
d=ilmin^
d-
kr(idims,^
d);
1146 jxc^l=ixc^l+
kr(idims,^
d);
1149 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1152 wrc(il^s,iw)=wrc(il^s,iw)-half*rdw(jxr^s)
1158 subroutine centdiff4mf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,w,fC,dx^D,x)
1167 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
1168 double precision,
intent(in) :: qdt, qtC, qt, dx^D
1169 double precision :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
1170 double precision,
intent(in) :: x(ixI^S,1:ndim)
1171 double precision :: fC(ixI^S,1:ndir,1:ndim)
1173 double precision :: v(ixI^S,ndim), f(ixI^S)
1174 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC
1175 double precision,
dimension(ixI^S) :: vLC, vRC,cmaxLC,cmaxRC
1176 double precision :: dxinv(1:ndim)
1177 integer :: idims, idir, idirmin,ix^D
1178 integer :: ix^L, hxO^L, ixC^L, jxC^L, hxC^L, kxC^L, kkxC^L, kkxR^L
1183 ix^l=ix^l^ladd2*
kr(idims,^d);
1186 if (ixi^l^ltix^l|.or.|.or.)
then
1187 call mpistop(
"Error in evolve_CentDiff4: Non-conforming input limits")
1189 ^d&dxinv(^d)=-qdt/dx^d;
1194 ix^l=ixo^l^ladd2*
kr(idims,^d);
1195 hxo^l=ixo^l-
kr(idims,^d);
1197 ixcmin^d=hxomin^d; ixcmax^d=ixomax^d;
1198 hxc^l=ixc^l-
kr(idims,^d);
1199 jxc^l=ixc^l+
kr(idims,^d);
1200 kxc^l=ixc^l+2*
kr(idims,^d);
1202 kkxcmin^d=iximin^d; kkxcmax^d=iximax^d-
kr(idims,^d);
1203 kkxr^l=kkxc^l+
kr(idims,^d);
1204 wrc(kkxc^s,1:nwflux)=wct(kkxr^s,1:nwflux)
1205 wlc(kkxc^s,1:nwflux)=wct(kkxc^s,1:nwflux)
1207 call upwindlrmf(ixi^l,ixc^l,ixc^l,idims,wct,wct,wlc,wrc,x)
1213 vlc(ixc^s)=max(cmaxrc(ixc^s),cmaxlc(ixc^s))
1217 call getfluxmf(wct,x,ixi^l,ix^l,idir,idims,f)
1220 fc(ixc^s,idir,idims)=(-f(kxc^s)+7.0d0*(f(jxc^s)+f(ixc^s))-f(hxc^s))/12.0d0
1225 *(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
1228 fc(ixc^s,idir,idims)=dxinv(idims)*fc(ixc^s,idir,idims)
1230 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+(fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1232 fc(ixc^s,idir,idims)=-qdt*
block%surfaceC(ixc^s,idims)*fc(ixc^s,idir,idims)
1233 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+ &
1234 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/
block%dvolume(ixo^s)
1249 integer,
intent(in) :: ixI^L, ixO^L
1250 double precision,
intent(in) :: x(ixI^S,1:ndim)
1251 double precision,
intent(inout) :: w(ixI^S,1:nw), dtnew
1253 double precision :: courantmax, dxinv(1:ndim)
1254 double precision :: cmax(ixI^S),tmp(ixI^S),alfven(ixI^S)
1265 tmp(ixo^s)=cmax(ixo^s)/
block%dx(ixo^s,idims)
1266 courantmax=max(courantmax,maxval(tmp(ixo^s)))
1268 tmp(ixo^s)=cmax(ixo^s)*dxinv(idims)
1269 courantmax=max(courantmax,maxval(tmp(ixo^s)))
1273 if (courantmax>smalldouble) dtnew=min(dtnew,
mf_cc/courantmax)
1280 logical :: new_cmax,needcmin
1281 integer,
intent(in) :: ixI^L, ixO^L, idims
1282 double precision,
intent(in) :: w(ixI^S,1:nw)
1283 double precision,
intent(out) :: cmax(ixI^S)
1287 cmax(ixo^s)=sqrt(sum((w(ixo^s,mag(:))+
block%b0(ixo^s,:,0))**2,dim=
ndim+1)/w(ixo^s,rho_))
1289 cmax(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2,dim=
ndim+1)/w(ixo^s,rho_))
1291 cmax(ixo^s)=cmax(ixo^s)+abs(w(ixo^s,mom(idims)))
1300 integer,
intent(in) :: ixI^L, ixO^L
1301 double precision,
intent(in) :: x(ixI^S,1:ndim),wCT(ixI^S,1:nw),qdt
1302 double precision,
intent(inout) :: w(ixI^S,1:nw)
1303 double precision :: divb(ixI^S),graddivb(ixI^S),bdivb(ixI^S,1:ndir)
1304 integer :: idims, ix^L, ixp^L, i^D, iside
1315 call gradient(divb,ixi^l,ixp^l,idims,graddivb)
1321 graddivb(ixp^s)=graddivb(ixp^s)*
mf_cdivb &
1322 /(^d&1.0d0/
block%dx(ixp^s,^d)**2+)
1329 w(ixp^s,mag(idims))=w(ixp^s,mag(idims))+&
1340 integer,
intent(in) :: ixI^L, ixO^L
1341 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
1342 double precision,
intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
1344 double precision :: tmp(ixI^S)
1346 integer :: mr_,mphi_
1347 integer :: br_,bphi_
1349 mr_=mom(1); mphi_=mom(1)-1+
phi_
1350 br_=mag(1); bphi_=mag(1)-1+
phi_
1356 tmp(ixo^s)=(wct(ixo^s,bphi_)*wct(ixo^s,mom(1)) &
1357 -wct(ixo^s,br_)*wct(ixo^s,mom(3)))
1358 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt*tmp(ixo^s)/x(ixo^s,1)
1364 tmp(ixo^s)= wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
1365 -wct(ixo^s,mom(2))*wct(ixo^s,mag(1))
1367 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*
block%b0(ixo^s,2,0) &
1368 -wct(ixo^s,mom(2))*
block%b0(ixo^s,1,0)
1371 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
1376 tmp(ixo^s)=wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
1377 -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)){^nooned &
1378 -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
1379 -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
1382 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*
block%b0(ixo^s,3,0) &
1383 -wct(ixo^s,mom(3))*
block%b0(ixo^s,1,0){^nooned &
1384 -(wct(ixo^s,mom(3))*
block%b0(ixo^s,2,0) &
1385 -wct(ixo^s,mom(2))*
block%b0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
1389 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
1402 integer :: ixO^L, idirmin, ixI^L
1403 double precision :: w(ixI^S,1:nw)
1406 double precision :: current(ixI^S,7-2*ndir:3),bvec(ixI^S,1:ndir)
1411 bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
1413 call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1415 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
1416 block%J0(ixo^s,idirmin0:3)
1425 integer,
intent(in) :: ixI^L, ixO^L
1426 double precision,
intent(in) :: w(ixI^S,1:nw)
1427 double precision :: divb(ixI^S)
1429 double precision :: bvec(ixI^S,1:ndir)
1431 bvec(ixi^s,:)=w(ixi^s,mag(:))
subroutine mask_inner(ixIL, ixOL, w, x)
double precision function integral_grid_mf(ixIL, ixOL, w, x, iw, patchwi)
subroutine, public resettree
reset AMR and (de)allocate boundary flux storage at level changes
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Module for flux conservation near refinement boundaries.
subroutine, public recvflux(idimLIM)
subroutine, public sendflux(idimLIM)
subroutine, public fix_conserve(psb, idimLIM, nw0, nwfluxin)
subroutine, public init_comm_fix_conserve(idimLIM, nwfluxin)
subroutine, public store_flux(igrid, fC, idimLIM, nwfluxin)
Module with geometry-related routines (e.g., divergence, curl)
subroutine divvectors(qvec, ixIL, ixOL, divq)
Calculate divergence of a vector qvec within ixL using limited extrapolation to cell edges.
integer, parameter spherical
integer, parameter cylindrical
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
subroutine curlvector(qvec, ixIL, ixOL, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine divvector(qvec, ixIL, ixOL, divq, fourthorder, sixthorder)
Calculate divergence of a vector qvec within ixL.
update ghost cells of all blocks including physical boundaries
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_f
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_p1
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_p2
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_p2
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_p1
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_p1
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_f
integer, dimension(:^d &,:^d &), pointer type_recv_srl
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_p2
subroutine create_bc_mpi_datatype(nwstart, nwbc)
integer, dimension(-1:2^d &,-1:1^d &), target type_recv_srl_p1
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_p2
integer, dimension(:^d &,:^d &), pointer type_send_srl
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_p2
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_f
integer, dimension(:^d &,:^d &), pointer type_recv_r
integer, dimension(:^d &,:^d &), pointer type_send_r
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_r_p1
integer, dimension(-1:1^d &, 0:3^d &), target type_recv_p_f
subroutine getbc(time, qdt, psb, nwstart, nwbc, req_diag)
do update ghost cells of all blocks including physical boundaries
integer, dimension(-1:1^d &, 0:3^d &), target type_send_p_p2
integer, dimension(-1:2^d &,-1:1^d &), target type_send_srl_f
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_p1
integer, dimension(:^d &,:^d &), pointer type_recv_p
integer, dimension(:^d &,:^d &), pointer type_send_p
integer, dimension(-1:1^d &,-1:1^d &), target type_send_r_f
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 tvdlfeps
integer, dimension(:), allocatable typepred1
The spatial discretization for the predictor step when using a two step PC method.
integer, parameter fs_tvdlf
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
integer, parameter unitpar
file handle for IO
double precision global_time
The global simulation time.
integer istep
Index of the sub-step in a multi-step time integrator.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer, parameter threestep
integer snapshotini
Resume from the snapshot with this index.
integer it
Number of time steps taken.
logical, dimension(:), allocatable loglimit
integer ditregrid
Reconstruct the AMR grid once every ditregrid iteration(s)
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
integer b0i
background magnetic field location indicator
integer mype
The rank of the current MPI task.
character(len=std_len) typediv
integer, parameter plevel_
double precision, dimension(:), allocatable, parameter d
double precision dt
global time step
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.
integer, dimension(:), allocatable flux_method
Which flux scheme of spatial discretization to use (per grid level)
logical slab
Cartesian geometry or not.
double precision unit_velocity
Physical scaling factor for velocity.
logical prolongprimitive
prolongate primitive variables in level-jump ghost cells
logical b0field
split magnetic field as background B0 field
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer, dimension(:), allocatable type_limiter
Type of slope limiter used for reconstructing variables on cell edges.
integer, parameter twostep
integer, parameter onestep
integer nghostcells
Number of ghost cells surrounding a grid.
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
integer t_stepper
time stepper type
integer, parameter fs_cd4
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
integer refine_max_level
Maximal number of AMR levels.
integer, parameter fs_hancock
integer, dimension(:,:), allocatable node
Module with slope/flux limiters.
integer, parameter limiter_ppm
subroutine dwlimiter2(dwC, ixIL, ixCL, idims, typelim, ldw, rdw, a2max)
Limit the centered dwC differences within ixC for iw in direction idim. The limiter is chosen accordi...
integer, parameter limiter_weno5
integer, parameter limiter_wenoz5
integer, parameter limiter_mp5
module mod_magnetofriction.t Purpose: use magnetofrictional method to relax 3D magnetic field to forc...
subroutine reconstructrmf(ixIL, iLL, idims, w, wRC)
subroutine advect1mf(method, dtin, dtfactor, idimLIM, qtC, psa, qt, psb)
subroutine getdtfff_courant(w, x, ixIL, ixOL, dtnew)
subroutine process1_gridmf(method, igrid, qdt, ixGL, idimLIM, qtC, wCT, qt, w)
subroutine divbclean(qdt, ixIL, ixOL, wCT, w, x)
Clean divergence of magnetic field by Janhunen's and Linde's source terms.
subroutine getfluxmf(w, x, ixIL, ixOL, idir, idim, f)
subroutine magnetofriction
subroutine frictional_velocity(w, x, ixIL, ixOL, qvmax, qdt)
subroutine get_divb(w, ixIL, ixOL, divb)
Calculate div B within ixO.
double precision mf_tvdlfeps
TVDLF dissipation coefficient controls the dissipation term.
subroutine magnetofriction_init()
Initialize the module.
subroutine vhat(w, x, ixIL, ixOL, vhatmaxgrid)
subroutine advectmf(idimLIM, qt, qdt)
double precision, public mf_vmax
maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
double precision tmf
time in magnetofriction process
double precision mf_cy_max
subroutine tvdlfmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, fC, dxD, x)
subroutine fdmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, fC, dxD, x)
subroutine centdiff4mf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, w, fC, dxD, x)
double precision mf_cdivb
divb cleaning coefficient controls diffusion speed of divb
subroutine getcmaxfff(w, ixIL, ixOL, idims, cmax)
subroutine mf_velocity_update(dtfff)
double precision cmax_mype
maximal speed for fd scheme
double precision mf_cy
frictional velocity coefficient
logical fix_conserve_at_step
double precision mf_tvdlfeps_min
double precision mf_cdivb_max
double precision mf_cc
stability coefficient controls numerical stability
subroutine mf_params_read(files)
Read this module"s parameters from a file.
subroutine hancockmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, dxD, x)
double precision cmax_global
maximal speed for fd scheme
subroutine upwindlrmf(ixIL, ixLL, ixRL, idim, w, wCT, wLC, wRC, x)
subroutine addgeometrymf(qdt, ixIL, ixOL, wCT, w, x)
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)...
subroutine reconstructlmf(ixIL, iLL, idims, w, wLC)
This module defines the procedures of a physics module. It contains function pointers for the various...
procedure(sub_convert), pointer phys_to_primitive
procedure(sub_convert), pointer phys_to_conserved
Module with all the methods that users can customize in AMRVAC.
procedure(p_no_args), pointer usr_before_main_loop