38 double precision ::
tmf
45 integer,
private,
protected :: rho_
48 integer,
allocatable,
private,
protected :: mom(:)
51 integer,
allocatable,
private,
protected :: mag(:)
63 character(len=*),
intent(in) :: files(:)
70 open(
unitpar, file=trim(files(n)), status=
"old")
115 double precision :: dvolume(ixG^T),dsurface(ixG^T),dvone
116 double precision :: dtfff,dtfff_pe,dtnew,dx^D
117 double precision :: cwsin_theta_new,cwsin_theta_old
118 double precision :: sum_jbb,sum_jbb_ipe,sum_j,sum_j_ipe,sum_l_ipe,sum_l
119 double precision :: f_i_ipe,f_i,volumepe,volume,tmpt,time_in
120 double precision,
external :: integral_grid
121 integer :: i,iigrid, igrid, idims,ix^D,hxM^LL,fhmf,tmpit,i^D
122 logical :: patchwi(ixG^T), stagger_flag
129 if(
mype==0)
write(*,*)
'Evolving to force-free field using magnetofricitonal method...'
161 do iigrid=1,igridstail; igrid=igrids(iigrid);
178 do iigrid=1,igridstail; igrid=igrids(iigrid);
180 pso(igrid)%w(ixg^t,mag(:))=ps(igrid)%w(ixg^t,mag(:))
183 dtfff_pe=min(dtfff_pe,dtnew)
185 call mpi_allreduce(dtfff_pe,dtfff,1,mpi_double_precision,mpi_min, &
210 if(mod(i,10)==0)
then
218 do iigrid=1,igridstail; igrid=igrids(iigrid);
224 do iigrid=1,igridstail; igrid=igrids(iigrid);
230 write(*,*)
'<CW sin theta>:',cwsin_theta_new
231 write(*,*)
'<f_i>:',f_i
232 write(*,*)
'----------------------------------------------------------'
238 if(mod(i,10)/=0)
then
244 write (*,*)
'Reach maximum iteration step!'
245 write (*,*)
'The total iteration step is:', i
259 do iigrid=1,igridstail; igrid=igrids(iigrid);
260 ps(igrid)%w(ixg^t,mom(:))=zero
269 if(
mype==0)
write(*,*)
'Magnetofriction phase took : ',mpi_wtime()-time_in,
' sec'
280 do iigrid=1,igridstail; igrid=igrids(iigrid);
284 hxm^ll=
ixm^ll-
kr(idims,^d);
285 dsurface(
ixm^t)=dsurface(
ixm^t)+
block%surfaceC(hxm^t,idims)
290 ps(igrid)%x,1,patchwi)
292 ps(igrid)%x,2,patchwi)
294 ps(igrid)%x,3,patchwi)
296 ps(igrid)%x,4,patchwi)
298 call mpi_allreduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
300 call mpi_allreduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,&
302 call mpi_allreduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,&
304 call mpi_allreduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,&
306 call mpi_allreduce(volumepe,volume,1,mpi_double_precision,mpi_sum,&
310 cwsin_theta_new = sum_jbb/sum_j
320 integer,
intent(in) :: ixI^L,ixO^L
321 double precision,
intent(in):: w(ixI^S,nw),x(ixI^S,1:ndim)
322 double precision :: xO^L
325 {xomin^d = xprobmin^d + 0.05d0*(xprobmax^d-xprobmin^d)\}
326 {xomax^d = xprobmax^d - 0.05d0*(xprobmax^d-xprobmin^d)\}
328 xomin^nd = xprobmin^nd
333 {
do ix^db=ixomin^db,ixomax^db\}
334 if({ x(ix^dd,^d) > xomin^d .and. x(ix^dd,^d) < xomax^d | .and. })
then
336 volumepe=volumepe+dvolume(ix^d)
338 patchwi(ix^d)=.false.
345 integer :: amode, status(MPI_STATUS_SIZE)
346 character(len=800) :: filename,filehead
347 character(len=2048) :: line,datastr
348 logical,
save :: logmfopened=.false.
351 if(.not.logmfopened)
then
353 write(filename,
"(a,a)") trim(base_filename),
"_mflog.csv"
355 amode=ior(mpi_mode_create,mpi_mode_wronly)
356 amode=ior(amode,mpi_mode_append)
357 call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,ierrmpi)
359 filehead=
" itmf, dt, <f_i>, <CW sin theta>, <Current>, <Lorenz force>"
360 call mpi_file_write(fhmf,filehead,len_trim(filehead), &
361 mpi_character,status,ierrmpi)
362 call mpi_file_write(fhmf,achar(10),1,mpi_character,status,ierrmpi)
365 write(datastr,
'(i6,a)') i,
','
366 line=trim(line)//trim(datastr)
367 write(datastr,
'(es13.6,a)') dtfff,
','
368 line=trim(line)//trim(datastr)
369 write(datastr,
'(es13.6,a)') f_i,
','
370 line=trim(line)//trim(datastr)
371 write(datastr,
'(es13.6,a)') cwsin_theta_new,
','
372 line=trim(line)//trim(datastr)
373 write(datastr,
'(es13.6,a)') sum_j,
','
374 line=trim(line)//trim(datastr)
375 write(datastr,
'(es13.6)') sum_l
376 line=trim(line)//trim(datastr)//new_line(
'A')
377 call mpi_file_write(fhmf,line,len_trim(line),mpi_character,status,ierrmpi)
385 integer,
intent(in) :: ixi^l,ixo^l,iw
386 double precision,
intent(in) :: x(ixi^s,1:ndim)
387 double precision,
intent(in) :: w(ixi^s,nw+nwauxio)
388 logical,
intent(in) :: patchwi(ixi^s)
390 double precision,
dimension(ixI^S,1:ndir) :: bvec,qvec,current
392 integer :: ix^d,i,idirmin,idir,jdir,kdir
399 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
401 bvec(ixi^s,:)=w(ixi^s,mag(:))
405 qvec(ixo^s,1:ndir)=zero
406 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
407 if(lvc(idir,jdir,kdir)/=0)
then
408 tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
409 if(lvc(idir,jdir,kdir)==1)
then
410 qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
412 qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
417 {
do ix^db=ixomin^db,ixomax^db\}
419 sum(bvec(ix^d,:)**2))*dvolume(ix^d)
424 {
do ix^db=ixomin^db,ixomax^db\}
432 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
434 bvec(ixi^s,:)=w(ixi^s,mag(:))
436 call divvector(bvec,ixi^l,ixo^l,tmp)
437 {
do ix^db=ixomin^db,ixomax^db\}
439 dvolume(ix^d)**2/sqrt(sum(bvec(ix^d,:)**2))/dsurface(ix^d)
444 bvec(ixi^s,:)=w(ixi^s,mag(:))+block%b0(ixi^s,mag(:),0)
446 bvec(ixi^s,:)=w(ixi^s,mag(:))
448 call curlvector(bvec,ixi^l,ixo^l,current,idirmin,1,ndir)
450 qvec(ixo^s,1:ndir)=zero
451 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
452 if(lvc(idir,jdir,kdir)/=0)
then
453 tmp(ixo^s)=current(ixo^s,jdir)*bvec(ixo^s,kdir)
454 if(lvc(idir,jdir,kdir)==1)
then
455 qvec(ixo^s,idir)=qvec(ixo^s,idir)+tmp(ixo^s)
457 qvec(ixo^s,idir)=qvec(ixo^s,idir)-tmp(ixo^s)
462 {
do ix^db=ixomin^db,ixomax^db\}
474 double precision,
intent(in) :: dtfff
475 integer :: i,iigrid, igrid
476 double precision :: vhatmax,vhatmax_pe,vhatmaxgrid
478 vhatmax_pe=smalldouble
479 do iigrid=1,igridstail; igrid=igrids(iigrid);
482 call vhat(ps(igrid)%w,ps(igrid)%x,ixg^
ll,
ixm^
ll,vhatmaxgrid)
483 vhatmax_pe=max(vhatmax_pe,vhatmaxgrid)
485 call mpi_allreduce(vhatmax_pe,vhatmax,1,mpi_double_precision,mpi_max, &
487 do iigrid=1,igridstail; igrid=igrids(iigrid);
496 subroutine vhat(w,x,ixI^L,ixO^L,vhatmaxgrid)
500 integer,
intent(in) :: ixI^L, ixO^L
501 double precision,
intent(inout) :: w(ixI^S,nw)
502 double precision,
intent(in) :: x(ixI^S,1:ndim)
503 double precision,
intent(out) :: vhatmaxgrid
505 double precision :: current(ixI^S,7-2*ndir:3),tmp(ixI^S),dxhm
506 double precision :: dxhms(ixO^S)
507 integer :: idirmin,idir,jdir,kdir
512 do idir=1,ndir;
do jdir=1,ndir;
do kdir=idirmin,3
513 if(
lvc(idir,jdir,kdir)/=0)
then
515 tmp(ixo^s)=current(ixo^s,jdir)*(w(ixo^s,mag(kdir))+
block%b0(ixo^s,kdir,0))
517 tmp(ixo^s)=current(ixo^s,jdir)*w(ixo^s,mag(kdir))
519 if(
lvc(idir,jdir,kdir)==1)
then
520 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+tmp(ixo^s)
522 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-tmp(ixo^s)
528 tmp(ixo^s)=sum((w(ixo^s,mag(:))+
block%b0(ixo^s,:,0))**2,dim=ndim+1)
530 tmp(ixo^s)=sum(w(ixo^s,mag(:))**2,dim=ndim+1)
536 w(ixo^s,mom(idir))=dxhm*w(ixo^s,mom(idir))/tmp(ixo^s)
539 dxhms(ixo^s)=dble(ndim)/sum(1.d0/
block%dx(ixo^s,:),dim=ndim+1)
541 w(ixo^s,mom(idir))=dxhms(ixo^s)*w(ixo^s,mom(idir))/tmp(ixo^s)
544 vhatmaxgrid=maxval(sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1)))
551 integer,
intent(in) :: ixI^L, ixO^L
552 double precision,
intent(in) :: x(ixI^S,1:ndim),qdt,qvmax
553 double precision,
intent(inout) :: w(ixI^S,1:nw)
555 double precision :: dxhm,disbd(6),bfzone^D
556 double precision :: dxhms(ixO^S)
557 integer :: ix^D, idir
561 dxhm=dble(ndim)/(^d&1.0d0/
dxlevel(^d)+)
564 w(ixo^s,mom(:))=w(ixo^s,mom(:))*dxhm
566 dxhms(ixo^s)=dble(ndim)/sum(1.d0/
block%dx(ixo^s,:),dim=ndim+1)
567 dxhms(ixo^s)=
mf_cc*
mf_cy/qvmax*dxhms(ixo^s)/qdt
570 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*dxhms(ixo^s)
575 bfzone1=0.05d0*(xprobmax1-xprobmin1)
576 bfzone2=0.05d0*(xprobmax2-xprobmin2)
577 bfzone3=0.05d0*(xprobmax3-xprobmin3)
578 {
do ix^db=ixomin^db,ixomax^db\}
579 disbd(1)=x(ix^d,1)-xprobmin1
580 disbd(2)=xprobmax1-x(ix^d,1)
581 disbd(3)=x(ix^d,2)-xprobmin2
582 disbd(4)=xprobmax2-x(ix^d,2)
583 disbd(5)=x(ix^d,3)-xprobmin1
584 disbd(6)=xprobmax3-x(ix^d,3)
587 if(disbd(1)<bfzone1)
then
588 w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(1))/bfzone1)**2)*w(ix^d,mom(:))
591 if(disbd(5)<bfzone3)
then
592 w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(5))/bfzone3)**2)*w(ix^d,mom(:))
595 if(disbd(2)<bfzone1)
then
596 w(ix^d,mom(:))=(1.d0-((bfzone1-disbd(2))/bfzone1)**2)*w(ix^d,mom(:))
598 if(disbd(3)<bfzone2)
then
599 w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(3))/bfzone2)**2)*w(ix^d,mom(:))
601 if(disbd(4)<bfzone2)
then
602 w(ix^d,mom(:))=(1.d0-((bfzone2-disbd(4))/bfzone2)**2)*w(ix^d,mom(:))
604 if(disbd(6)<bfzone3)
then
605 w(ix^d,mom(:))=(1.d0-((bfzone3-disbd(6))/bfzone3)**2)*w(ix^d,mom(:))
618 integer,
intent(in) :: idim^LIM
619 double precision,
intent(in) :: qt, qdt
621 integer :: iigrid, igrid
627 do iigrid=1,igridstail; igrid=igrids(iigrid);
628 ps1(igrid)%w=ps(igrid)%w
647 do iigrid=1,igridstail; igrid=igrids(iigrid);
648 ps2(igrid)%w(ixg^t,1:nwflux)=0.75d0*ps(igrid)%w(ixg^t,1:nwflux)+0.25d0*&
649 ps1(igrid)%w(ixg^t,1:nwflux)
650 if (nw>nwflux) ps2(igrid)%w(ixg^t,nwflux+1:nw) = &
651 ps(igrid)%w(ixg^t,nwflux+1:nw)
656 do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
657 ps(igrid)%w(ixg^t,1:nwflux)=1.0d0/3.0d0*ps(igrid)%w(ixg^t,1:nwflux)+&
658 2.0d0/3.0d0*ps2(igrid)%w(ixg^t,1:nwflux)
663 call mpistop(
"unkown time_stepper in advectmf")
668 subroutine advect1mf(method,dtin,dtfactor,idim^LIM,qtC,psa,qt,psb)
676 integer,
intent(in) :: idim^LIM
677 type(state) :: psa(max_blocks)
678 type(state) :: psb(max_blocks)
679 double precision,
intent(in) :: dtin,dtfactor, qtC, qt
680 integer,
intent(in) :: method(nlevelshi)
682 double precision :: qdt
683 integer :: iigrid, igrid, level, i^D
690 do iigrid=1,igridstail; igrid=igrids(iigrid);
695 psa(igrid)%w,qt,psb(igrid)%w,pso(igrid)%w)
734 integer,
intent(in) :: method
735 integer,
intent(in) :: igrid, ixG^L, idim^LIM
736 double precision,
intent(in) :: qdt, qtC, qt
737 double precision :: wCT(ixG^S,1:nw), w(ixG^S,1:nw), wold(ixG^S,1:nw)
738 double precision :: dx^D, fC(ixG^S,1:ndir,1:ndim)
745 ixo^l=ixg^l^lsubnghostcells;
751 call centdiff4mf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
756 call tvdlfmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
759 call hancockmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,dx^d,ps(igrid)%x)
764 call fdmf(qdt,ixg^l,ixo^l,idim^lim,qtc,wct,qt,w,wold,fc,dx^d,ps(igrid)%x)
766 call mpistop(
"unknown flux scheme in advect1_gridmf")
775 subroutine upwindlrmf(ixI^L,ixL^L,ixR^L,idim,w,wCT,wLC,wRC,x)
781 integer,
intent(in) :: ixI^L, ixL^L, ixR^L, idim
782 double precision,
dimension(ixI^S,1:nw) :: w, wCT
783 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC
784 double precision,
dimension(ixI^S,1:ndim) :: x
786 integer :: jxR^L, ixC^L, jxC^L, iw
787 double precision :: ldw(ixI^S), rdw(ixI^S), dwC(ixI^S)
790 call mp5limiter(ixi^l,ixl^l,idim,w,wlc,wrc)
792 call ppmlimiter(ixi^l,
ixm^
ll,idim,w,wct,wlc,wrc)
794 jxr^l=ixr^l+
kr(idim,^
d);
795 ixcmax^
d=jxrmax^
d; ixcmin^
d=ixlmin^
d-
kr(idim,^
d);
796 jxc^l=ixc^l+
kr(idim,^
d);
800 w(ixcmin^
d:jxcmax^
d,iw)=dlog10(w(ixcmin^
d:jxcmax^
d,iw))
801 wlc(ixl^s,iw)=dlog10(wlc(ixl^s,iw))
802 wrc(ixr^s,iw)=dlog10(wrc(ixr^s,iw))
805 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
809 wlc(ixl^s,iw)=wlc(ixl^s,iw)+half*ldw(ixl^s)
810 wrc(ixr^s,iw)=wrc(ixr^s,iw)-half*rdw(jxr^s)
813 w(ixcmin^
d:jxcmax^
d,iw)=10.0d0**w(ixcmin^
d:jxcmax^
d,iw)
814 wlc(ixl^s,iw)=10.0d0**wlc(ixl^s,iw)
815 wrc(ixr^s,iw)=10.0d0**wrc(ixr^s,iw)
827 integer,
intent(in) :: ixI^L, ixO^L, idir, idim
828 double precision,
intent(in) :: w(ixI^S,nw)
829 double precision,
intent(in) :: x(ixI^S,1:ndim)
830 double precision,
intent(out) :: f(ixI^S)
837 f(ixo^s)=w(ixo^s,mom(idim))*w(ixo^s,mag(idir))-w(ixo^s,mag(idim))*w(ixo^s,mom(idir))
840 +w(ixo^s,mom(idim))*
block%B0(ixo^s,idir,idim)&
841 -w(ixo^s,mom(idir))*
block%B0(ixo^s,idim,idim)
847 subroutine tvdlfmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,wold,fC,dx^D,x)
852 double precision,
intent(in) :: qdt, qtC, qt, dx^D
853 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
854 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
855 double precision,
dimension(ixI^S,1:nw) :: wCT, wnew, wold
856 double precision,
dimension(ixI^S,1:ndir,1:ndim) :: fC
858 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC, wmean
859 double precision,
dimension(ixI^S) :: fLC, fRC
860 double precision,
dimension(ixI^S) :: cmaxC
861 double precision :: dxinv(1:ndim), inv_volume(ixO^S)
862 integer :: idims, idir, ix^L, hxO^L, ixC^L, ixCR^L, jxC^L, kxC^L, kxR^L
868 ix^l=ix^l^ladd2*
kr(idims,^d);
870 if (ixi^l^ltix^l|.or.|.or.) &
871 call mpistop(
"Error in tvdlfmf: Nonconforming input limits")
873 ^d&dxinv(^d)=-qdt/dx^d;
878 hxo^l=ixo^l-
kr(idims,^d);
880 ixcmax^d=ixomax^d; ixcmin^d=hxomin^d;
882 jxc^l=ixc^l+
kr(idims,^d);
883 kxcmin^d=iximin^d; kxcmax^d=iximax^d-
kr(idims,^d);
884 kxr^l=kxc^l+
kr(idims,^d);
887 wrc(kxc^s,1:nwflux)=wct(kxr^s,1:nwflux)
888 wlc(kxc^s,1:nwflux)=wct(kxc^s,1:nwflux)
890 call upwindlrmf(ixi^l,ixcr^l,ixcr^l,idims,wct,wct,wlc,wrc,x)
895 wmean=0.5d0*(wlc+wrc)
903 flc(ixc^s)=half*(flc(ixc^s)+frc(ixc^s))
906 if (idir==idims)
then
907 flc(ixc^s)=flc(ixc^s)-
mf_tvdlfeps*
tvdlfeps*cmaxc(ixc^s)*half*(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
910 fc(ixc^s,idir,idims)=flc(ixc^s)
912 fc(ixc^s,idir,idims)=
block%surfaceC(ixc^s,idims)*flc(ixc^s)
921 hxo^l=ixo^l-
kr(idims,^d);
924 fc(ixi^s,:,idims)=dxinv(idims)*fc(ixi^s,:,idims)
925 wnew(ixo^s,mag(:))=wnew(ixo^s,mag(:)) &
926 + (fc(ixo^s,:,idims)-fc(hxo^s,:,idims))
928 inv_volume = 1.0d0/
block%dvolume(ixo^s)
929 fc(ixi^s,:,idims)=-qdt*fc(ixi^s,:,idims)
932 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir)) + (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims)) * &
940 call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
944 subroutine hancockmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,dx^D,x)
952 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
953 double precision,
intent(in) :: qdt, qtC, qt, dx^D, x(ixI^S,1:ndim)
954 double precision,
intent(inout) :: wCT(ixI^S,1:nw), wnew(ixI^S,1:nw)
956 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC
957 double precision,
dimension(ixI^S) :: fLC, fRC
958 double precision :: dxinv(1:ndim)
959 integer :: idims, idir, ix^L, hxO^L, ixtest^L
964 ix^l=ix^l^laddkr(idims,^d);
966 if (ixi^l^ltix^l|.or.|.or.) &
967 call mpistop(
"Error in Hancockmf: Nonconforming input limits")
969 ^d&dxinv(^d)=-qdt/dx^d;
975 hxo^l=ixo^l-
kr(idims,^d);
977 wrc(hxo^s,1:nwflux)=wct(ixo^s,1:nwflux)
978 wlc(ixo^s,1:nwflux)=wct(ixo^s,1:nwflux)
980 call upwindlrmf(ixi^l,ixo^l,hxo^l,idims,wct,wct,wlc,wrc,x)
985 call getfluxmf(wrc,x,ixi^l,hxo^l,idir,idims,frc)
986 call getfluxmf(wlc,x,ixi^l,ixo^l,idir,idims,flc)
989 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+dxinv(idims)* &
990 (flc(ixo^s)-frc(hxo^s))
992 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))-qdt/
block%dvolume(ixo^s) &
993 *(
block%surfaceC(ixo^s,idims)*flc(ixo^s) &
994 -
block%surfaceC(hxo^s,idims)*frc(hxo^s))
1004 subroutine fdmf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,wnew,wold,fC,dx^D,x)
1006 double precision,
intent(in) :: qdt, qtC, qt, dx^D
1007 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
1008 double precision,
dimension(ixI^S,1:ndim),
intent(in) :: x
1010 double precision,
dimension(ixI^S,1:nw),
intent(inout) :: wCT, wnew, wold
1011 double precision,
dimension(ixI^S,1:ndir,1:ndim),
intent(out) :: fC
1013 double precision,
dimension(ixI^S) :: fCT
1014 double precision,
dimension(ixI^S,1:nw) :: fm, fp, fmR, fpL
1015 double precision,
dimension(ixI^S) :: v
1016 double precision :: dxinv(1:ndim)
1017 integer :: idims, idir, ixC^L, ix^L, hxO^L, ixCR^L
1019 ^d&dxinv(^d)=-qdt/dx^d;
1026 hxo^l=ixo^l-
kr(idims,^d);
1028 ixmax^d=ixomax^d; ixmin^d=hxomin^d;
1032 call getfluxmf(wct,x,ixg^
ll,ixcr^l,idir,idims,fct)
1044 fc(ix^s,idir,idims) = dxinv(idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
1045 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
1046 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1048 fc(ix^s,idir,idims)=-qdt*
block%surfaceC(ix^s,idims) * (fpl(ix^s,mag(idir)) + fmr(ix^s,mag(idir)))
1049 wnew(ixo^s,mag(idir))=wnew(ixo^s,mag(idir))+ &
1050 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/
block%dvolume(ixo^s)
1057 call divbclean(qdt,ixi^l,ixo^l,wct,wnew,x)
1065 integer,
intent(in) :: ixI^L, iL^L, idims
1066 double precision,
intent(in) :: w(ixI^S,1:nw)
1068 double precision,
intent(out) :: wLC(ixI^S,1:nw)
1070 double precision :: ldw(ixI^S), dwC(ixI^S)
1071 integer :: jxR^L, ixC^L, jxC^L, kxC^L, iw
1075 call mp5limiterl(ixi^l,il^l,idims,w,wlc)
1077 call weno5limiterl(ixi^l,il^l,idims,w,wlc,1)
1079 call weno5limiterl(ixi^l,il^l,idims,w,wlc,2)
1082 kxcmin^
d=iximin^
d; kxcmax^
d=iximax^
d-
kr(idims,^
d);
1084 wlc(kxc^s,1:nwflux) = w(kxc^s,1:nwflux)
1086 jxr^l=il^l+
kr(idims,^
d);
1088 ixcmax^
d=jxrmax^
d; ixcmin^
d=ilmin^
d-
kr(idims,^
d);
1089 jxc^l=ixc^l+
kr(idims,^
d);
1092 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1096 wlc(il^s,iw)=wlc(il^s,iw)+half*ldw(il^s)
1106 integer,
intent(in) :: ixI^L, iL^L, idims
1107 double precision,
intent(in) :: w(ixI^S,1:nw)
1109 double precision,
intent(out) :: wRC(ixI^S,1:nw)
1111 double precision :: rdw(ixI^S), dwC(ixI^S)
1112 integer :: jxR^L, ixC^L, jxC^L, kxC^L, kxR^L, iw
1116 call mp5limiterr(ixi^l,il^l,idims,w,wrc)
1118 call weno5limiterr(ixi^l,il^l,idims,w,wrc,1)
1120 call weno5limiterr(ixi^l,il^l,idims,w,wrc,2)
1123 kxcmin^
d=iximin^
d; kxcmax^
d=iximax^
d-
kr(idims,^
d);
1124 kxr^l=kxc^l+
kr(idims,^
d);
1126 wrc(kxc^s,1:nwflux)=w(kxr^s,1:nwflux)
1128 jxr^l=il^l+
kr(idims,^
d);
1129 ixcmax^
d=jxrmax^
d; ixcmin^
d=ilmin^
d-
kr(idims,^
d);
1130 jxc^l=ixc^l+
kr(idims,^
d);
1133 dwc(ixc^s)=w(jxc^s,iw)-w(ixc^s,iw)
1136 wrc(il^s,iw)=wrc(il^s,iw)-half*rdw(jxr^s)
1142 subroutine centdiff4mf(qdt,ixI^L,ixO^L,idim^LIM,qtC,wCT,qt,w,wold,fC,dx^D,x)
1150 integer,
intent(in) :: ixI^L, ixO^L, idim^LIM
1151 double precision,
intent(in) :: qdt, qtC, qt, dx^D
1152 double precision :: wCT(ixI^S,1:nw), w(ixI^S,1:nw), wold(ixI^S,1:nw)
1153 double precision,
intent(in) :: x(ixI^S,1:ndim)
1154 double precision :: fC(ixI^S,1:ndir,1:ndim)
1156 double precision :: v(ixI^S,ndim), f(ixI^S)
1157 double precision,
dimension(ixI^S,1:nw) :: wLC, wRC
1158 double precision,
dimension(ixI^S) :: vLC, vRC,cmaxLC,cmaxRC
1159 double precision :: dxinv(1:ndim)
1160 integer :: idims, idir, idirmin,ix^D
1161 integer :: ix^L, hxO^L, ixC^L, jxC^L, hxC^L, kxC^L, kkxC^L, kkxR^L
1166 ix^l=ix^l^ladd2*
kr(idims,^d);
1169 if (ixi^l^ltix^l|.or.|.or.)
then
1170 call mpistop(
"Error in evolve_CentDiff4: Non-conforming input limits")
1172 ^d&dxinv(^d)=-qdt/dx^d;
1177 ix^l=ixo^l^ladd2*
kr(idims,^d);
1178 hxo^l=ixo^l-
kr(idims,^d);
1180 ixcmin^d=hxomin^d; ixcmax^d=ixomax^d;
1181 hxc^l=ixc^l-
kr(idims,^d);
1182 jxc^l=ixc^l+
kr(idims,^d);
1183 kxc^l=ixc^l+2*
kr(idims,^d);
1185 kkxcmin^d=iximin^d; kkxcmax^d=iximax^d-
kr(idims,^d);
1186 kkxr^l=kkxc^l+
kr(idims,^d);
1187 wrc(kkxc^s,1:nwflux)=wct(kkxr^s,1:nwflux)
1188 wlc(kkxc^s,1:nwflux)=wct(kkxc^s,1:nwflux)
1190 call upwindlrmf(ixi^l,ixc^l,ixc^l,idims,wct,wct,wlc,wrc,x)
1196 vlc(ixc^s)=max(cmaxrc(ixc^s),cmaxlc(ixc^s))
1200 call getfluxmf(wct,x,ixi^l,ix^l,idir,idims,f)
1203 fc(ixc^s,idir,idims)=(-f(kxc^s)+7.0d0*(f(jxc^s)+f(ixc^s))-f(hxc^s))/12.0d0
1208 *(wrc(ixc^s,mag(idir))-wlc(ixc^s,mag(idir)))
1211 fc(ixc^s,idir,idims)=dxinv(idims)*fc(ixc^s,idir,idims)
1213 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+(fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))
1215 fc(ixc^s,idir,idims)=-qdt*
block%surfaceC(ixc^s,idims)*fc(ixc^s,idir,idims)
1216 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+ &
1217 (fc(ixo^s,idir,idims)-fc(hxo^s,idir,idims))/
block%dvolume(ixo^s)
1232 integer,
intent(in) :: ixI^L, ixO^L
1233 double precision,
intent(in) :: x(ixI^S,1:ndim)
1234 double precision,
intent(inout) :: w(ixI^S,1:nw), dtnew
1236 double precision :: courantmax, dxinv(1:ndim)
1237 double precision :: cmax(ixI^S),tmp(ixI^S),alfven(ixI^S)
1248 tmp(ixo^s)=cmax(ixo^s)/
block%dx(ixo^s,idims)
1249 courantmax=max(courantmax,maxval(tmp(ixo^s)))
1251 tmp(ixo^s)=cmax(ixo^s)*dxinv(idims)
1252 courantmax=max(courantmax,maxval(tmp(ixo^s)))
1256 if (courantmax>smalldouble) dtnew=min(dtnew,
mf_cc/courantmax)
1263 logical :: new_cmax,needcmin
1264 integer,
intent(in) :: ixI^L, ixO^L, idims
1265 double precision,
intent(in) :: w(ixI^S,1:nw)
1266 double precision,
intent(out) :: cmax(ixI^S)
1270 cmax(ixo^s)=sqrt(sum((w(ixo^s,mag(:))+
block%b0(ixo^s,:,0))**2,dim=
ndim+1)/w(ixo^s,rho_))
1272 cmax(ixo^s)=sqrt(sum(w(ixo^s,mag(:))**2,dim=
ndim+1)/w(ixo^s,rho_))
1274 cmax(ixo^s)=cmax(ixo^s)+abs(w(ixo^s,mom(idims)))
1283 integer,
intent(in) :: ixI^L, ixO^L
1284 double precision,
intent(in) :: x(ixI^S,1:ndim),wCT(ixI^S,1:nw),qdt
1285 double precision,
intent(inout) :: w(ixI^S,1:nw)
1286 integer :: idims, ix^L, ixp^L, i^D, iside
1287 double precision :: divb(ixI^S),graddivb(ixI^S),bdivb(ixI^S,1:ndir)
1298 call gradient(divb,ixi^l,ixp^l,idims,graddivb)
1304 graddivb(ixp^s)=graddivb(ixp^s)*
mf_cdivb &
1305 /(^d&1.0d0/
block%dx(ixp^s,^d)**2+)
1312 w(ixp^s,mag(idims))=w(ixp^s,mag(idims))+&
1323 integer,
intent(in) :: ixI^L, ixO^L
1324 double precision,
intent(in) :: qdt, x(ixI^S,1:ndim)
1325 double precision,
intent(inout) :: wCT(ixI^S,1:nw), w(ixI^S,1:nw)
1327 double precision :: tmp(ixI^S)
1329 integer :: mr_,mphi_
1330 integer :: br_,bphi_
1332 mr_=mom(1); mphi_=mom(1)-1+
phi_
1333 br_=mag(1); bphi_=mag(1)-1+
phi_
1339 tmp(ixo^s)=(wct(ixo^s,bphi_)*wct(ixo^s,mom(1)) &
1340 -wct(ixo^s,br_)*wct(ixo^s,mom(3)))
1341 w(ixo^s,bphi_)=w(ixo^s,bphi_)+qdt*tmp(ixo^s)/x(ixo^s,1)
1347 tmp(ixo^s)= wct(ixo^s,mom(1))*wct(ixo^s,mag(2)) &
1348 -wct(ixo^s,mom(2))*wct(ixo^s,mag(1))
1350 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*
block%b0(ixo^s,2,0) &
1351 -wct(ixo^s,mom(2))*
block%b0(ixo^s,1,0)
1354 w(ixo^s,mag(2))=w(ixo^s,mag(2))+qdt*tmp(ixo^s)/x(ixo^s,1)
1359 tmp(ixo^s)=wct(ixo^s,mom(1))*wct(ixo^s,mag(3)) &
1360 -wct(ixo^s,mom(3))*wct(ixo^s,mag(1)){^nooned &
1361 -(wct(ixo^s,mom(3))*wct(ixo^s,mag(2)) &
1362 -wct(ixo^s,mom(2))*wct(ixo^s,mag(3)))*dcos(x(ixo^s,2)) &
1365 tmp(ixo^s)=tmp(ixo^s)+wct(ixo^s,mom(1))*
block%b0(ixo^s,3,0) &
1366 -wct(ixo^s,mom(3))*
block%b0(ixo^s,1,0){^nooned &
1367 -(wct(ixo^s,mom(3))*
block%b0(ixo^s,2,0) &
1368 -wct(ixo^s,mom(2))*
block%b0(ixo^s,3,0))*dcos(x(ixo^s,2)) &
1372 w(ixo^s,mag(3))=w(ixo^s,mag(3))+qdt*tmp(ixo^s)/x(ixo^s,1)
1385 integer :: ixO^L, idirmin, ixI^L
1386 double precision :: w(ixI^S,1:nw)
1390 double precision :: current(ixI^S,7-2*ndir:3),bvec(ixI^S,1:ndir)
1394 bvec(ixi^s,1:ndir)=w(ixi^s,mag(1:ndir))
1396 call curlvector(bvec,ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1398 if(
b0field) current(ixo^s,idirmin0:3)=current(ixo^s,idirmin0:3)+&
1399 block%J0(ixo^s,idirmin0:3)
1408 integer,
intent(in) :: ixI^L, ixO^L
1409 double precision,
intent(in) :: w(ixI^S,1:nw)
1410 double precision :: divb(ixI^S)
1412 double precision :: bvec(ixI^S,1:ndir)
1414 bvec(ixi^s,:)=w(ixi^s,mag(:))
subroutine resettree
reset AMR and (de)allocate boundary flux storage at level changes
subroutine mpistop(message)
Exit MPI-AMRVAC with an error message.
subroutine mask_inner(ixIL, ixOL, w, x)
double precision function integral_grid_mf(ixIL, ixOL, w, x, iw, patchwi)
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_
integer, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range (within a block with 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.
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.
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
double precision, dimension(ndim) dxlevel
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 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 process1_gridmf(method, igrid, qdt, ixGL, idimLIM, qtC, wCT, qt, w, wold)
subroutine magnetofriction
subroutine frictional_velocity(w, x, ixIL, ixOL, qvmax, qdt)
subroutine fdmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, wold, fC, dxD, x)
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 tmf
time in magnetofriction process
double precision mf_cy_max
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
subroutine tvdlfmf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, wnew, wold, fC, dxD, x)
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 centdiff4mf(qdt, ixIL, ixOL, idimLIM, qtC, wCT, qt, w, wold, fC, dxD, x)
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