MPI-AMRVAC 3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
Loading...
Searching...
No Matches
mod_mf_phys.t
Go to the documentation of this file.
1!> Magnetofriction module
3 use mod_global_parameters, only: std_len
5 use mod_comm_lib, only: mpistop
6
7
8 implicit none
9 private
10
11 !> viscosity coefficient in s cm^-2 for solar corona (Cheung 2012 ApJ)
12 double precision, public :: mf_nu = 1.d-15
13
14 !> maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
15 double precision, public :: mf_vmax = 3.d6
16
17 !> decay scale of frictional velocity near boundaries
18 double precision, public :: mf_decay_scale(2*^nd)=0.d0
19
20 !> GLM-MHD parameter: ratio of the diffusive and advective time scales for div b
21 !> taking values within [0, 1]
22 double precision, public :: mf_glm_alpha = 0.5d0
23
24 !> The resistivity
25 double precision, public :: mf_eta = 0.0d0
26
27 !> The hyper-resistivity
28 double precision, public :: mf_eta_hyper = 0.0d0
29
30 !> Coefficient of diffusive divB cleaning
31 double precision :: divbdiff = 0.8d0
32
33 !> Helium abundance over Hydrogen
34 double precision, public, protected :: he_abundance=0.1d0
35
36 !> Method type in a integer for good performance
37 integer :: type_divb
38
39 !> Indices of the momentum density
40 integer, allocatable, public, protected :: mom(:)
41
42 !> Indices of the momentum density for the form of better vectorization
43 integer, public, protected :: ^c&m^C_
44
45 !> Indices of the magnetic field for the form of better vectorization
46 integer, public, protected :: ^c&b^C_
47
48 !> Indices of the GLM psi
49 integer, public, protected :: psi_
50
51 !> To skip * layer of ghost cells during divB=0 fix for boundary
52 integer, public, protected :: boundary_divbfix_skip(2*^nd)=0
53
54 ! DivB cleaning methods
55 integer, parameter :: divb_none = 0
56 integer, parameter :: divb_multigrid = -1
57 integer, parameter :: divb_glm = 1
58 integer, parameter :: divb_powel = 2
59 integer, parameter :: divb_janhunen = 3
60 integer, parameter :: divb_linde = 4
61 integer, parameter :: divb_lindejanhunen = 5
62 integer, parameter :: divb_lindepowel = 6
63 integer, parameter :: divb_lindeglm = 7
64 integer, parameter :: divb_ct = 8
65
66 !> Whether particles module is added
67 logical, public, protected :: mf_particles = .false.
68
69 !> Whether GLM-MHD is used
70 logical, public, protected :: mf_glm = .false.
71
72 !> Whether divB cleaning sources are added splitting from fluid solver
73 logical, public, protected :: source_split_divb = .false.
74
75 !> MHD fourth order
76 logical, public, protected :: mf_4th_order = .false.
77
78 !> set to true if need to record electric field on cell edges
79 logical, public, protected :: mf_record_electric_field = .false.
80
81 !> Whether divB is computed with a fourth order approximation
82 integer, public, protected :: mf_divb_nth = 1
83
84 !> To control divB=0 fix for boundary
85 logical, public, protected :: boundary_divbfix(2*^nd)=.true.
86
87 !> Use a compact way to add resistivity
88 logical :: compactres = .false.
89
90 !> Add divB wave in Roe solver
91 logical, public :: divbwave = .true.
92
93 !> clean divb in the initial condition
94 logical, public, protected :: clean_initial_divb=.false.
95
96 !> Method type to clean divergence of B
97 character(len=std_len), public, protected :: typedivbfix = 'ct'
98
99 !> Method type of constrained transport
100 character(len=std_len), public, protected :: type_ct = 'average'
101
102 !> Update all equations due to divB cleaning
103 character(len=std_len) :: typedivbdiff = 'all'
104
105
106 ! Public methods
107 public :: mf_phys_init
108 public :: mf_to_conserved
109 public :: mf_to_primitive
110 public :: mf_face_to_center
111 public :: get_divb
112 public :: get_current
113 public :: get_normalized_divb
116 {^nooned
118 }
119
120contains
121
122 !> Read this module's parameters from a file
123 subroutine mf_read_params(files)
125 use mod_particles, only: particles_eta
126 character(len=*), intent(in) :: files(:)
127 integer :: n
128
129 namelist /mf_list/ mf_nu, mf_vmax, mf_decay_scale, &
131 particles_eta, mf_record_electric_field,&
133 typedivbdiff, type_ct, compactres, divbwave, he_abundance, si_unit, &
136
137 do n = 1, size(files)
138 open(unitpar, file=trim(files(n)), status="old")
139 read(unitpar, mf_list, end=111)
140111 close(unitpar)
141 end do
142
143 end subroutine mf_read_params
144
145 !> Write this module's parameters to a snapsoht
146 subroutine mf_write_info(fh)
148 integer, intent(in) :: fh
149 integer, parameter :: n_par = 1
150 double precision :: values(n_par)
151 character(len=name_len) :: names(n_par)
152 integer, dimension(MPI_STATUS_SIZE) :: st
153 integer :: er
154
155 call mpi_file_write(fh, n_par, 1, mpi_integer, st, er)
156
157 names(1) = "nu"
158 values(1) = mf_nu
159 call mpi_file_write(fh, values, n_par, mpi_double_precision, st, er)
160 call mpi_file_write(fh, names, n_par * name_len, mpi_character, st, er)
161 end subroutine mf_write_info
162
163 subroutine mf_phys_init()
165 use mod_physics
166 use mod_particles, only: particles_init, particles_eta, particles_etah
167 {^nooned
169 }
170
171 integer :: itr, idir
172
173 call mf_read_params(par_files)
174
175 physics_type = "mf"
176 phys_energy = .false.
178
179 if(ndim==1) typedivbfix='none'
180 select case (typedivbfix)
181 case ('none')
182 type_divb = divb_none
183 {^nooned
184 case ('multigrid')
185 type_divb = divb_multigrid
186 use_multigrid = .true.
187 mg%operator_type = mg_laplacian
189 }
190 case ('glm')
191 mf_glm = .true.
192 need_global_cmax = .true.
193 type_divb = divb_glm
194 case ('powel', 'powell')
195 type_divb = divb_powel
196 case ('janhunen')
197 type_divb = divb_janhunen
198 case ('linde')
199 type_divb = divb_linde
200 case ('lindejanhunen')
201 type_divb = divb_lindejanhunen
202 case ('lindepowel')
203 type_divb = divb_lindepowel
204 case ('lindeglm')
205 mf_glm = .true.
206 need_global_cmax = .true.
207 type_divb = divb_lindeglm
208 case ('ct')
209 type_divb = divb_ct
210 stagger_grid = .true.
211 case default
212 call mpistop('Unknown divB fix')
213 end select
214
215 ! set velocity field as flux variables
216 allocate(mom(ndir))
217 mom(:) = var_set_momentum(ndir)
218 m^c_=mom(^c);
219
220 ! set magnetic field as flux variables
221 allocate(mag(ndir))
222 mag(:) = var_set_bfield(ndir)
223 b^c_=mag(^c);
224
225 ! start with magnetic field and skip velocity when update ghostcells
226 iwstart=mag(1)
227 allocate(start_indices(number_species),stop_indices(number_species))
228 ! set the index of the first flux variable for species 1
229 start_indices(1)=mag(1)
230
231 if (mf_glm) then
232 psi_ = var_set_fluxvar('psi', 'psi', need_bc=.false.)
233 else
234 psi_ = -1
235 end if
236
237 ! set number of variables which need update ghostcells
238 nwgc=nwflux-ndir
239
240 ! set the index of the last flux variable for species 1
241 stop_indices(1)=nwflux
242
243 ! determine number of stagger variables
244 nws=ndim
245
246 nvector = 2 ! No. vector vars
247 allocate(iw_vector(nvector))
248 iw_vector(1) = mom(1) - 1 ! TODO: why like this?
249 iw_vector(2) = mag(1) - 1 ! TODO: why like this?
250
251 ! Check whether custom flux types have been defined
252 if (.not. allocated(flux_type)) then
253 allocate(flux_type(ndir, nw))
255 else if (any(shape(flux_type) /= [ndir, nw])) then
256 call mpistop("phys_check error: flux_type has wrong shape")
257 end if
258 do idir=1,ndir
259 if(ndim>1) flux_type(idir,mag(idir))=flux_tvdlf
260 end do
261 if(mf_glm .and. ndim>1) flux_type(:,psi_)=flux_tvdlf
262
263 phys_get_dt => mf_get_dt
264 phys_get_cmax => mf_get_cmax
265 phys_get_cbounds => mf_get_cbounds
266 phys_get_flux => mf_get_flux
267 phys_add_source_geom => mf_add_source_geom
268 phys_add_source => mf_add_source
271 phys_check_params => mf_check_params
272 phys_write_info => mf_write_info
273 phys_special_advance => mf_velocity_update
274
275 if(type_divb==divb_glm) then
276 phys_modify_wlr => mf_modify_wlr
277 end if
278
279 ! pass to global variable to record electric field
281
282 ! Initialize particles module
283 if(mf_particles) then
284 call particles_init()
285 ! never allow Hall effects in particles when doing magnetofrictional
286 particles_etah=0.0d0
287 if(mype==0) then
288 write(*,*) '*****Using particles: with mf_eta, mf_eta_hyper :', mf_eta, mf_eta_hyper
289 write(*,*) '*****Using particles: particles_eta, particles_etah :', particles_eta, particles_etah
290 end if
291 end if
292
293 ! if using ct stagger grid, boundary divb=0 is not done here
294 if(stagger_grid) then
295 phys_get_ct_velocity => mf_get_ct_velocity
296 phys_update_faces => mf_update_faces
298 else if(ndim>1) then
299 phys_boundary_adjust => mf_boundary_adjust
300 end if
301
302 {^nooned
303 ! clean initial divb
305 }
306
307 ! derive units from basic units
308 call mf_physical_units()
309
310 end subroutine mf_phys_init
311
312 subroutine mf_check_params
314
315 end subroutine mf_check_params
316
317 subroutine mf_physical_units()
319 double precision :: mp,kb,miu0,c_lightspeed
320 double precision :: a,b
321 ! Derive scaling units
322 if(si_unit) then
323 mp=mp_si
324 kb=kb_si
325 miu0=miu0_si
326 c_lightspeed=c_si
327 else
328 mp=mp_cgs
329 kb=kb_cgs
330 miu0=4.d0*dpi
331 c_lightspeed=const_c
332 end if
333 a=1.d0+4.d0*he_abundance
334 b=2.d0+3.d0*he_abundance
335 if(unit_density/=1.d0 .or. unit_numberdensity/=1.d0) then
336 if(unit_density/=1.d0) then
338 else if(unit_numberdensity/=1.d0) then
340 end if
341 if(unit_temperature/=1.d0) then
345 if(unit_length/=1.d0) then
347 else if(unit_time/=1.d0) then
349 end if
350 else if(unit_magneticfield/=1.d0) then
354 if(unit_length/=1.d0) then
356 else if(unit_time/=1.d0) then
358 end if
359 else if(unit_pressure/=1.d0) then
363 if(unit_length/=1.d0) then
365 else if(unit_time/=1.d0) then
367 end if
368 else if(unit_velocity/=1.d0) then
372 if(unit_length/=1.d0) then
374 else if(unit_time/=1.d0) then
376 end if
377 else if(unit_time/=1.d0) then
382 end if
383 else if(unit_temperature/=1.d0) then
384 ! units of temperature and velocity are dependent
385 if(unit_magneticfield/=1.d0) then
390 if(unit_length/=1.d0) then
392 else if(unit_time/=1.d0) then
394 end if
395 else if(unit_pressure/=1.d0) then
400 if(unit_length/=1.d0) then
402 else if(unit_time/=1.d0) then
404 end if
405 end if
406 else if(unit_magneticfield/=1.d0) then
407 ! units of magnetic field and pressure are dependent
408 if(unit_velocity/=1.d0) then
413 if(unit_length/=1.d0) then
415 else if(unit_time/=1.d0) then
417 end if
418 else if(unit_time/=0.d0) then
424 end if
425 else if(unit_pressure/=1.d0) then
426 if(unit_velocity/=1.d0) then
431 if(unit_length/=1.d0) then
433 else if(unit_time/=1.d0) then
435 end if
436 else if(unit_time/=0.d0) then
442 end if
443 end if
444 ! Additional units needed for the particles
445 c_norm=c_lightspeed/unit_velocity
447 if (.not. si_unit) unit_charge = unit_charge*const_c
449
450 ! get dimensionless mf nu
452
453 ! get dimensionless maximal mf velocity limit
455
456 end subroutine mf_physical_units
457
458 !> Transform primitive variables into conservative ones
459 subroutine mf_to_conserved(ixI^L,ixO^L,w,x)
461 integer, intent(in) :: ixi^l, ixo^l
462 double precision, intent(inout) :: w(ixi^s, nw)
463 double precision, intent(in) :: x(ixi^s, 1:ndim)
464
465 ! nothing to do for mf
466 end subroutine mf_to_conserved
467
468 !> Transform conservative variables into primitive ones
469 subroutine mf_to_primitive(ixI^L,ixO^L,w,x)
471 integer, intent(in) :: ixi^l, ixo^l
472 double precision, intent(inout) :: w(ixi^s, nw)
473 double precision, intent(in) :: x(ixi^s, 1:ndim)
474
475 ! nothing to do for mf
476 end subroutine mf_to_primitive
477
478 !> Calculate cmax_idim=csound+abs(v_idim) within ixO^L
479 subroutine mf_get_cmax(w,x,ixI^L,ixO^L,idim,cmax)
481
482 integer, intent(in) :: ixi^l, ixo^l, idim
483 double precision, intent(in) :: w(ixi^s, nw), x(ixi^s,1:ndim)
484 double precision, intent(inout) :: cmax(ixi^s)
485
486 cmax(ixo^s)=abs(w(ixo^s,mom(idim)))+one
487
488 end subroutine mf_get_cmax
489
490 !> Estimating bounds for the minimum and maximum signal velocities
491 subroutine mf_get_cbounds(wLC,wRC,wLp,wRp,x,ixI^L,ixO^L,idim,Hspeed,cmax,cmin)
493 use mod_variables
494
495 integer, intent(in) :: ixi^l, ixo^l, idim
496 double precision, intent(in) :: wlc(ixi^s, nw), wrc(ixi^s, nw)
497 double precision, intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
498 double precision, intent(in) :: x(ixi^s,1:ndim)
499 double precision, intent(inout) :: cmax(ixi^s,1:number_species)
500 double precision, intent(inout), optional :: cmin(ixi^s,1:number_species)
501 double precision, intent(in) :: hspeed(ixi^s,1:number_species)
502
503 double precision, dimension(ixI^S) :: tmp1
504
505 tmp1(ixo^s)=0.5d0*(wlc(ixo^s,mom(idim))+wrc(ixo^s,mom(idim)))
506 if(present(cmin)) then
507 cmax(ixo^s,1)=max(tmp1(ixo^s)+one,zero)
508 cmin(ixo^s,1)=min(tmp1(ixo^s)-one,zero)
509 else
510 cmax(ixo^s,1)=abs(tmp1(ixo^s))+one
511 end if
512
513 end subroutine mf_get_cbounds
514
515 !> prepare velocities for ct methods
516 subroutine mf_get_ct_velocity(vcts,wLp,wRp,ixI^L,ixO^L,idim,cmax,cmin)
518
519 integer, intent(in) :: ixi^l, ixo^l, idim
520 double precision, intent(in) :: wlp(ixi^s, nw), wrp(ixi^s, nw)
521 double precision, intent(in) :: cmax(ixi^s)
522 double precision, intent(in), optional :: cmin(ixi^s)
523 type(ct_velocity), intent(inout):: vcts
524
525 integer :: idime,idimn
526
527 ! calculate velocities related to different UCT schemes
528 select case(type_ct)
529 case('average')
530 case('uct_contact')
531 if(.not.allocated(vcts%vnorm)) allocate(vcts%vnorm(ixi^s,1:ndim))
532 ! get average normal velocity at cell faces
533 vcts%vnorm(ixo^s,idim)=0.5d0*(wlp(ixo^s,mom(idim))+wrp(ixo^s,mom(idim)))
534 case('uct_hll')
535 if(.not.allocated(vcts%vbarC)) then
536 allocate(vcts%vbarC(ixi^s,1:ndir,2),vcts%vbarLC(ixi^s,1:ndir,2),vcts%vbarRC(ixi^s,1:ndir,2))
537 allocate(vcts%cbarmin(ixi^s,1:ndim),vcts%cbarmax(ixi^s,1:ndim))
538 end if
539 ! Store magnitude of characteristics
540 if(present(cmin)) then
541 vcts%cbarmin(ixo^s,idim)=max(-cmin(ixo^s),zero)
542 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
543 else
544 vcts%cbarmax(ixo^s,idim)=max( cmax(ixo^s),zero)
545 vcts%cbarmin(ixo^s,idim)=vcts%cbarmax(ixo^s,idim)
546 end if
547
548 idimn=mod(idim,ndir)+1 ! 'Next' direction
549 idime=mod(idim+1,ndir)+1 ! Electric field direction
550 ! Store velocities
551 vcts%vbarLC(ixo^s,idim,1)=wlp(ixo^s,mom(idimn))
552 vcts%vbarRC(ixo^s,idim,1)=wrp(ixo^s,mom(idimn))
553 vcts%vbarC(ixo^s,idim,1)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,1) &
554 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
555 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
556
557 vcts%vbarLC(ixo^s,idim,2)=wlp(ixo^s,mom(idime))
558 vcts%vbarRC(ixo^s,idim,2)=wrp(ixo^s,mom(idime))
559 vcts%vbarC(ixo^s,idim,2)=(vcts%cbarmax(ixo^s,idim)*vcts%vbarLC(ixo^s,idim,2) &
560 +vcts%cbarmin(ixo^s,idim)*vcts%vbarRC(ixo^s,idim,1))&
561 /(vcts%cbarmax(ixo^s,idim)+vcts%cbarmin(ixo^s,idim))
562 case default
563 call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
564 end select
565
566 end subroutine mf_get_ct_velocity
567
568 !> Calculate fluxes within ixO^L.
569 subroutine mf_get_flux(wC,w,x,ixI^L,ixO^L,idim,f)
572 use mod_geometry
573
574 integer, intent(in) :: ixi^l, ixo^l, idim
575 ! conservative w
576 double precision, intent(in) :: wc(ixi^s,nw)
577 ! primitive w
578 double precision, intent(in) :: w(ixi^s,nw)
579 double precision, intent(in) :: x(ixi^s,1:ndim)
580 double precision,intent(out) :: f(ixi^s,nwflux)
581
582 integer :: ix^d
583
584 {do ix^db=ixomin^db,ixomax^db\}
585 ! compute flux of magnetic field
586 ! f_i[b_k]=v_i*b_k-v_k*b_i
587 ^c&f(ix^d,b^c_)=w(ix^d,mom(idim))*w(ix^d,b^c_)-w(ix^d,mag(idim))*w(ix^d,m^c_)\
588 {end do\}
589 if(mf_glm) then
590 {do ix^db=ixomin^db,ixomax^db\}
591 f(ix^d,mag(idim))=w(ix^d,psi_)
592 !f_i[psi]=Ch^2*b_{i} Eq. 24e and Eq. 38c Dedner et al 2002 JCP, 175, 645
593 f(ix^d,psi_)=cmax_global**2*w(ix^d,mag(idim))
594 {end do\}
595 end if
596
597 end subroutine mf_get_flux
598
599 !> Add global source terms to update frictional velocity on complete domain
600 subroutine mf_velocity_update(qt,psa)
603 double precision, intent(in) :: qt !< Current time
604 type(state), target :: psa(max_blocks) !< Compute based on this state
605
606 integer :: iigrid,igrid
607 logical :: stagger_flag
608 logical :: firstmf=.true.
609
610 if(firstmf) then
611 ! point bc mpi datatype to partial type for velocity field
619 firstmf=.false.
620 end if
621
622 !$OMP PARALLEL DO PRIVATE(igrid)
623 do iigrid=1,igridstail_active; igrid=igrids_active(iigrid);
624 block=>psa(igrid)
625 call frictional_velocity(psa(igrid)%w,psa(igrid)%x,ixg^ll,ixm^ll)
626 end do
627 !$OMP END PARALLEL DO
628
629 ! only update mf velocity in ghost cells
630 stagger_flag=stagger_grid
637 stagger_grid=.false.
638 bcphys=.false.
639 call getbc(qt,0.d0,psa,mom(1),ndir)
640 bcphys=.true.
647 stagger_grid=stagger_flag
648
649 end subroutine mf_velocity_update
650
651 !> w[iws]=w[iws]+qdt*S[iws,wCT] where S is the source based on wCT within ixO
652 subroutine mf_add_source(qdt,dtfactor,ixI^L,ixO^L,wCT,wCTprim,w,x,qsourcesplit,active)
654
655 integer, intent(in) :: ixi^l, ixo^l
656 double precision, intent(in) :: qdt, dtfactor
657 double precision, intent(in) :: wct(ixi^s,1:nw),wctprim(ixi^s,1:nw), x(ixi^s,1:ndim)
658 double precision, intent(inout) :: w(ixi^s,1:nw)
659 logical, intent(in) :: qsourcesplit
660 logical, intent(inout) :: active
661
662 if (.not. qsourcesplit) then
663
664 ! Sources for resistivity in eqs. for e, B1, B2 and B3
665 if (abs(mf_eta)>smalldouble)then
666 active = .true.
667 call add_source_res2(qdt,ixi^l,ixo^l,wct,w,x)
668 end if
669
670 if (mf_eta_hyper>0.d0)then
671 active = .true.
672 call add_source_hyperres(qdt,ixi^l,ixo^l,wct,w,x)
673 end if
674
675 end if
676
677 {^nooned
678 if(source_split_divb .eqv. qsourcesplit) then
679 ! Sources related to div B
680 select case (type_divb)
681 case (divb_none)
682 ! Do nothing
683 case (divb_glm)
684 active = .true.
685 call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
686 case (divb_powel)
687 active = .true.
688 call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
689 case (divb_janhunen)
690 active = .true.
691 call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
692 case (divb_linde)
693 active = .true.
694 call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
695 case (divb_lindejanhunen)
696 active = .true.
697 call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
698 call add_source_janhunen(qdt,ixi^l,ixo^l,wct,w,x)
699 case (divb_lindepowel)
700 active = .true.
701 call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
702 call add_source_powel(qdt,ixi^l,ixo^l,wct,w,x)
703 case (divb_lindeglm)
704 active = .true.
705 call add_source_linde(qdt,ixi^l,ixo^l,wct,w,x)
706 call add_source_glm(qdt,ixi^l,ixo^l,wct,w,x)
707 case (divb_ct)
708 continue ! Do nothing
709 case (divb_multigrid)
710 continue ! Do nothing
711 case default
712 call mpistop('Unknown divB fix')
713 end select
714 end if
715 }
716
717 end subroutine mf_add_source
718
719 subroutine frictional_velocity(w,x,ixI^L,ixO^L)
721
722 integer, intent(in) :: ixi^l, ixo^l
723 double precision, intent(in) :: x(ixi^s,1:ndim)
724 double precision, intent(inout) :: w(ixi^s,1:nw)
725
726 double precision :: xmin(ndim),xmax(ndim)
727 double precision :: decay(ixo^s)
728 double precision :: current(ixi^s,7-2*ndir:3),tmp(ixo^s)
729 integer :: ix^d, idirmin,idir,jdir,kdir
730
731 call get_current(w,ixi^l,ixo^l,idirmin,current)
732 ! extrapolate current for the outmost layer,
733 ! because extrapolation of B at boundaries introduces artificial current
734{
735 if(block%is_physical_boundary(2*^d-1)) then
736 if(slab) then
737 ! exclude boundary 5 in 3D Cartesian
738 if(2*^d-1==5.and.ndim==3) then
739 else
740 current(ixomin^d^d%ixO^s,:)=current(ixomin^d+1^d%ixO^s,:)
741 end if
742 else
743 ! exclude boundary 1 spherical/cylindrical
744 if(2*^d-1>1) then
745 current(ixomin^d^d%ixO^s,:)=current(ixomin^d+1^d%ixO^s,:)
746 end if
747 end if
748 end if
749 if(block%is_physical_boundary(2*^d)) then
750 current(ixomax^d^d%ixO^s,:)=current(ixomax^d-1^d%ixO^s,:)
751 end if
752\}
753 w(ixo^s,mom(:))=0.d0
754 ! calculate Lorentz force
755 do idir=1,ndir; do jdir=idirmin,3; do kdir=1,ndir
756 if(lvc(idir,jdir,kdir)/=0) then
757 if(lvc(idir,jdir,kdir)==1) then
758 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))+current(ixo^s,jdir)*w(ixo^s,mag(kdir))
759 else
760 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))-current(ixo^s,jdir)*w(ixo^s,mag(kdir))
761 end if
762 end if
763 end do; end do; end do
764
765 tmp(ixo^s)=sum(w(ixo^s,mag(:))**2,dim=ndim+1)
766 ! frictional coefficient
767 where(tmp(ixo^s)/=0.d0)
768 tmp(ixo^s)=1.d0/(tmp(ixo^s)*mf_nu)
769 endwhere
770
771 do idir=1,ndir
772 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*tmp(ixo^s)
773 end do
774
775 ! decay frictional velocity near selected boundaries
776 ^d&xmin(^d)=xprobmin^d\
777 ^d&xmax(^d)=xprobmax^d\
778 decay(ixo^s)=1.d0
779 do idir=1,ndim
780 if(mf_decay_scale(2*idir-1)>0.d0) decay(ixo^s)=decay(ixo^s)*&
781 (1.d0-exp(-(x(ixo^s,idir)-xmin(idir))/mf_decay_scale(2*idir-1)))
782 if(mf_decay_scale(2*idir)>0.d0) decay(ixo^s)=decay(ixo^s)*&
783 (1.d0-exp((x(ixo^s,idir)-xmax(idir))/mf_decay_scale(2*idir)))
784 end do
785
786 ! saturate mf velocity at mf_vmax
787 tmp(ixo^s)=sqrt(sum(w(ixo^s,mom(:))**2,dim=ndim+1))/mf_vmax+1.d-12
788 tmp(ixo^s)=dtanh(tmp(ixo^s))/tmp(ixo^s)
789 do idir=1,ndir
790 w(ixo^s,mom(idir))=w(ixo^s,mom(idir))*tmp(ixo^s)*decay(ixo^s)
791 end do
792
793 end subroutine frictional_velocity
794
795 !> Add resistive source to w within ixO Uses 3 point stencil (1 neighbour) in
796 !> each direction, non-conservative. If the fourthorder precompiler flag is
797 !> set, uses fourth order central difference for the laplacian. Then the
798 !> stencil is 5 (2 neighbours).
799 subroutine add_source_res1(qdt,ixI^L,ixO^L,wCT,w,x)
802 use mod_geometry
803
804 integer, intent(in) :: ixi^l, ixo^l
805 double precision, intent(in) :: qdt
806 double precision, intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:ndim)
807 double precision, intent(inout) :: w(ixi^s,1:nw)
808
809 ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
810 double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
811 double precision :: gradeta(ixi^s,1:ndim), bf(ixi^s,1:ndir)
812 double precision :: tmp(ixi^s),tmp2(ixi^s)
813 integer :: ixa^l,idir,jdir,kdir,idirmin,idim,jxo^l,hxo^l,ix
814 integer :: lxo^l, kxo^l
815
816 ! Calculating resistive sources involve one extra layer
817 if (mf_4th_order) then
818 ixa^l=ixo^l^ladd2;
819 else
820 ixa^l=ixo^l^ladd1;
821 end if
822
823 if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
824 call mpistop("Error in add_source_res1: Non-conforming input limits")
825
826 ! Calculate current density and idirmin
827 call get_current(wct,ixi^l,ixo^l,idirmin,current)
828
829 if (mf_eta>zero)then
830 eta(ixa^s)=mf_eta
831 gradeta(ixo^s,1:ndim)=zero
832 else
833 call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
834 ! assumes that eta is not function of current?
835 do idim=1,ndim
836 call gradient(eta,ixi^l,ixo^l,idim,tmp)
837 gradeta(ixo^s,idim)=tmp(ixo^s)
838 end do
839 end if
840
841 bf(ixi^s,1:ndir)=wct(ixi^s,mag(1:ndir))
842
843 do idir=1,ndir
844 ! Put B_idir into tmp2 and eta*Laplace B_idir into tmp
845 if (mf_4th_order) then
846 tmp(ixo^s)=zero
847 tmp2(ixi^s)=bf(ixi^s,idir)
848 do idim=1,ndim
849 lxo^l=ixo^l+2*kr(idim,^d);
850 jxo^l=ixo^l+kr(idim,^d);
851 hxo^l=ixo^l-kr(idim,^d);
852 kxo^l=ixo^l-2*kr(idim,^d);
853 tmp(ixo^s)=tmp(ixo^s)+&
854 (-tmp2(lxo^s)+16.0d0*tmp2(jxo^s)-30.0d0*tmp2(ixo^s)+16.0d0*tmp2(hxo^s)-tmp2(kxo^s)) &
855 /(12.0d0 * dxlevel(idim)**2)
856 end do
857 else
858 tmp(ixo^s)=zero
859 tmp2(ixi^s)=bf(ixi^s,idir)
860 do idim=1,ndim
861 jxo^l=ixo^l+kr(idim,^d);
862 hxo^l=ixo^l-kr(idim,^d);
863 tmp(ixo^s)=tmp(ixo^s)+&
864 (tmp2(jxo^s)-2.0d0*tmp2(ixo^s)+tmp2(hxo^s))/dxlevel(idim)**2
865 end do
866 end if
867
868 ! Multiply by eta
869 tmp(ixo^s)=tmp(ixo^s)*eta(ixo^s)
870
871 ! Subtract grad(eta) x J = eps_ijk d_j eta J_k if eta is non-constant
872 if (mf_eta<zero)then
873 do jdir=1,ndim; do kdir=idirmin,3
874 if (lvc(idir,jdir,kdir)/=0)then
875 if (lvc(idir,jdir,kdir)==1)then
876 tmp(ixo^s)=tmp(ixo^s)-gradeta(ixo^s,jdir)*current(ixo^s,kdir)
877 else
878 tmp(ixo^s)=tmp(ixo^s)+gradeta(ixo^s,jdir)*current(ixo^s,kdir)
879 end if
880 end if
881 end do; end do
882 end if
883
884 ! Add sources related to eta*laplB-grad(eta) x J to B and e
885 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))+qdt*tmp(ixo^s)
886
887 end do ! idir
888
889 end subroutine add_source_res1
890
891 !> Add resistive source to w within ixO
892 !> Uses 5 point stencil (2 neighbours) in each direction, conservative
893 subroutine add_source_res2(qdt,ixI^L,ixO^L,wCT,w,x)
896 use mod_geometry
897
898 integer, intent(in) :: ixi^l, ixo^l
899 double precision, intent(in) :: qdt
900 double precision, intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:ndim)
901 double precision, intent(inout) :: w(ixi^s,1:nw)
902
903 ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
904 double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s),curlj(ixi^s,1:3)
905 double precision :: tmpvec(ixi^s,1:3)
906 integer :: ixa^l,idir,idirmin,idirmin1
907
908 ! resistive term already added as an electric field component
909 if(stagger_grid .and. ndim==ndir) return
910
911 ixa^l=ixo^l^ladd2;
912
913 if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
914 call mpistop("Error in add_source_res2: Non-conforming input limits")
915
916 ixa^l=ixo^l^ladd1;
917 ! Calculate current density within ixL: J=curl B, thus J_i=eps_ijk*d_j B_k
918 ! Determine exact value of idirmin while doing the loop.
919 call get_current(wct,ixi^l,ixa^l,idirmin,current)
920
921 if (mf_eta>zero)then
922 eta(ixa^s)=mf_eta
923 else
924 call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,current,eta)
925 end if
926
927 ! dB/dt= -curl(J*eta), thus B_i=B_i-eps_ijk d_j Jeta_k
928 tmpvec(ixa^s,1:ndir)=zero
929 do idir=idirmin,3
930 tmpvec(ixa^s,idir)=current(ixa^s,idir)*eta(ixa^s)
931 end do
932 curlj=0.d0
933 call curlvector(tmpvec,ixi^l,ixo^l,curlj,idirmin1,1,3)
934 if(stagger_grid) then
935 if(ndim==2.and.ndir==3) then
936 ! if 2.5D
937 w(ixo^s,mag(ndir)) = w(ixo^s,mag(ndir))-qdt*curlj(ixo^s,ndir)
938 end if
939 else
940 w(ixo^s,mag(1:ndir)) = w(ixo^s,mag(1:ndir))-qdt*curlj(ixo^s,1:ndir)
941 end if
942
943 end subroutine add_source_res2
944
945 !> Add Hyper-resistive source to w within ixO
946 !> Uses 9 point stencil (4 neighbours) in each direction.
947 subroutine add_source_hyperres(qdt,ixI^L,ixO^L,wCT,w,x)
949 use mod_geometry
950
951 integer, intent(in) :: ixi^l, ixo^l
952 double precision, intent(in) :: qdt
953 double precision, intent(in) :: wct(ixi^s,1:nw), x(ixi^s,1:ndim)
954 double precision, intent(inout) :: w(ixi^s,1:nw)
955 !.. local ..
956 double precision :: current(ixi^s,7-2*ndir:3)
957 double precision :: tmpvec(ixi^s,1:3),tmpvec2(ixi^s,1:3),tmp(ixi^s),ehyper(ixi^s,1:3)
958 integer :: ixa^l,idir,jdir,kdir,idirmin,idirmin1
959
960 ixa^l=ixo^l^ladd3;
961 if (iximin^d>ixamin^d.or.iximax^d<ixamax^d|.or.) &
962 call mpistop("Error in add_source_hyperres: Non-conforming input limits")
963
964 call get_current(wct,ixi^l,ixa^l,idirmin,current)
965 tmpvec(ixa^s,1:ndir)=zero
966 do jdir=idirmin,3
967 tmpvec(ixa^s,jdir)=current(ixa^s,jdir)
968 end do
969
970 ixa^l=ixo^l^ladd2;
971 call curlvector(tmpvec,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
972
973 ixa^l=ixo^l^ladd1;
974 tmpvec(ixa^s,1:ndir)=zero
975 call curlvector(tmpvec2,ixi^l,ixa^l,tmpvec,idirmin1,1,3)
976 ehyper(ixa^s,1:ndir) = - tmpvec(ixa^s,1:ndir)*mf_eta_hyper
977
978 ixa^l=ixo^l;
979 tmpvec2(ixa^s,1:ndir)=zero
980 call curlvector(ehyper,ixi^l,ixa^l,tmpvec2,idirmin1,1,3)
981
982 do idir=1,ndir
983 w(ixo^s,mag(idir)) = w(ixo^s,mag(idir))-tmpvec2(ixo^s,idir)*qdt
984 end do
985
986 end subroutine add_source_hyperres
987
988 subroutine add_source_glm(qdt,ixI^L,ixO^L,wCT,w,x)
989 ! Add divB related sources to w within ixO
990 ! corresponding to Dedner JCP 2002, 175, 645 _equation 24_
991 ! giving the EGLM-MHD scheme
993 use mod_geometry
994
995 integer, intent(in) :: ixi^l, ixo^l
996 double precision, intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:ndim)
997 double precision, intent(inout) :: w(ixi^s,1:nw)
998 double precision:: divb(ixi^s)
999 double precision :: gradpsi(ixi^s)
1000 integer :: idim,idir
1001
1002 ! We calculate now div B
1003 call get_divb(wct,ixi^l,ixo^l,divb,mf_divb_nth)
1004
1005 ! dPsi/dt = - Ch^2/Cp^2 Psi
1006 if (mf_glm_alpha < zero) then
1007 w(ixo^s,psi_) = abs(mf_glm_alpha)*wct(ixo^s,psi_)
1008 else
1009 ! implicit update of Psi variable
1010 ! equation (27) in Mignone 2010 J. Com. Phys. 229, 2117
1011 if(slab_uniform) then
1012 w(ixo^s,psi_) = dexp(-qdt*cmax_global*mf_glm_alpha/minval(dxlevel(:)))*w(ixo^s,psi_)
1013 else
1014 w(ixo^s,psi_) = dexp(-qdt*cmax_global*mf_glm_alpha/minval(block%ds(ixo^s,:),dim=ndim+1))*w(ixo^s,psi_)
1015 end if
1016 end if
1017
1018 ! gradient of Psi
1019 do idim=1,ndim
1020 select case(typegrad)
1021 case("central")
1022 call gradient(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
1023 case("limited")
1024 call gradientl(wct(ixi^s,psi_),ixi^l,ixo^l,idim,gradpsi)
1025 end select
1026 end do
1027 end subroutine add_source_glm
1028
1029 !> Add divB related sources to w within ixO corresponding to Powel
1030 subroutine add_source_powel(qdt,ixI^L,ixO^L,wCT,w,x)
1032
1033 integer, intent(in) :: ixi^l, ixo^l
1034 double precision, intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:ndim)
1035 double precision, intent(inout) :: w(ixi^s,1:nw)
1036
1037 double precision :: divb(ixi^s)
1038 integer :: idir
1039
1040 ! We calculate now div B
1041 call get_divb(wct,ixi^l,ixo^l,divb,mf_divb_nth)
1042
1043 ! b = b - qdt v * div b
1044 do idir=1,ndir
1045 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,mom(idir))*divb(ixo^s)
1046 end do
1047
1048 end subroutine add_source_powel
1049
1050 subroutine add_source_janhunen(qdt,ixI^L,ixO^L,wCT,w,x)
1051 ! Add divB related sources to w within ixO
1052 ! corresponding to Janhunen, just the term in the induction equation.
1054
1055 integer, intent(in) :: ixi^l, ixo^l
1056 double precision, intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:ndim)
1057 double precision, intent(inout) :: w(ixi^s,1:nw)
1058 double precision :: divb(ixi^s)
1059 integer :: idir
1060
1061 ! We calculate now div B
1062 call get_divb(wct,ixi^l,ixo^l,divb,mf_divb_nth)
1063
1064 ! b = b - qdt v * div b
1065 do idir=1,ndir
1066 w(ixo^s,mag(idir))=w(ixo^s,mag(idir))-qdt*wct(ixo^s,mom(idir))*divb(ixo^s)
1067 end do
1068
1069 end subroutine add_source_janhunen
1070
1071 subroutine add_source_linde(qdt,ixI^L,ixO^L,wCT,w,x)
1072 ! Add Linde's divB related sources to wnew within ixO
1074 use mod_geometry
1075
1076 integer, intent(in) :: ixi^l, ixo^l
1077 double precision, intent(in) :: qdt, wct(ixi^s,1:nw), x(ixi^s,1:ndim)
1078 double precision, intent(inout) :: w(ixi^s,1:nw)
1079 double precision :: divb(ixi^s),graddivb(ixi^s)
1080 integer :: idim, idir, ixp^l, i^d, iside
1081 logical, dimension(-1:1^D&) :: leveljump
1082
1083 ! Calculate div B
1084 ixp^l=ixo^l^ladd1;
1085 call get_divb(wct,ixi^l,ixp^l,divb,mf_divb_nth)
1086
1087 ! for AMR stability, retreat one cell layer from the boarders of level jump
1088 {do i^db=-1,1\}
1089 if(i^d==0|.and.) cycle
1090 if(neighbor_type(i^d,block%igrid)==2 .or. neighbor_type(i^d,block%igrid)==4) then
1091 leveljump(i^d)=.true.
1092 else
1093 leveljump(i^d)=.false.
1094 end if
1095 {end do\}
1096
1097 ixp^l=ixo^l;
1098 do idim=1,ndim
1099 select case(idim)
1100 {case(^d)
1101 do iside=1,2
1102 i^dd=kr(^dd,^d)*(2*iside-3);
1103 if (leveljump(i^dd)) then
1104 if (iside==1) then
1105 ixpmin^d=ixomin^d-i^d
1106 else
1107 ixpmax^d=ixomax^d-i^d
1108 end if
1109 end if
1110 end do
1111 \}
1112 end select
1113 end do
1114
1115 ! Add Linde's diffusive terms
1116 do idim=1,ndim
1117 ! Calculate grad_idim(divb)
1118 select case(typegrad)
1119 case("central")
1120 call gradient(divb,ixi^l,ixp^l,idim,graddivb)
1121 case("limited")
1122 call gradientl(divb,ixi^l,ixp^l,idim,graddivb)
1123 end select
1124
1125 ! Multiply by Linde's eta*dt = divbdiff*(c_max*dx)*dt = divbdiff*dx**2
1126 if (slab_uniform) then
1127 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff/(^d&1.0d0/dxlevel(^d)**2+)
1128 else
1129 graddivb(ixp^s)=graddivb(ixp^s)*divbdiff &
1130 /(^d&1.0d0/block%ds(ixp^s,^d)**2+)
1131 end if
1132
1133 w(ixp^s,mag(idim))=w(ixp^s,mag(idim))+graddivb(ixp^s)
1134 end do
1135
1136 end subroutine add_source_linde
1137
1138
1139 !> get dimensionless div B = |divB| * volume / area / |B|
1140 subroutine get_normalized_divb(w,ixI^L,ixO^L,divb)
1141
1143
1144 integer, intent(in) :: ixi^l, ixo^l
1145 double precision, intent(in) :: w(ixi^s,1:nw)
1146 double precision :: divb(ixi^s), dsurface(ixi^s)
1147
1148 double precision :: invb(ixo^s)
1149 integer :: ixa^l,idims
1150
1151 call get_divb(w,ixi^l,ixo^l,divb)
1152 invb(ixo^s)=sqrt(sum(w(ixo^s, mag(:))**2, dim=ndim+1))
1153 where(invb(ixo^s)/=0.d0)
1154 invb(ixo^s)=1.d0/invb(ixo^s)
1155 end where
1156 if(slab_uniform) then
1157 divb(ixo^s)=0.5d0*abs(divb(ixo^s))*invb(ixo^s)/sum(1.d0/dxlevel(:))
1158 else
1159 ixamin^d=ixomin^d-1;
1160 ixamax^d=ixomax^d-1;
1161 dsurface(ixo^s)= sum(block%surfaceC(ixo^s,:),dim=ndim+1)
1162 do idims=1,ndim
1163 ixa^l=ixo^l-kr(idims,^d);
1164 dsurface(ixo^s)=dsurface(ixo^s)+block%surfaceC(ixa^s,idims)
1165 end do
1166 divb(ixo^s)=abs(divb(ixo^s))*invb(ixo^s)*&
1167 block%dvolume(ixo^s)/dsurface(ixo^s)
1168 end if
1169
1170 end subroutine get_normalized_divb
1171
1172 !> Calculate idirmin and the idirmin:3 components of the common current array
1173 !> make sure that dxlevel(^D) is set correctly.
1174 subroutine get_current(w,ixI^L,ixO^L,idirmin,current,fourthorder)
1176 use mod_geometry
1177
1178 integer, intent(in) :: ixo^l, ixi^l
1179 double precision, intent(in) :: w(ixi^s,1:nw)
1180 integer, intent(out) :: idirmin
1181 logical, intent(in), optional :: fourthorder
1182
1183 ! For ndir=2 only 3rd component of J can exist, ndir=1 is impossible for MHD
1184 double precision :: current(ixi^s,7-2*ndir:3)
1185 integer :: idir, idirmin0
1186
1187 idirmin0 = 7-2*ndir
1188
1189 if(present(fourthorder)) then
1190 call curlvector(w(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir,fourthorder)
1191 else
1192 call curlvector(w(ixi^s,mag(1:ndir)),ixi^l,ixo^l,current,idirmin,idirmin0,ndir)
1193 end if
1194
1195 end subroutine get_current
1196
1197 !> If resistivity is not zero, check diffusion time limit for dt
1198 subroutine mf_get_dt(w,ixI^L,ixO^L,dtnew,dx^D,x)
1200 use mod_usr_methods
1201
1202 integer, intent(in) :: ixi^l, ixo^l
1203 double precision, intent(inout) :: dtnew
1204 double precision, intent(in) :: dx^d
1205 double precision, intent(in) :: w(ixi^s,1:nw)
1206 double precision, intent(in) :: x(ixi^s,1:ndim)
1207
1208 double precision :: dxarr(ndim)
1209 double precision :: current(ixi^s,7-2*ndir:3),eta(ixi^s)
1210 integer :: idirmin,idim
1211
1212 dtnew = bigdouble
1213
1214 ^d&dxarr(^d)=dx^d;
1215 if (mf_eta>zero)then
1216 dtnew=dtdiffpar*minval(dxarr(1:ndim))**2/mf_eta
1217 else if (mf_eta<zero)then
1218 call get_current(w,ixi^l,ixo^l,idirmin,current)
1219 call usr_special_resistivity(w,ixi^l,ixo^l,idirmin,x,current,eta)
1220 dtnew=bigdouble
1221 do idim=1,ndim
1222 if(slab_uniform) then
1223 dtnew=min(dtnew,&
1224 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/dxarr(idim)**2)))
1225 else
1226 dtnew=min(dtnew,&
1227 dtdiffpar/(smalldouble+maxval(eta(ixo^s)/block%ds(ixo^s,idim)**2)))
1228 end if
1229 end do
1230 end if
1231
1232 if(mf_eta_hyper>zero) then
1233 if(slab_uniform) then
1234 dtnew=min(dtdiffpar*minval(dxarr(1:ndim))**4/mf_eta_hyper,dtnew)
1235 else
1236 dtnew=min(dtdiffpar*minval(block%ds(ixo^s,1:ndim))**4/mf_eta_hyper,dtnew)
1237 end if
1238 end if
1239 end subroutine mf_get_dt
1240
1241 ! Add geometrical source terms to w
1242 subroutine mf_add_source_geom(qdt,dtfactor, ixI^L,ixO^L,wCT,wprim,w,x)
1244 use mod_geometry
1245
1246 integer, intent(in) :: ixi^l, ixo^l
1247 double precision, intent(in) :: qdt, dtfactor, x(ixi^s,1:ndim)
1248 double precision, intent(inout) :: wct(ixi^s,1:nw), wprim(ixi^s,1:nw), w(ixi^s,1:nw)
1249
1250 double precision :: tmp,invr,cs2
1251 integer :: iw,idir,ix^d
1252
1253 integer :: mr_,mphi_ ! Polar var. names
1254 integer :: br_,bphi_
1255
1256 mr_=mom(1); mphi_=mom(1)-1+phi_ ! Polar var. names
1257 br_=mag(1); bphi_=mag(1)-1+phi_
1258
1259 select case (coordinate)
1260 case (cylindrical)
1261 if(phi_>0) then
1262 if(.not.stagger_grid) then
1263 {do ix^db=ixomin^db,ixomax^db\}
1264 w(ix^d,bphi_)=w(ix^d,bphi_)+qdt/x(ix^d,1)*&
1265 (wct(ix^d,bphi_)*wct(ix^d,mr_) &
1266 -wct(ix^d,br_)*wct(ix^d,mphi_))
1267 {end do\}
1268 end if
1269 end if
1270 if(mf_glm) w(ixo^s,br_)=w(ixo^s,br_)+qdt*wct(ixo^s,psi_)/x(ixo^s,1)
1271 case (spherical)
1272 {do ix^db=ixomin^db,ixomax^db\}
1273 invr=qdt/x(ix^d,1)
1274 ! b1
1275 if(mf_glm) then
1276 w(ix^d,mag(1))=w(ix^d,mag(1))+invr*2.0d0*wct(ix^d,psi_)
1277 end if
1278
1279 {^nooned
1280 ! b2
1281 if(.not.stagger_grid) then
1282 tmp=wct(ix^d,mom(1))*wct(ix^d,mag(2))-wct(ix^d,mom(2))*wct(ix^d,mag(1))
1283 cs2=1.d0/tan(x(ix^d,2))
1284 if(mf_glm) then
1285 tmp=tmp+cs2*wct(ix^d,psi_)
1286 end if
1287 w(ix^d,mag(2))=w(ix^d,mag(2))+tmp*invr
1288 end if
1289 }
1290
1291 if(ndir==3) then
1292 {^ifoned
1293 ! b3
1294 if(.not.stagger_grid) then
1295 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
1296 (wct(ix^d,mom(1))*wct(ix^d,mag(3)) &
1297 -wct(ix^d,mom(3))*wct(ix^d,mag(1)))
1298 end if
1299 }
1300 {^nooned
1301 ! b3
1302 if(.not.stagger_grid) then
1303 w(ix^d,mag(3))=w(ix^d,mag(3))+invr*&
1304 (wct(ix^d,mom(1))*wct(ix^d,mag(3)) &
1305 -wct(ix^d,mom(3))*wct(ix^d,mag(1)) &
1306 -(wct(ix^d,mom(3))*wct(ix^d,mag(2)) &
1307 -wct(ix^d,mom(2))*wct(ix^d,mag(3)))*cs2)
1308 end if
1309 }
1310 end if
1311 {end do\}
1312 end select
1313
1314 end subroutine mf_add_source_geom
1315
1316 subroutine mf_modify_wlr(ixI^L,ixO^L,qt,wLC,wRC,wLp,wRp,s,idir)
1318 use mod_usr_methods
1319 integer, intent(in) :: ixi^l, ixo^l, idir
1320 double precision, intent(in) :: qt
1321 double precision, intent(inout) :: wlc(ixi^s,1:nw), wrc(ixi^s,1:nw)
1322 double precision, intent(inout) :: wlp(ixi^s,1:nw), wrp(ixi^s,1:nw)
1323 type(state) :: s
1324 double precision :: db(ixi^s), dpsi(ixi^s)
1325
1326 ! Solve the Riemann problem for the linear 2x2 system for normal
1327 ! B-field and GLM_Psi according to Dedner 2002:
1328 ! This implements eq. (42) in Dedner et al. 2002 JcP 175
1329 ! Gives the Riemann solution on the interface
1330 ! for the normal B component and Psi in the GLM-MHD system.
1331 ! 23/04/2013 Oliver Porth
1332 db(ixo^s) = wrp(ixo^s,mag(idir)) - wlp(ixo^s,mag(idir))
1333 dpsi(ixo^s) = wrp(ixo^s,psi_) - wlp(ixo^s,psi_)
1334
1335 wlp(ixo^s,mag(idir)) = 0.5d0 * (wrp(ixo^s,mag(idir)) + wlp(ixo^s,mag(idir))) &
1336 - 0.5d0/cmax_global * dpsi(ixo^s)
1337 wlp(ixo^s,psi_) = 0.5d0 * (wrp(ixo^s,psi_) + wlp(ixo^s,psi_)) &
1338 - 0.5d0*cmax_global * db(ixo^s)
1339
1340 wrp(ixo^s,mag(idir)) = wlp(ixo^s,mag(idir))
1341 wrp(ixo^s,psi_) = wlp(ixo^s,psi_)
1342
1343 end subroutine mf_modify_wlr
1344
1345 subroutine mf_boundary_adjust(igrid,psb)
1347 integer, intent(in) :: igrid
1348 type(state), target :: psb(max_blocks)
1349
1350 integer :: ib, idims, iside, ixo^l, i^d
1351
1352 block=>ps(igrid)
1353 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
1354 do idims=1,ndim
1355 ! to avoid using as yet unknown corner info in more than 1D, we
1356 ! fill only interior mesh ranges of the ghost cell ranges at first,
1357 ! and progressively enlarge the ranges to include corners later
1358 do iside=1,2
1359 i^d=kr(^d,idims)*(2*iside-3);
1360 if (neighbor_type(i^d,igrid)/=1) cycle
1361 ib=(idims-1)*2+iside
1362 if(.not.boundary_divbfix(ib)) cycle
1363 if(any(typeboundary(:,ib)==bc_special)) then
1364 ! MF nonlinear force-free B field extrapolation and data driven
1365 ! require normal B of the first ghost cell layer to be untouched by
1366 ! fixdivB=0 process, set boundary_divbfix_skip(iB)=1 in par file
1367 select case (idims)
1368 {case (^d)
1369 if (iside==2) then
1370 ! maximal boundary
1371 ixomin^dd=ixghi^d+1-nghostcells+boundary_divbfix_skip(2*^d)^d%ixOmin^dd=ixglo^dd;
1372 ixomax^dd=ixghi^dd;
1373 else
1374 ! minimal boundary
1375 ixomin^dd=ixglo^dd;
1376 ixomax^dd=ixglo^d-1+nghostcells-boundary_divbfix_skip(2*^d-1)^d%ixOmax^dd=ixghi^dd;
1377 end if \}
1378 end select
1379 call fixdivb_boundary(ixg^ll,ixo^l,psb(igrid)%w,psb(igrid)%x,ib)
1380 end if
1381 end do
1382 end do
1383
1384 end subroutine mf_boundary_adjust
1385
1386 subroutine fixdivb_boundary(ixG^L,ixO^L,w,x,iB)
1388
1389 integer, intent(in) :: ixg^l,ixo^l,ib
1390 double precision, intent(inout) :: w(ixg^s,1:nw)
1391 double precision, intent(in) :: x(ixg^s,1:ndim)
1392
1393 double precision :: dx1x2,dx1x3,dx2x1,dx2x3,dx3x1,dx3x2
1394 integer :: ix^d,ixf^l
1395
1396 select case(ib)
1397 case(1)
1398 ! 2nd order CD for divB=0 to set normal B component better
1399 {^iftwod
1400 ixfmin1=ixomin1+1
1401 ixfmax1=ixomax1+1
1402 ixfmin2=ixomin2+1
1403 ixfmax2=ixomax2-1
1404 if(slab_uniform) then
1405 dx1x2=dxlevel(1)/dxlevel(2)
1406 do ix1=ixfmax1,ixfmin1,-1
1407 w(ix1-1,ixfmin2:ixfmax2,mag(1))=w(ix1+1,ixfmin2:ixfmax2,mag(1)) &
1408 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
1409 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
1410 enddo
1411 else
1412 do ix1=ixfmax1,ixfmin1,-1
1413 w(ix1-1,ixfmin2:ixfmax2,mag(1))=( (w(ix1+1,ixfmin2:ixfmax2,mag(1))+&
1414 w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1,ixfmin2:ixfmax2,1)&
1415 +(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
1416 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
1417 -(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
1418 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
1419 /block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
1420 end do
1421 end if
1422 }
1423 {^ifthreed
1424 ixfmin1=ixomin1+1
1425 ixfmax1=ixomax1+1
1426 ixfmin2=ixomin2+1
1427 ixfmax2=ixomax2-1
1428 ixfmin3=ixomin3+1
1429 ixfmax3=ixomax3-1
1430 if(slab_uniform) then
1431 dx1x2=dxlevel(1)/dxlevel(2)
1432 dx1x3=dxlevel(1)/dxlevel(3)
1433 do ix1=ixfmax1,ixfmin1,-1
1434 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1435 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
1436 +dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
1437 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
1438 +dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
1439 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
1440 end do
1441 else
1442 do ix1=ixfmax1,ixfmin1,-1
1443 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1444 ( (w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
1445 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
1446 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
1447 +(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
1448 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
1449 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
1450 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
1451 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
1452 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
1453 +(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
1454 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
1455 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
1456 -(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
1457 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1458 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
1459 /block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
1460 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
1461 end do
1462 end if
1463 }
1464 case(2)
1465 {^iftwod
1466 ixfmin1=ixomin1-1
1467 ixfmax1=ixomax1-1
1468 ixfmin2=ixomin2+1
1469 ixfmax2=ixomax2-1
1470 if(slab_uniform) then
1471 dx1x2=dxlevel(1)/dxlevel(2)
1472 do ix1=ixfmin1,ixfmax1
1473 w(ix1+1,ixfmin2:ixfmax2,mag(1))=w(ix1-1,ixfmin2:ixfmax2,mag(1)) &
1474 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))-&
1475 w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))
1476 enddo
1477 else
1478 do ix1=ixfmin1,ixfmax1
1479 w(ix1+1,ixfmin2:ixfmax2,mag(1))=( (w(ix1-1,ixfmin2:ixfmax2,mag(1))+&
1480 w(ix1,ixfmin2:ixfmax2,mag(1)))*block%surfaceC(ix1-1,ixfmin2:ixfmax2,1)&
1481 -(w(ix1,ixfmin2+1:ixfmax2+1,mag(2))+w(ix1,ixfmin2:ixfmax2,mag(2)))*&
1482 block%surfaceC(ix1,ixfmin2:ixfmax2,2)&
1483 +(w(ix1,ixfmin2:ixfmax2,mag(2))+w(ix1,ixfmin2-1:ixfmax2-1,mag(2)))*&
1484 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,2) )&
1485 /block%surfaceC(ix1,ixfmin2:ixfmax2,1)-w(ix1,ixfmin2:ixfmax2,mag(1))
1486 end do
1487 end if
1488 }
1489 {^ifthreed
1490 ixfmin1=ixomin1-1
1491 ixfmax1=ixomax1-1
1492 ixfmin2=ixomin2+1
1493 ixfmax2=ixomax2-1
1494 ixfmin3=ixomin3+1
1495 ixfmax3=ixomax3-1
1496 if(slab_uniform) then
1497 dx1x2=dxlevel(1)/dxlevel(2)
1498 dx1x3=dxlevel(1)/dxlevel(3)
1499 do ix1=ixfmin1,ixfmax1
1500 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1501 w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)) &
1502 -dx1x2*(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))-&
1503 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2))) &
1504 -dx1x3*(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))-&
1505 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))
1506 end do
1507 else
1508 do ix1=ixfmin1,ixfmax1
1509 w(ix1+1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))=&
1510 ( (w(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))+&
1511 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1)))*&
1512 block%surfaceC(ix1-1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)&
1513 -(w(ix1,ixfmin2+1:ixfmax2+1,ixfmin3:ixfmax3,mag(2))+&
1514 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2)))*&
1515 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,2)&
1516 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(2))+&
1517 w(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,mag(2)))*&
1518 block%surfaceC(ix1,ixfmin2-1:ixfmax2-1,ixfmin3:ixfmax3,2)&
1519 -(w(ix1,ixfmin2:ixfmax2,ixfmin3+1:ixfmax3+1,mag(3))+&
1520 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3)))*&
1521 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,3)&
1522 +(w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(3))+&
1523 w(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1524 block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3-1:ixfmax3-1,3) )&
1525 /block%surfaceC(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,1)-&
1526 w(ix1,ixfmin2:ixfmax2,ixfmin3:ixfmax3,mag(1))
1527 end do
1528 end if
1529 }
1530 case(3)
1531 {^iftwod
1532 ixfmin1=ixomin1+1
1533 ixfmax1=ixomax1-1
1534 ixfmin2=ixomin2+1
1535 ixfmax2=ixomax2+1
1536 if(slab_uniform) then
1537 dx2x1=dxlevel(2)/dxlevel(1)
1538 do ix2=ixfmax2,ixfmin2,-1
1539 w(ixfmin1:ixfmax1,ix2-1,mag(2))=w(ixfmin1:ixfmax1,ix2+1,mag(2)) &
1540 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
1541 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
1542 enddo
1543 else
1544 do ix2=ixfmax2,ixfmin2,-1
1545 w(ixfmin1:ixfmax1,ix2-1,mag(2))=( (w(ixfmin1:ixfmax1,ix2+1,mag(2))+&
1546 w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2,2)&
1547 +(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
1548 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
1549 -(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
1550 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
1551 /block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
1552 end do
1553 end if
1554 }
1555 {^ifthreed
1556 ixfmin1=ixomin1+1
1557 ixfmax1=ixomax1-1
1558 ixfmin3=ixomin3+1
1559 ixfmax3=ixomax3-1
1560 ixfmin2=ixomin2+1
1561 ixfmax2=ixomax2+1
1562 if(slab_uniform) then
1563 dx2x1=dxlevel(2)/dxlevel(1)
1564 dx2x3=dxlevel(2)/dxlevel(3)
1565 do ix2=ixfmax2,ixfmin2,-1
1566 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
1567 ix2+1,ixfmin3:ixfmax3,mag(2)) &
1568 +dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
1569 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
1570 +dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
1571 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
1572 end do
1573 else
1574 do ix2=ixfmax2,ixfmin2,-1
1575 w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))=&
1576 ( (w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))+&
1577 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
1578 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)&
1579 +(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
1580 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1581 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
1582 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
1583 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1584 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
1585 +(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
1586 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
1587 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
1588 -(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
1589 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1590 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
1591 /block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)-&
1592 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
1593 end do
1594 end if
1595 }
1596 case(4)
1597 {^iftwod
1598 ixfmin1=ixomin1+1
1599 ixfmax1=ixomax1-1
1600 ixfmin2=ixomin2-1
1601 ixfmax2=ixomax2-1
1602 if(slab_uniform) then
1603 dx2x1=dxlevel(2)/dxlevel(1)
1604 do ix2=ixfmin2,ixfmax2
1605 w(ixfmin1:ixfmax1,ix2+1,mag(2))=w(ixfmin1:ixfmax1,ix2-1,mag(2)) &
1606 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))-&
1607 w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))
1608 end do
1609 else
1610 do ix2=ixfmin2,ixfmax2
1611 w(ixfmin1:ixfmax1,ix2+1,mag(2))=( (w(ixfmin1:ixfmax1,ix2-1,mag(2))+&
1612 w(ixfmin1:ixfmax1,ix2,mag(2)))*block%surfaceC(ixfmin1:ixfmax1,ix2-1,2)&
1613 -(w(ixfmin1+1:ixfmax1+1,ix2,mag(1))+w(ixfmin1:ixfmax1,ix2,mag(1)))*&
1614 block%surfaceC(ixfmin1:ixfmax1,ix2,1)&
1615 +(w(ixfmin1:ixfmax1,ix2,mag(1))+w(ixfmin1-1:ixfmax1-1,ix2,mag(1)))*&
1616 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,1) )&
1617 /block%surfaceC(ixfmin1:ixfmax1,ix2,2)-w(ixfmin1:ixfmax1,ix2,mag(2))
1618 end do
1619 end if
1620 }
1621 {^ifthreed
1622 ixfmin1=ixomin1+1
1623 ixfmax1=ixomax1-1
1624 ixfmin3=ixomin3+1
1625 ixfmax3=ixomax3-1
1626 ixfmin2=ixomin2-1
1627 ixfmax2=ixomax2-1
1628 if(slab_uniform) then
1629 dx2x1=dxlevel(2)/dxlevel(1)
1630 dx2x3=dxlevel(2)/dxlevel(3)
1631 do ix2=ixfmin2,ixfmax2
1632 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=w(ixfmin1:ixfmax1,&
1633 ix2-1,ixfmin3:ixfmax3,mag(2)) &
1634 -dx2x1*(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))-&
1635 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1))) &
1636 -dx2x3*(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))-&
1637 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))
1638 end do
1639 else
1640 do ix2=ixfmin2,ixfmax2
1641 w(ixfmin1:ixfmax1,ix2+1,ixfmin3:ixfmax3,mag(2))=&
1642 ( (w(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,mag(2))+&
1643 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2)))*&
1644 block%surfaceC(ixfmin1:ixfmax1,ix2-1,ixfmin3:ixfmax3,2)&
1645 -(w(ixfmin1+1:ixfmax1+1,ix2,ixfmin3:ixfmax3,mag(1))+&
1646 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1647 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,1)&
1648 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(1))+&
1649 w(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,mag(1)))*&
1650 block%surfaceC(ixfmin1-1:ixfmax1-1,ix2,ixfmin3:ixfmax3,1)&
1651 -(w(ixfmin1:ixfmax1,ix2,ixfmin3+1:ixfmax3+1,mag(3))+&
1652 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3)))*&
1653 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,3)&
1654 +(w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(3))+&
1655 w(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,mag(3)))*&
1656 block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3-1:ixfmax3-1,3) )&
1657 /block%surfaceC(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,2)-&
1658 w(ixfmin1:ixfmax1,ix2,ixfmin3:ixfmax3,mag(2))
1659 end do
1660 end if
1661 }
1662 {^ifthreed
1663 case(5)
1664 ixfmin1=ixomin1+1
1665 ixfmax1=ixomax1-1
1666 ixfmin2=ixomin2+1
1667 ixfmax2=ixomax2-1
1668 ixfmin3=ixomin3+1
1669 ixfmax3=ixomax3+1
1670 if(slab_uniform) then
1671 dx3x1=dxlevel(3)/dxlevel(1)
1672 dx3x2=dxlevel(3)/dxlevel(2)
1673 do ix3=ixfmax3,ixfmin3,-1
1674 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=w(ixfmin1:ixfmax1,&
1675 ixfmin2:ixfmax2,ix3+1,mag(3)) &
1676 +dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
1677 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
1678 +dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
1679 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
1680 end do
1681 else
1682 do ix3=ixfmax3,ixfmin3,-1
1683 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))=&
1684 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))+&
1685 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
1686 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)&
1687 +(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
1688 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1689 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
1690 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
1691 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1692 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
1693 +(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
1694 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
1695 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
1696 -(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
1697 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
1698 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
1699 /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)-&
1700 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
1701 end do
1702 end if
1703 case(6)
1704 ixfmin1=ixomin1+1
1705 ixfmax1=ixomax1-1
1706 ixfmin2=ixomin2+1
1707 ixfmax2=ixomax2-1
1708 ixfmin3=ixomin3-1
1709 ixfmax3=ixomax3-1
1710 if(slab_uniform) then
1711 dx3x1=dxlevel(3)/dxlevel(1)
1712 dx3x2=dxlevel(3)/dxlevel(2)
1713 do ix3=ixfmin3,ixfmax3
1714 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=w(ixfmin1:ixfmax1,&
1715 ixfmin2:ixfmax2,ix3-1,mag(3)) &
1716 -dx3x1*(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))-&
1717 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1))) &
1718 -dx3x2*(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))-&
1719 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))
1720 end do
1721 else
1722 do ix3=ixfmin3,ixfmax3
1723 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3+1,mag(3))=&
1724 ( (w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,mag(3))+&
1725 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3)))*&
1726 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3-1,3)&
1727 -(w(ixfmin1+1:ixfmax1+1,ixfmin2:ixfmax2,ix3,mag(1))+&
1728 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1729 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,1)&
1730 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(1))+&
1731 w(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,mag(1)))*&
1732 block%surfaceC(ixfmin1-1:ixfmax1-1,ixfmin2:ixfmax2,ix3,1)&
1733 -(w(ixfmin1:ixfmax1,ixfmin2+1:ixfmax2+1,ix3,mag(2))+&
1734 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2)))*&
1735 block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,2)&
1736 +(w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(2))+&
1737 w(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,mag(2)))*&
1738 block%surfaceC(ixfmin1:ixfmax1,ixfmin2-1:ixfmax2-1,ix3,2) )&
1739 /block%surfaceC(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,3)-&
1740 w(ixfmin1:ixfmax1,ixfmin2:ixfmax2,ix3,mag(3))
1741 end do
1742 end if
1743 }
1744 case default
1745 call mpistop("Special boundary is not defined for this region")
1746 end select
1747
1748 end subroutine fixdivb_boundary
1749
1750 {^nooned
1751 subroutine mf_clean_divb_multigrid(qdt, qt, active)
1752 use mod_forest
1755 use mod_geometry
1756
1757 double precision, intent(in) :: qdt !< Current time step
1758 double precision, intent(in) :: qt !< Current time
1759 logical, intent(inout) :: active !< Output if the source is active
1760
1761 double precision :: tmp(ixg^t), grad(ixg^t, ndim)
1762 double precision :: res
1763 double precision, parameter :: max_residual = 1d-3
1764 double precision, parameter :: residual_reduction = 1d-10
1765 integer :: id
1766 integer, parameter :: max_its = 50
1767 double precision :: residual_it(max_its), max_divb
1768 integer :: iigrid, igrid
1769 integer :: n, nc, lvl, ix^l, ixc^l, idim
1770 type(tree_node), pointer :: pnode
1771
1772 mg%operator_type = mg_laplacian
1773
1774 ! Set boundary conditions
1775 do n = 1, 2*ndim
1776 idim = (n+1)/2
1777 select case (typeboundary(mag(idim), n))
1778 case (bc_symm)
1779 ! d/dx B = 0, take phi = 0
1780 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1781 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1782 case (bc_asymm)
1783 ! B = 0, so grad(phi) = 0
1784 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
1785 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1786 case (bc_cont)
1787 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1788 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1789 case (bc_special)
1790 ! Assume Dirichlet boundary conditions, derivative zero
1791 mg%bc(n, mg_iphi)%bc_type = mg_bc_neumann
1792 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1793 case (bc_periodic)
1794 ! Nothing to do here
1795 case default
1796 write(*,*) "mf_clean_divb_multigrid warning: unknown boundary type"
1797 mg%bc(n, mg_iphi)%bc_type = mg_bc_dirichlet
1798 mg%bc(n, mg_iphi)%bc_value = 0.0_dp
1799 end select
1800 end do
1801
1802 ix^l=ixm^ll^ladd1;
1803 max_divb = 0.0d0
1804
1805 ! Store divergence of B as right-hand side
1806 do iigrid = 1, igridstail
1807 igrid = igrids(iigrid);
1808 pnode => igrid_to_node(igrid, mype)%node
1809 id = pnode%id
1810 lvl = mg%boxes(id)%lvl
1811 nc = mg%box_size_lvl(lvl)
1812
1813 ! Geometry subroutines expect this to be set
1814 block => ps(igrid)
1815 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
1816
1817 call get_divb(ps(igrid)%w(ixg^t, 1:nw), ixg^ll, ixm^ll, tmp, &
1819 mg%boxes(id)%cc({1:nc}, mg_irhs) = tmp(ixm^t)
1820 max_divb = max(max_divb, maxval(abs(tmp(ixm^t))))
1821 end do
1822
1823 ! Solve laplacian(phi) = divB
1824 if(stagger_grid) then
1825 call mpi_allreduce(mpi_in_place, max_divb, 1, mpi_double_precision, &
1826 mpi_max, icomm, ierrmpi)
1827
1828 if (mype == 0) print *, "Performing multigrid divB cleaning"
1829 if (mype == 0) print *, "iteration vs residual"
1830 ! Solve laplacian(phi) = divB
1831 do n = 1, max_its
1832 call mg_fas_fmg(mg, n>1, max_res=residual_it(n))
1833 if (mype == 0) write(*, "(I4,E11.3)") n, residual_it(n)
1834 if (residual_it(n) < residual_reduction * max_divb) exit
1835 end do
1836 if (mype == 0 .and. n > max_its) then
1837 print *, "divb_multigrid warning: not fully converged"
1838 print *, "current amplitude of divb: ", residual_it(max_its)
1839 print *, "multigrid smallest grid: ", &
1840 mg%domain_size_lvl(:, mg%lowest_lvl)
1841 print *, "note: smallest grid ideally has <= 8 cells"
1842 print *, "multigrid dx/dy/dz ratio: ", mg%dr(:, 1)/mg%dr(1, 1)
1843 print *, "note: dx/dy/dz should be similar"
1844 end if
1845 else
1846 do n = 1, max_its
1847 call mg_fas_vcycle(mg, max_res=res)
1848 if (res < max_residual) exit
1849 end do
1850 if (res > max_residual) call mpistop("divb_multigrid: no convergence")
1851 end if
1852
1853
1854 ! Correct the magnetic field
1855 do iigrid = 1, igridstail
1856 igrid = igrids(iigrid);
1857 pnode => igrid_to_node(igrid, mype)%node
1858 id = pnode%id
1859
1860 ! Geometry subroutines expect this to be set
1861 block => ps(igrid)
1862 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
1863
1864 ! Compute the gradient of phi
1865 tmp(ix^s) = mg%boxes(id)%cc({:,}, mg_iphi)
1866
1867 if(stagger_grid) then
1868 do idim =1, ndim
1869 ixcmin^d=ixmlo^d-kr(idim,^d);
1870 ixcmax^d=ixmhi^d;
1871 call gradientf(tmp,ps(igrid)%x,ixg^ll,ixc^l,idim,grad(ixg^t,idim))
1872 ps(igrid)%ws(ixc^s,idim)=ps(igrid)%ws(ixc^s,idim)-grad(ixc^s,idim)
1873 end do
1874 call mf_face_to_center(ixm^ll,ps(igrid))
1875 else
1876 do idim = 1, ndim
1877 call gradient(tmp,ixg^ll,ixm^ll,idim,grad(ixg^t, idim))
1878 end do
1879 ! Apply the correction B* = B - gradient(phi)
1880 ps(igrid)%w(ixm^t, mag(1:ndim)) = &
1881 ps(igrid)%w(ixm^t, mag(1:ndim)) - grad(ixm^t, :)
1882 end if
1883
1884 end do
1885
1886 active = .true.
1887
1888 end subroutine mf_clean_divb_multigrid
1889 }
1890
1891 subroutine mf_update_faces(ixI^L,ixO^L,qt,qdt,wprim,fC,fE,sCT,s,vcts)
1893
1894 integer, intent(in) :: ixi^l, ixo^l
1895 double precision, intent(in) :: qt, qdt
1896 ! cell-center primitive variables
1897 double precision, intent(in) :: wprim(ixi^s,1:nw)
1898 type(state) :: sct, s
1899 type(ct_velocity) :: vcts
1900 double precision, intent(in) :: fc(ixi^s,1:nwflux,1:ndim)
1901 double precision, intent(inout) :: fe(ixi^s,sdim:3)
1902
1903 select case(type_ct)
1904 case('average')
1905 call update_faces_average(ixi^l,ixo^l,qt,qdt,fc,fe,sct,s)
1906 case('uct_contact')
1907 call update_faces_contact(ixi^l,ixo^l,qt,qdt,wprim,fc,fe,sct,s,vcts)
1908 case('uct_hll')
1909 call update_faces_hll(ixi^l,ixo^l,qt,qdt,fe,sct,s,vcts)
1910 case default
1911 call mpistop('choose average, uct_contact,or uct_hll for type_ct!')
1912 end select
1913
1914 end subroutine mf_update_faces
1915
1916 !> get electric field though averaging neighors to update faces in CT
1917 subroutine update_faces_average(ixI^L,ixO^L,qt,qdt,fC,fE,sCT,s)
1919 use mod_usr_methods
1920
1921 integer, intent(in) :: ixi^l, ixo^l
1922 double precision, intent(in) :: qt, qdt
1923 type(state) :: sct, s
1924 double precision, intent(in) :: fc(ixi^s,1:nwflux,1:ndim)
1925 double precision, intent(inout) :: fe(ixi^s,sdim:3)
1926
1927 double precision :: circ(ixi^s,1:ndim)
1928 double precision :: e_resi(ixi^s,sdim:3)
1929 integer :: ix^d,ixc^l,ixa^l,i1kr^d,i2kr^d
1930 integer :: idim1,idim2,idir,iwdim1,iwdim2
1931
1932 associate(bfaces=>s%ws,x=>s%x)
1933
1934 ! Calculate contribution to FEM of each edge,
1935 ! that is, estimate value of line integral of
1936 ! electric field in the positive idir direction.
1937
1938 ! if there is resistivity, get eta J
1939 if(mf_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,e_resi)
1940
1941 do idim1=1,ndim
1942 iwdim1 = mag(idim1)
1943 i1kr^d=kr(idim1,^d);
1944 do idim2=1,ndim
1945 iwdim2 = mag(idim2)
1946 i2kr^d=kr(idim2,^d);
1947 do idir=sdim,3! Direction of line integral
1948 ! Allow only even permutations
1949 if (lvc(idim1,idim2,idir)==1) then
1950 ixcmax^d=ixomax^d;
1951 ixcmin^d=ixomin^d+kr(idir,^d)-1;
1952 ! average cell-face electric field to cell edges
1953 {do ix^db=ixcmin^db,ixcmax^db\}
1954 fe(ix^d,idir)=quarter*&
1955 (fc(ix^d,iwdim1,idim2)+fc({ix^d+i1kr^d},iwdim1,idim2)&
1956 -fc(ix^d,iwdim2,idim1)-fc({ix^d+i2kr^d},iwdim2,idim1))
1957 ! add resistive electric field at cell edges E=-vxB+eta J
1958 if(mf_eta/=zero) fe(ix^d,idir)=fe(ix^d,idir)+e_resi(ix^d,idir)
1959 if(record_electric_field) s%we(ix^d,idir)=fe(ix^d,idir)
1960 ! times time step and edge length
1961 fe(ix^d,idir)=fe(ix^d,idir)*qdt*s%dsC(ix^d,idir)
1962 {end do\}
1963 end if
1964 end do
1965 end do
1966 end do
1967
1968 ! allow user to change inductive electric field, especially for boundary driven applications
1969 if(associated(usr_set_electric_field)) &
1970 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
1971
1972 circ(ixi^s,1:ndim)=zero
1973
1974 ! Calculate circulation on each face
1975
1976 do idim1=1,ndim ! Coordinate perpendicular to face
1977 ixcmax^d=ixomax^d;
1978 ixcmin^d=ixomin^d-kr(idim1,^d);
1979 do idim2=1,ndim
1980 ixa^l=ixc^l-kr(idim2,^d);
1981 do idir=sdim,3 ! Direction of line integral
1982 ! Assemble indices
1983 if(lvc(idim1,idim2,idir)==1) then
1984 ! Add line integrals in direction idir
1985 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
1986 +(fe(ixc^s,idir)&
1987 -fe(ixa^s,idir))
1988 else if(lvc(idim1,idim2,idir)==-1) then
1989 ! Add line integrals in direction idir
1990 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
1991 -(fe(ixc^s,idir)&
1992 -fe(ixa^s,idir))
1993 end if
1994 end do
1995 end do
1996 ! Divide by the area of the face to get dB/dt
1997 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
1998 ! Time update
1999 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2000 end do
2001
2002 end associate
2003
2004 end subroutine update_faces_average
2005
2006 !> update faces using UCT contact mode by Gardiner and Stone 2005 JCP 205, 509
2007 subroutine update_faces_contact(ixI^L,ixO^L,qt,qdt,wp,fC,fE,sCT,s,vcts)
2009 use mod_usr_methods
2010
2011 integer, intent(in) :: ixi^l, ixo^l
2012 double precision, intent(in) :: qt, qdt
2013 ! cell-center primitive variables
2014 double precision, intent(in) :: wp(ixi^s,1:nw)
2015 type(state) :: sct, s
2016 type(ct_velocity) :: vcts
2017 double precision, intent(in) :: fc(ixi^s,1:nwflux,1:ndim)
2018 double precision, intent(inout) :: fe(ixi^s,sdim:3)
2019
2020 double precision :: circ(ixi^s,1:ndim)
2021 ! electric field at cell centers
2022 double precision :: ecc(ixi^s,sdim:3)
2023 ! gradient of E at left and right side of a cell face
2024 double precision :: el(ixi^s),er(ixi^s)
2025 ! gradient of E at left and right side of a cell corner
2026 double precision :: elc(ixi^s),erc(ixi^s)
2027 ! current on cell edges
2028 double precision :: jce(ixi^s,sdim:3)
2029 integer :: hxc^l,ixc^l,jxc^l,ixa^l,ixb^l
2030 integer :: idim1,idim2,idir,iwdim1,iwdim2
2031
2032 associate(bfaces=>s%ws,x=>s%x,w=>s%w,vnorm=>vcts%vnorm)
2033
2034 ecc=0.d0
2035 ! Calculate electric field at cell centers
2036 do idim1=1,ndim; do idim2=1,ndim; do idir=sdim,3
2037 if(lvc(idim1,idim2,idir)==1)then
2038 ecc(ixi^s,idir)=ecc(ixi^s,idir)+wp(ixi^s,mag(idim1))*wp(ixi^s,mom(idim2))
2039 else if(lvc(idim1,idim2,idir)==-1) then
2040 ecc(ixi^s,idir)=ecc(ixi^s,idir)-wp(ixi^s,mag(idim1))*wp(ixi^s,mom(idim2))
2041 endif
2042 enddo; enddo; enddo
2043
2044 ! if there is resistivity, get eta J
2045 if(mf_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
2046
2047 ! Calculate contribution to FEM of each edge,
2048 ! that is, estimate value of line integral of
2049 ! electric field in the positive idir direction.
2050 ! evaluate electric field along cell edges according to equation (41)
2051 do idim1=1,ndim
2052 iwdim1 = mag(idim1)
2053 do idim2=1,ndim
2054 iwdim2 = mag(idim2)
2055 do idir=sdim,3 ! Direction of line integral
2056 ! Allow only even permutations
2057 if (lvc(idim1,idim2,idir)==1) then
2058 ixcmax^d=ixomax^d;
2059 ixcmin^d=ixomin^d+kr(idir,^d)-1;
2060 ! Assemble indices
2061 jxc^l=ixc^l+kr(idim1,^d);
2062 hxc^l=ixc^l+kr(idim2,^d);
2063 ! average cell-face electric field to cell edges
2064 fe(ixc^s,idir)=quarter*&
2065 (fc(ixc^s,iwdim1,idim2)+fc(jxc^s,iwdim1,idim2)&
2066 -fc(ixc^s,iwdim2,idim1)-fc(hxc^s,iwdim2,idim1))
2067
2068 ! add slope in idim2 direction from equation (50)
2069 ixamin^d=ixcmin^d;
2070 ixamax^d=ixcmax^d+kr(idim1,^d);
2071 el(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(ixa^s,idir)
2072 hxc^l=ixa^l+kr(idim2,^d);
2073 er(ixa^s)=fc(ixa^s,iwdim1,idim2)-ecc(hxc^s,idir)
2074 where(vnorm(ixc^s,idim1)>0.d0)
2075 elc(ixc^s)=el(ixc^s)
2076 else where(vnorm(ixc^s,idim1)<0.d0)
2077 elc(ixc^s)=el(jxc^s)
2078 else where
2079 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
2080 end where
2081 hxc^l=ixc^l+kr(idim2,^d);
2082 where(vnorm(hxc^s,idim1)>0.d0)
2083 erc(ixc^s)=er(ixc^s)
2084 else where(vnorm(hxc^s,idim1)<0.d0)
2085 erc(ixc^s)=er(jxc^s)
2086 else where
2087 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
2088 end where
2089 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
2090
2091 ! add slope in idim1 direction from equation (50)
2092 jxc^l=ixc^l+kr(idim2,^d);
2093 ixamin^d=ixcmin^d;
2094 ixamax^d=ixcmax^d+kr(idim2,^d);
2095 el(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(ixa^s,idir)
2096 hxc^l=ixa^l+kr(idim1,^d);
2097 er(ixa^s)=-fc(ixa^s,iwdim2,idim1)-ecc(hxc^s,idir)
2098 where(vnorm(ixc^s,idim2)>0.d0)
2099 elc(ixc^s)=el(ixc^s)
2100 else where(vnorm(ixc^s,idim2)<0.d0)
2101 elc(ixc^s)=el(jxc^s)
2102 else where
2103 elc(ixc^s)=0.5d0*(el(ixc^s)+el(jxc^s))
2104 end where
2105 hxc^l=ixc^l+kr(idim1,^d);
2106 where(vnorm(hxc^s,idim2)>0.d0)
2107 erc(ixc^s)=er(ixc^s)
2108 else where(vnorm(hxc^s,idim2)<0.d0)
2109 erc(ixc^s)=er(jxc^s)
2110 else where
2111 erc(ixc^s)=0.5d0*(er(ixc^s)+er(jxc^s))
2112 end where
2113 fe(ixc^s,idir)=fe(ixc^s,idir)+0.25d0*(elc(ixc^s)+erc(ixc^s))
2114
2115 ! add current component of electric field at cell edges E=-vxB+eta J
2116 if(mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2117
2118 if(record_electric_field) s%we(ixc^s,idir)=fe(ixc^s,idir)
2119 ! times time step and edge length
2120 fe(ixc^s,idir)=fe(ixc^s,idir)*qdt*s%dsC(ixc^s,idir)
2121 if (.not.slab) then
2122 where(abs(x(ixc^s,r_)+half*dxlevel(r_))<1.0d-9)
2123 fe(ixc^s,idir)=zero
2124 end where
2125 end if
2126 end if
2127 end do
2128 end do
2129 end do
2130
2131 ! allow user to change inductive electric field, especially for boundary driven applications
2132 if(associated(usr_set_electric_field)) &
2133 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
2134
2135 circ(ixi^s,1:ndim)=zero
2136
2137 ! Calculate circulation on each face
2138 do idim1=1,ndim ! Coordinate perpendicular to face
2139 ixcmax^d=ixomax^d;
2140 ixcmin^d=ixomin^d-kr(idim1,^d);
2141 do idim2=1,ndim
2142 do idir=sdim,3 ! Direction of line integral
2143 ! Assemble indices
2144 if(lvc(idim1,idim2,idir)/=0) then
2145 hxc^l=ixc^l-kr(idim2,^d);
2146 ! Add line integrals in direction idir
2147 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2148 +lvc(idim1,idim2,idir)&
2149 *(fe(ixc^s,idir)&
2150 -fe(hxc^s,idir))
2151 end if
2152 end do
2153 end do
2154 ! Divide by the area of the face to get dB/dt
2155 where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
2156 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2157 elsewhere
2158 circ(ixc^s,idim1)=zero
2159 end where
2160 ! Time update cell-face magnetic field component
2161 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2162 end do
2163
2164 end associate
2165
2166 end subroutine update_faces_contact
2167
2168 !> update faces
2169 subroutine update_faces_hll(ixI^L,ixO^L,qt,qdt,fE,sCT,s,vcts)
2172 use mod_usr_methods
2173
2174 integer, intent(in) :: ixi^l, ixo^l
2175 double precision, intent(in) :: qt, qdt
2176 double precision, intent(inout) :: fe(ixi^s,sdim:3)
2177 type(state) :: sct, s
2178 type(ct_velocity) :: vcts
2179
2180 double precision :: vtill(ixi^s,2)
2181 double precision :: vtilr(ixi^s,2)
2182 double precision :: btill(ixi^s,ndim)
2183 double precision :: btilr(ixi^s,ndim)
2184 double precision :: cp(ixi^s,2)
2185 double precision :: cm(ixi^s,2)
2186 double precision :: circ(ixi^s,1:ndim)
2187 ! current on cell edges
2188 double precision :: jce(ixi^s,sdim:3)
2189 integer :: hxc^l,ixc^l,ixcp^l,jxc^l,ixcm^l
2190 integer :: idim1,idim2,idir
2191
2192 associate(bfaces=>s%ws,bfacesct=>sct%ws,x=>s%x,vbarc=>vcts%vbarC,cbarmin=>vcts%cbarmin,&
2193 cbarmax=>vcts%cbarmax)
2194
2195 ! Calculate contribution to FEM of each edge,
2196 ! that is, estimate value of line integral of
2197 ! electric field in the positive idir direction.
2198
2199 ! Loop over components of electric field
2200
2201 ! idir: electric field component we need to calculate
2202 ! idim1: directions in which we already performed the reconstruction
2203 ! idim2: directions in which we perform the reconstruction
2204
2205 ! if there is resistivity, get eta J
2206 if(mf_eta/=zero) call get_resistive_electric_field(ixi^l,ixo^l,sct,s,jce)
2207
2208 do idir=sdim,3
2209 ! Indices
2210 ! idir: electric field component
2211 ! idim1: one surface
2212 ! idim2: the other surface
2213 ! cyclic permutation: idim1,idim2,idir=1,2,3
2214 ! Velocity components on the surface
2215 ! follow cyclic premutations:
2216 ! Sx(1),Sx(2)=y,z ; Sy(1),Sy(2)=z,x ; Sz(1),Sz(2)=x,y
2217
2218 ixcmax^d=ixomax^d;
2219 ixcmin^d=ixomin^d-1+kr(idir,^d);
2220
2221 ! Set indices and directions
2222 idim1=mod(idir,3)+1
2223 idim2=mod(idir+1,3)+1
2224
2225 jxc^l=ixc^l+kr(idim1,^d);
2226 ixcp^l=ixc^l+kr(idim2,^d);
2227
2228 ! Reconstruct transverse transport velocities
2229 call reconstruct(ixi^l,ixc^l,idim2,vbarc(ixi^s,idim1,1),&
2230 vtill(ixi^s,2),vtilr(ixi^s,2))
2231
2232 call reconstruct(ixi^l,ixc^l,idim1,vbarc(ixi^s,idim2,2),&
2233 vtill(ixi^s,1),vtilr(ixi^s,1))
2234
2235 ! Reconstruct magnetic fields
2236 ! Eventhough the arrays are larger, reconstruct works with
2237 ! the limits ixG.
2238 call reconstruct(ixi^l,ixc^l,idim2,bfacesct(ixi^s,idim1),&
2239 btill(ixi^s,idim1),btilr(ixi^s,idim1))
2240
2241 call reconstruct(ixi^l,ixc^l,idim1,bfacesct(ixi^s,idim2),&
2242 btill(ixi^s,idim2),btilr(ixi^s,idim2))
2243
2244 ! Take the maximum characteristic
2245
2246 cm(ixc^s,1)=max(cbarmin(ixcp^s,idim1),cbarmin(ixc^s,idim1))
2247 cp(ixc^s,1)=max(cbarmax(ixcp^s,idim1),cbarmax(ixc^s,idim1))
2248
2249 cm(ixc^s,2)=max(cbarmin(jxc^s,idim2),cbarmin(ixc^s,idim2))
2250 cp(ixc^s,2)=max(cbarmax(jxc^s,idim2),cbarmax(ixc^s,idim2))
2251
2252
2253 ! Calculate eletric field
2254 fe(ixc^s,idir)=-(cp(ixc^s,1)*vtill(ixc^s,1)*btill(ixc^s,idim2) &
2255 + cm(ixc^s,1)*vtilr(ixc^s,1)*btilr(ixc^s,idim2) &
2256 - cp(ixc^s,1)*cm(ixc^s,1)*(btilr(ixc^s,idim2)-btill(ixc^s,idim2)))&
2257 /(cp(ixc^s,1)+cm(ixc^s,1)) &
2258 +(cp(ixc^s,2)*vtill(ixc^s,2)*btill(ixc^s,idim1) &
2259 + cm(ixc^s,2)*vtilr(ixc^s,2)*btilr(ixc^s,idim1) &
2260 - cp(ixc^s,2)*cm(ixc^s,2)*(btilr(ixc^s,idim1)-btill(ixc^s,idim1)))&
2261 /(cp(ixc^s,2)+cm(ixc^s,2))
2262
2263 ! add current component of electric field at cell edges E=-vxB+eta J
2264 if(mf_eta/=zero) fe(ixc^s,idir)=fe(ixc^s,idir)+jce(ixc^s,idir)
2265
2266 if(record_electric_field) s%we(ixc^s,idir)=fe(ixc^s,idir)
2267 ! times time step and edge length
2268 fe(ixc^s,idir)=qdt*s%dsC(ixc^s,idir)*fe(ixc^s,idir)
2269
2270 if (.not.slab) then
2271 where(abs(x(ixc^s,r_)+half*dxlevel(r_)).lt.1.0d-9)
2272 fe(ixc^s,idir)=zero
2273 end where
2274 end if
2275
2276 end do
2277
2278 ! allow user to change inductive electric field, especially for boundary driven applications
2279 if(associated(usr_set_electric_field)) &
2280 call usr_set_electric_field(ixi^l,ixo^l,qt,qdt,fe,sct)
2281
2282 circ(ixi^s,1:ndim)=zero
2283
2284 ! Calculate circulation on each face: interal(fE dot dl)
2285
2286 do idim1=1,ndim ! Coordinate perpendicular to face
2287 ixcmax^d=ixomax^d;
2288 ixcmin^d=ixomin^d-kr(idim1,^d);
2289 do idim2=1,ndim
2290 do idir=sdim,3 ! Direction of line integral
2291 ! Assemble indices
2292 if(lvc(idim1,idim2,idir)/=0) then
2293 hxc^l=ixc^l-kr(idim2,^d);
2294 ! Add line integrals in direction idir
2295 circ(ixc^s,idim1)=circ(ixc^s,idim1)&
2296 +lvc(idim1,idim2,idir)&
2297 *(fe(ixc^s,idir)&
2298 -fe(hxc^s,idir))
2299 end if
2300 end do
2301 end do
2302 ! Divide by the area of the face to get dB/dt
2303 where(s%surfaceC(ixc^s,idim1) > 1.0d-9*s%dvolume(ixc^s))
2304 circ(ixc^s,idim1)=circ(ixc^s,idim1)/s%surfaceC(ixc^s,idim1)
2305 elsewhere
2306 circ(ixc^s,idim1)=zero
2307 end where
2308 ! Time update
2309 bfaces(ixc^s,idim1)=bfaces(ixc^s,idim1)-circ(ixc^s,idim1)
2310 end do
2311
2312 end associate
2313 end subroutine update_faces_hll
2314
2315 !> calculate eta J at cell edges
2316 subroutine get_resistive_electric_field(ixI^L,ixO^L,sCT,s,jce)
2318 use mod_usr_methods
2319 use mod_geometry
2320
2321 integer, intent(in) :: ixi^l, ixo^l
2322 type(state), intent(in) :: sct, s
2323 ! current on cell edges
2324 double precision :: jce(ixi^s,sdim:3)
2325
2326 ! current on cell centers
2327 double precision :: jcc(ixi^s,7-2*ndir:3)
2328 ! location at cell faces
2329 double precision :: xs(ixgs^t,1:ndim)
2330 ! resistivity
2331 double precision :: eta(ixi^s)
2332 double precision :: gradi(ixgs^t)
2333 integer :: ix^d,ixc^l,ixa^l,ixb^l,idir,idirmin,idim1,idim2
2334
2335 associate(x=>s%x,dx=>s%dx,w=>s%w,wct=>sct%w,wcts=>sct%ws)
2336 ! calculate current density at cell edges
2337 jce=0.d0
2338 do idim1=1,ndim
2339 do idim2=1,ndim
2340 do idir=sdim,3
2341 if (lvc(idim1,idim2,idir)==0) cycle
2342 ixcmax^d=ixomax^d+1;
2343 ixcmin^d=ixomin^d+kr(idir,^d)-2;
2344 ixbmax^d=ixcmax^d-kr(idir,^d)+1;
2345 ixbmin^d=ixcmin^d;
2346 ! current at transverse faces
2347 xs(ixb^s,:)=x(ixb^s,:)
2348 xs(ixb^s,idim2)=x(ixb^s,idim2)+half*dx(ixb^s,idim2)
2349 call gradientf(wcts(ixgs^t,idim2),xs,ixgs^ll,ixc^l,idim1,gradi)
2350 if (lvc(idim1,idim2,idir)==1) then
2351 jce(ixc^s,idir)=jce(ixc^s,idir)+gradi(ixc^s)
2352 else
2353 jce(ixc^s,idir)=jce(ixc^s,idir)-gradi(ixc^s)
2354 end if
2355 end do
2356 end do
2357 end do
2358 ! get resistivity
2359 if(mf_eta>zero)then
2360 jce(ixi^s,:)=jce(ixi^s,:)*mf_eta
2361 else
2362 ixa^l=ixo^l^ladd1;
2363 call get_current(wct,ixi^l,ixa^l,idirmin,jcc)
2364 call usr_special_resistivity(wct,ixi^l,ixa^l,idirmin,x,jcc,eta)
2365 ! calcuate eta on cell edges
2366 do idir=sdim,3
2367 ixcmax^d=ixomax^d;
2368 ixcmin^d=ixomin^d+kr(idir,^d)-1;
2369 jcc(ixc^s,idir)=0.d0
2370 {do ix^db=0,1\}
2371 if({ ix^d==1 .and. ^d==idir | .or.}) cycle
2372 ixamin^d=ixcmin^d+ix^d;
2373 ixamax^d=ixcmax^d+ix^d;
2374 jcc(ixc^s,idir)=jcc(ixc^s,idir)+eta(ixa^s)
2375 {end do\}
2376 jcc(ixc^s,idir)=jcc(ixc^s,idir)*0.25d0
2377 jce(ixc^s,idir)=jce(ixc^s,idir)*jcc(ixc^s,idir)
2378 enddo
2379 end if
2380
2381 end associate
2382 end subroutine get_resistive_electric_field
2383
2384 !> calculate cell-center values from face-center values
2385 subroutine mf_face_to_center(ixO^L,s)
2387 ! Non-staggered interpolation range
2388 integer, intent(in) :: ixo^l
2389 type(state) :: s
2390
2391 integer :: ix^d
2392
2393 ! calculate cell-center values from face-center values in 2nd order
2394 {!dir$ ivdep
2395 do ix^db=ixomin^db,ixomax^db\}
2396 {^ifthreed
2397 s%w(ix^d,mag(1))=half/s%surface(ix^d,1)*(s%ws(ix^d,1)*s%surfaceC(ix^d,1)&
2398 +s%ws(ix1-1,ix2,ix3,1)*s%surfaceC(ix1-1,ix2,ix3,1))
2399 s%w(ix^d,mag(2))=half/s%surface(ix^d,2)*(s%ws(ix^d,2)*s%surfaceC(ix^d,2)&
2400 +s%ws(ix1,ix2-1,ix3,2)*s%surfaceC(ix1,ix2-1,ix3,2))
2401 s%w(ix^d,mag(3))=half/s%surface(ix^d,3)*(s%ws(ix^d,3)*s%surfaceC(ix^d,3)&
2402 +s%ws(ix1,ix2,ix3-1,3)*s%surfaceC(ix1,ix2,ix3-1,3))
2403 }
2404 {^iftwod
2405 s%w(ix^d,mag(1))=half/s%surface(ix^d,1)*(s%ws(ix^d,1)*s%surfaceC(ix^d,1)&
2406 +s%ws(ix1-1,ix2,1)*s%surfaceC(ix1-1,ix2,1))
2407 s%w(ix^d,mag(2))=half/s%surface(ix^d,2)*(s%ws(ix^d,2)*s%surfaceC(ix^d,2)&
2408 +s%ws(ix1,ix2-1,2)*s%surfaceC(ix1,ix2-1,2))
2409 }
2410 {end do\}
2411
2412 ! calculate cell-center values from face-center values in 4th order
2413 !do idim=1,ndim
2414 ! gxO^L=ixO^L-2*kr(idim,^D);
2415 ! hxO^L=ixO^L-kr(idim,^D);
2416 ! jxO^L=ixO^L+kr(idim,^D);
2417
2418 ! ! Interpolate to cell barycentre using fourth order central formula
2419 ! w(ixO^S,mag(idim))=(0.0625d0/s%surface(ixO^S,idim))*&
2420 ! ( -ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
2421 ! +9.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
2422 ! +9.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
2423 ! -ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) )
2424 !end do
2425
2426 ! calculate cell-center values from face-center values in 6th order
2427 !do idim=1,ndim
2428 ! fxO^L=ixO^L-3*kr(idim,^D);
2429 ! gxO^L=ixO^L-2*kr(idim,^D);
2430 ! hxO^L=ixO^L-kr(idim,^D);
2431 ! jxO^L=ixO^L+kr(idim,^D);
2432 ! kxO^L=ixO^L+2*kr(idim,^D);
2433
2434 ! ! Interpolate to cell barycentre using sixth order central formula
2435 ! w(ixO^S,mag(idim))=(0.00390625d0/s%surface(ixO^S,idim))* &
2436 ! ( +3.0d0*ws(fxO^S,idim)*s%surfaceC(fxO^S,idim) &
2437 ! -25.0d0*ws(gxO^S,idim)*s%surfaceC(gxO^S,idim) &
2438 ! +150.0d0*ws(hxO^S,idim)*s%surfaceC(hxO^S,idim) &
2439 ! +150.0d0*ws(ixO^S,idim)*s%surfaceC(ixO^S,idim) &
2440 ! -25.0d0*ws(jxO^S,idim)*s%surfaceC(jxO^S,idim) &
2441 ! +3.0d0*ws(kxO^S,idim)*s%surfaceC(kxO^S,idim) )
2442 !end do
2443
2444 end subroutine mf_face_to_center
2445
2446 !> calculate magnetic field from vector potential
2447 subroutine b_from_vector_potential(ixIs^L, ixI^L, ixO^L, ws, x)
2450
2451 integer, intent(in) :: ixis^l, ixi^l, ixo^l
2452 double precision, intent(inout) :: ws(ixis^s,1:nws)
2453 double precision, intent(in) :: x(ixi^s,1:ndim)
2454
2455 double precision :: adummy(ixis^s,1:3)
2456
2457 call b_from_vector_potentiala(ixis^l, ixi^l, ixo^l, ws, x, adummy)
2458
2459 end subroutine b_from_vector_potential
2460
2463
2464 double precision :: sum_jbb_ipe, sum_j_ipe, sum_l_ipe, f_i_ipe, volume_pe
2465 double precision :: sum_jbb, sum_j, sum_l, f_i, volume, cw_sin_theta
2466 integer :: iigrid, igrid, ix^d
2467 integer :: amode, istatus(mpi_status_size)
2468 integer, save :: fhmf
2469 logical :: patchwi(ixg^t)
2470 logical, save :: logmfopened=.false.
2471 character(len=800) :: filename,filehead
2472 character(len=800) :: line,datastr
2473
2474 sum_jbb_ipe = 0.d0
2475 sum_j_ipe = 0.d0
2476 sum_l_ipe = 0.d0
2477 f_i_ipe = 0.d0
2478 volume_pe=0.d0
2479 do iigrid=1,igridstail; igrid=igrids(iigrid);
2480 block=>ps(igrid)
2481 ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
2482 call mask_inner(ixg^ll,ixm^ll,ps(igrid)%w,ps(igrid)%x,patchwi,volume_pe)
2483 sum_jbb_ipe = sum_jbb_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
2484 ps(igrid)%x,1,patchwi)
2485 sum_j_ipe = sum_j_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
2486 ps(igrid)%x,2,patchwi)
2487 f_i_ipe=f_i_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
2488 ps(igrid)%x,3,patchwi)
2489 sum_l_ipe = sum_l_ipe+integral_grid_mf(ixg^ll,ixm^ll,ps(igrid)%w,&
2490 ps(igrid)%x,4,patchwi)
2491 end do
2492 call mpi_reduce(sum_jbb_ipe,sum_jbb,1,mpi_double_precision,&
2493 mpi_sum,0,icomm,ierrmpi)
2494 call mpi_reduce(sum_j_ipe,sum_j,1,mpi_double_precision,mpi_sum,0,&
2495 icomm,ierrmpi)
2496 call mpi_reduce(f_i_ipe,f_i,1,mpi_double_precision,mpi_sum,0,&
2497 icomm,ierrmpi)
2498 call mpi_reduce(sum_l_ipe,sum_l,1,mpi_double_precision,mpi_sum,0,&
2499 icomm,ierrmpi)
2500 call mpi_reduce(volume_pe,volume,1,mpi_double_precision,mpi_sum,0,&
2501 icomm,ierrmpi)
2502
2503 if(mype==0) then
2504 ! current- and volume-weighted average of the sine of the angle
2505 ! between the magnetic field B and the current density J
2506 cw_sin_theta = sum_jbb/sum_j
2507 ! volume-weighted average of the absolute value of the fractional
2508 ! magnetic flux change
2509 f_i = f_i/volume
2510 sum_j=sum_j/volume
2511 sum_l=sum_l/volume
2512 if(.not.logmfopened) then
2513 ! generate filename
2514 write(filename,"(a,a)") trim(base_filename), "_mf_metrics.csv"
2515
2516 ! Delete the log when not doing a restart run
2517 if (restart_from_file == undefined) then
2518 open(unit=20,file=trim(filename),status='replace')
2519 close(20, status='delete')
2520 end if
2521
2522 amode=ior(mpi_mode_create,mpi_mode_wronly)
2523 amode=ior(amode,mpi_mode_append)
2524 call mpi_file_open(mpi_comm_self,filename,amode,mpi_info_null,fhmf,ierrmpi)
2525 logmfopened=.true.
2526 filehead=" it,time,CW_sin_theta,current,Lorenz_force,f_i"
2527 if (it == 0) then
2528 call mpi_file_write(fhmf,filehead,len_trim(filehead), &
2529 mpi_character,istatus,ierrmpi)
2530 call mpi_file_write(fhmf,achar(10),1,mpi_character,istatus,ierrmpi)
2531 end if
2532 end if
2533 line=''
2534 write(datastr,'(i6,a)') it,','
2535 line=trim(line)//trim(datastr)
2536 write(datastr,'(es13.6,a)') global_time,','
2537 line=trim(line)//trim(datastr)
2538 write(datastr,'(es13.6,a)') cw_sin_theta,','
2539 line=trim(line)//trim(datastr)
2540 write(datastr,'(es13.6,a)') sum_j,','
2541 line=trim(line)//trim(datastr)
2542 write(datastr,'(es13.6,a)') sum_l,','
2543 line=trim(line)//trim(datastr)
2544 write(datastr,'(es13.6)') f_i
2545 line=trim(line)//trim(datastr)//new_line('A')
2546 call mpi_file_write(fhmf,line,len_trim(line),mpi_character,istatus,ierrmpi)
2547 if(.not.time_advance) then
2548 call mpi_file_close(fhmf,ierrmpi)
2549 logmfopened=.false.
2550 end if
2551 end if
2552
2553 end subroutine record_force_free_metrics
2554
2555 subroutine mask_inner(ixI^L,ixO^L,w,x,patchwi,volume_pe)
2557
2558 integer, intent(in) :: ixi^l,ixo^l
2559 double precision, intent(in):: w(ixi^s,nw),x(ixi^s,1:ndim)
2560 double precision, intent(inout) :: volume_pe
2561 logical, intent(out) :: patchwi(ixi^s)
2562
2563 double precision :: xo^l
2564 integer :: ix^d
2565
2566 {xomin^d = xprobmin^d + 0.05d0*(xprobmax^d-xprobmin^d)\}
2567 {xomax^d = xprobmax^d - 0.05d0*(xprobmax^d-xprobmin^d)\}
2568 if(slab) then
2569 xomin^nd = xprobmin^nd
2570 else
2571 xomin1 = xprobmin1
2572 end if
2573
2574 {do ix^db=ixomin^db,ixomax^db\}
2575 if({ x(ix^dd,^d) > xomin^d .and. x(ix^dd,^d) < xomax^d | .and. }) then
2576 patchwi(ix^d)=.true.
2577 volume_pe=volume_pe+block%dvolume(ix^d)
2578 else
2579 patchwi(ix^d)=.false.
2580 endif
2581 {end do\}
2582
2583 end subroutine mask_inner
2584
2585 function integral_grid_mf(ixI^L,ixO^L,w,x,iw,patchwi)
2587
2588 integer, intent(in) :: ixi^l,ixo^l,iw
2589 double precision, intent(in) :: x(ixi^s,1:ndim)
2590 double precision :: w(ixi^s,nw+nwauxio)
2591 logical, intent(in) :: patchwi(ixi^s)
2592
2593 double precision, dimension(ixI^S,7-2*ndir:3) :: current
2594 double precision, dimension(ixI^S,1:ndir) :: bvec,qvec
2595 double precision :: integral_grid_mf,tmp(ixi^s),bm2
2596 integer :: ix^d,idirmin0,idirmin,idir,jdir,kdir
2597
2598 integral_grid_mf=0.d0
2599 select case(iw)
2600 case(1)
2601 ! Sum(|JxB|/|B|*dvolume)
2602 bvec(ixo^s,:)=w(ixo^s,mag(:))
2603 call get_current(w,ixi^l,ixo^l,idirmin,current)
2604 ! calculate Lorentz force
2605 qvec(ixo^s,1:ndir)=zero
2606 do idir=1,ndir; do jdir=idirmin,3; do kdir=1,ndir
2607 if(lvc(idir,jdir,kdir)/=0)then
2608 if(lvc(idir,jdir,kdir)==1)then
2609 qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
2610 else
2611 qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
2612 end if
2613 end if
2614 end do; end do; end do
2615
2616 {do ix^db=ixomin^db,ixomax^db\}
2617 if(patchwi(ix^d)) then
2618 bm2=sum(bvec(ix^d,:)**2)
2619 if(bm2/=0.d0) bm2=1.d0/bm2
2620 integral_grid_mf=integral_grid_mf+sqrt(sum(qvec(ix^d,:)**2)*&
2621 bm2)*block%dvolume(ix^d)
2622 end if
2623 {end do\}
2624 case(2)
2625 ! Sum(|J|*dvolume)
2626 call get_current(w,ixi^l,ixo^l,idirmin,current)
2627 {do ix^db=ixomin^db,ixomax^db\}
2628 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+sqrt(sum(current(ix^d,:)**2))*&
2629 block%dvolume(ix^d)
2630 {end do\}
2631 case(3)
2632 ! f_i solenoidal property of B: (dvolume |div B|)/(dsurface |B|)
2633 ! Sum(f_i*dvolume)
2634 call get_normalized_divb(w,ixi^l,ixo^l,tmp)
2635 {do ix^db=ixomin^db,ixomax^db\}
2636 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+tmp(ix^d)*block%dvolume(ix^d)
2637 {end do\}
2638 case(4)
2639 ! Sum(|JxB|*dvolume)
2640 bvec(ixo^s,:)=w(ixo^s,mag(:))
2641 call get_current(w,ixi^l,ixo^l,idirmin,current)
2642 ! calculate Lorentz force
2643 qvec(ixo^s,1:ndir)=zero
2644 do idir=1,ndir; do jdir=idirmin,3; do kdir=1,ndir
2645 if(lvc(idir,jdir,kdir)/=0)then
2646 if(lvc(idir,jdir,kdir)==1)then
2647 qvec(ixo^s,idir)=qvec(ixo^s,idir)+current(ixo^s,jdir)*bvec(ixo^s,kdir)
2648 else
2649 qvec(ixo^s,idir)=qvec(ixo^s,idir)-current(ixo^s,jdir)*bvec(ixo^s,kdir)
2650 end if
2651 end if
2652 end do; end do; end do
2653
2654 {do ix^db=ixomin^db,ixomax^db\}
2655 if(patchwi(ix^d)) integral_grid_mf=integral_grid_mf+&
2656 sqrt(sum(qvec(ix^d,:)**2))*block%dvolume(ix^d)
2657 {end do\}
2658 end select
2659 return
2660 end function integral_grid_mf
2661
2662end module mod_mf_phys
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
subroutine reconstruct(ixil, ixcl, idir, q, ql, qr)
Reconstruct scalar q within ixO^L to 1/2 dx in direction idir Return both left and right reconstructe...
subroutine b_from_vector_potentiala(ixisl, ixil, ixol, ws, x, a)
calculate magnetic field from vector potential A at cell edges
Module with basic grid data structures.
Definition mod_forest.t:2
type(tree_node_ptr), dimension(:,:), allocatable, save igrid_to_node
Array to go from an [igrid, ipe] index to a node pointer.
Definition mod_forest.t:32
subroutine, public get_divb(w, ixil, ixol, divb, nth_in)
Calculate div B within ixO.
integer, dimension(:), allocatable, public mag
Indices of the magnetic field.
Module with geometry-related routines (e.g., divergence, curl)
Definition mod_geometry.t:2
integer coordinate
Definition mod_geometry.t:7
integer, parameter cylindrical
subroutine curlvector(qvec, ixil, ixol, curlvec, idirmin, idirmin0, ndir0, fourthorder)
Calculate curl of a vector qvec within ixL Options to employ standard second order CD evaluations use...
subroutine gradient(q, ixil, ixol, idir, gradq, nth_in)
subroutine gradientf(q, x, ixil, ixol, idir, gradq, nth_in, pm_in)
subroutine gradientl(q, ixil, ixol, idir, gradq)
update ghost cells of all blocks including physical boundaries
integer, dimension(0:3^d &), target type_recv_p_p1
integer, dimension( :^d &), pointer type_recv_r
integer, dimension(-1:1^d &), target type_send_srl_p1
subroutine getbc(time, qdt, psb, nwstart, nwbc)
do update ghost cells of all blocks including physical boundaries
integer, dimension(0:3^d &), target type_send_p_f
integer, dimension( :^d &), pointer type_send_p
integer, dimension( :^d &), pointer type_send_srl
integer, dimension(-1:1^d &), target type_send_r_p1
subroutine create_bc_mpi_datatype(nwstart, nwbc)
integer, dimension(0:3^d &), target type_recv_r_f
integer, dimension(-1:1^d &), target type_recv_srl_f
integer, dimension(-1:1^d &), target type_send_r_f
integer, dimension(-1:1^d &), target type_send_srl_f
integer, dimension(0:3^d &), target type_send_p_p1
integer, dimension( :^d &), pointer type_send_r
integer, dimension(0:3^d &), target type_recv_r_p1
integer, dimension( :^d &), pointer type_recv_p
integer, dimension(-1:1^d &), target type_recv_srl_p1
integer, dimension(0:3^d &), target type_recv_p_f
integer, dimension( :^d &), pointer type_recv_srl
This module contains definitions of global parameters and variables and some generic functions/subrou...
type(state), pointer block
Block pointer for using one block and its previous state.
double precision dtdiffpar
For resistive MHD, the time step is also limited by the diffusion time: .
character(len=std_len) typegrad
double precision unit_charge
Physical scaling factor for charge.
integer ixghi
Upper index of grid block arrays.
integer, dimension(3, 3, 3) lvc
Levi-Civita tensor.
double precision unit_time
Physical scaling factor for time.
double precision unit_density
Physical scaling factor for density.
integer, parameter unitpar
file handle for IO
double precision global_time
The global simulation time.
double precision unit_mass
Physical scaling factor for mass.
integer, dimension(3, 3) kr
Kronecker delta tensor.
integer it
Number of time steps taken.
integer, dimension(:, :), allocatable typeboundary
Array indicating the type of boundary condition per variable and per physical boundary.
double precision unit_numberdensity
Physical scaling factor for number density.
double precision unit_pressure
Physical scaling factor for pressure.
integer, parameter ndim
Number of spatial dimensions for grid variables.
double precision unit_length
Physical scaling factor for length.
logical stagger_grid
True for using stagger grid.
double precision cmax_global
global fastest wave speed needed in fd scheme and glm method
logical use_particles
Use particles module or not.
character(len=std_len), dimension(:), allocatable par_files
Which par files are used as input.
integer icomm
The MPI communicator.
double precision bdip
amplitude of background dipolar, quadrupolar, octupolar, user's field
integer mype
The rank of the current MPI task.
double precision, dimension(:), allocatable, parameter d
integer ndir
Number of spatial dimensions (components) for vector variables.
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
logical slab
Cartesian geometry or not.
integer, parameter bc_periodic
integer, parameter bc_special
boundary condition types
double precision unit_magneticfield
Physical scaling factor for magnetic field.
integer nwauxio
Number of auxiliary variables that are only included in output.
double precision unit_velocity
Physical scaling factor for velocity.
logical time_advance
do time evolving
character(len=std_len) restart_from_file
If not 'unavailable', resume from snapshot with this base file name.
double precision c_norm
Normalised speed of light.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
double precision unit_temperature
Physical scaling factor for temperature.
logical si_unit
Use SI units (.true.) or use cgs units (.false.)
double precision, dimension(:,:), allocatable dx
integer nghostcells
Number of ghost cells surrounding a grid.
integer, parameter sdim
starting dimension for electric field
character(len= *), parameter undefined
logical need_global_cmax
need global maximal wave speed
double precision, dimension(^nd) dxlevel
store unstretched cell size of current level
logical use_multigrid
Use multigrid (only available in 2D and 3D)
logical slab_uniform
uniform Cartesian geometry or not (stretched Cartesian)
character(len=std_len) base_filename
Base file name for simulation output, which will be followed by a number.
integer max_blocks
The maximum number of grid blocks in a processor.
integer r_
Indices for cylindrical coordinates FOR TESTS, negative value when not used:
logical record_electric_field
True for record electric field.
integer, parameter ixglo
Lower index of grid block arrays (always 1)
Magnetofriction module.
Definition mod_mf_phys.t:2
integer, public, protected b
Definition mod_mf_phys.t:46
double precision, public mf_nu
viscosity coefficient in s cm^-2 for solar corona (Cheung 2012 ApJ)
Definition mod_mf_phys.t:12
subroutine, public mf_phys_init()
integer, public, protected mf_divb_nth
Whether divB is computed with a fourth order approximation.
Definition mod_mf_phys.t:82
subroutine, public b_from_vector_potential(ixisl, ixil, ixol, ws, x)
calculate magnetic field from vector potential
logical, public, protected mf_record_electric_field
set to true if need to record electric field on cell edges
Definition mod_mf_phys.t:79
logical, public, protected clean_initial_divb
clean divb in the initial condition
Definition mod_mf_phys.t:94
logical, public divbwave
Add divB wave in Roe solver.
Definition mod_mf_phys.t:91
subroutine, public mf_face_to_center(ixol, s)
calculate cell-center values from face-center values
double precision, public mf_vmax
maximal limit of magnetofrictional velocity in cm s^-1 (Pomoell 2019 A&A)
Definition mod_mf_phys.t:15
logical, public, protected mf_glm
Whether GLM-MHD is used.
Definition mod_mf_phys.t:70
double precision, public mf_eta_hyper
The hyper-resistivity.
Definition mod_mf_phys.t:28
double precision, public mf_glm_alpha
GLM-MHD parameter: ratio of the diffusive and advective time scales for div b taking values within [0...
Definition mod_mf_phys.t:22
integer, public, protected c
Indices of the momentum density for the form of better vectorization.
Definition mod_mf_phys.t:43
integer, public, protected c_
Definition mod_mf_phys.t:43
logical, dimension(2 *^nd), public, protected boundary_divbfix
To control divB=0 fix for boundary.
Definition mod_mf_phys.t:85
character(len=std_len), public, protected type_ct
Method type of constrained transport.
character(len=std_len), public, protected typedivbfix
Method type to clean divergence of B.
Definition mod_mf_phys.t:97
subroutine, public get_normalized_divb(w, ixil, ixol, divb)
get dimensionless div B = |divB| * volume / area / |B|
integer, dimension(2 *^nd), public, protected boundary_divbfix_skip
To skip * layer of ghost cells during divB=0 fix for boundary.
Definition mod_mf_phys.t:52
double precision, dimension(2 *^nd), public mf_decay_scale
decay scale of frictional velocity near boundaries
Definition mod_mf_phys.t:18
subroutine, public record_force_free_metrics()
subroutine, public mf_to_conserved(ixil, ixol, w, x)
Transform primitive variables into conservative ones.
logical, public, protected mf_4th_order
MHD fourth order.
Definition mod_mf_phys.t:76
subroutine, public get_current(w, ixil, ixol, idirmin, current, fourthorder)
Calculate idirmin and the idirmin:3 components of the common current array make sure that dxlevel(^D)...
integer, public, protected psi_
Indices of the GLM psi.
Definition mod_mf_phys.t:49
integer, dimension(:), allocatable, public, protected mom
Indices of the momentum density.
Definition mod_mf_phys.t:40
subroutine, public mf_clean_divb_multigrid(qdt, qt, active)
logical, public, protected mf_particles
Whether particles module is added.
Definition mod_mf_phys.t:67
integer, public, protected m
Definition mod_mf_phys.t:43
subroutine, public mf_to_primitive(ixil, ixol, w, x)
Transform conservative variables into primitive ones.
double precision, public mf_eta
The resistivity.
Definition mod_mf_phys.t:25
logical, public, protected source_split_divb
Whether divB cleaning sources are added splitting from fluid solver.
Definition mod_mf_phys.t:73
double precision, public, protected he_abundance
Helium abundance over Hydrogen.
Definition mod_mf_phys.t:34
Module to couple the octree-mg library to AMRVAC. This file uses the VACPP preprocessor,...
type(mg_t) mg
Data structure containing the multigrid tree.
Module containing all the particle routines.
subroutine particles_init()
Initialize particle data and parameters.
This module defines the procedures of a physics module. It contains function pointers for the various...
Definition mod_physics.t:4
procedure(sub_get_ct_velocity), pointer phys_get_ct_velocity
Definition mod_physics.t:79
procedure(sub_convert), pointer phys_to_primitive
Definition mod_physics.t:55
integer, parameter flux_tvdlf
Indicates the flux should be treated with tvdlf.
Definition mod_physics.t:26
procedure(sub_write_info), pointer phys_write_info
Definition mod_physics.t:77
procedure(sub_get_flux), pointer phys_get_flux
Definition mod_physics.t:64
procedure(sub_get_cbounds), pointer phys_get_cbounds
Definition mod_physics.t:63
procedure(sub_get_dt), pointer phys_get_dt
Definition mod_physics.t:67
procedure(sub_add_source_geom), pointer phys_add_source_geom
Definition mod_physics.t:68
procedure(sub_check_params), pointer phys_check_params
Definition mod_physics.t:52
integer, parameter flux_default
Indicates a normal flux.
Definition mod_physics.t:24
integer, dimension(:, :), allocatable flux_type
Array per direction per variable, which can be used to specify that certain fluxes have to be treated...
Definition mod_physics.t:21
procedure(sub_update_faces), pointer phys_update_faces
Definition mod_physics.t:80
procedure(sub_convert), pointer phys_to_conserved
Definition mod_physics.t:54
character(len=name_len) physics_type
String describing the physics type of the simulation.
Definition mod_physics.t:50
procedure(sub_clean_divb), pointer phys_clean_divb
Definition mod_physics.t:84
procedure(sub_boundary_adjust), pointer phys_boundary_adjust
Definition mod_physics.t:76
procedure(sub_global_source), pointer phys_global_source_after
Definition mod_physics.t:70
procedure(sub_add_source), pointer phys_add_source
Definition mod_physics.t:69
procedure(sub_face_to_center), pointer phys_face_to_center
Definition mod_physics.t:81
procedure(sub_modify_wlr), pointer phys_modify_wlr
Definition mod_physics.t:56
procedure(sub_get_cmax), pointer phys_get_cmax
Definition mod_physics.t:57
logical phys_energy
Solve energy equation or not.
Definition mod_physics.t:35
procedure(sub_special_advance), pointer phys_special_advance
Definition mod_physics.t:71
Module with all the methods that users can customize in AMRVAC.
procedure(special_resistivity), pointer usr_special_resistivity
procedure(set_electric_field), pointer usr_set_electric_field
integer nw
Total number of variables.
integer number_species
number of species: each species has different characterictic speeds and should be used accordingly in...
The data structure that contains information about a tree node/grid block.
Definition mod_forest.t:11