MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_lfff.t
Go to the documentation of this file.
1 !> Program to extrapolate linear force-free fields in 3D Cartesian coordinates,
2 !> based on exact Green function method (Chiu & Hilton 1977 ApJ 212,873).
3 !>
4 !> Usage:
5 !> 1 In the subroutine usr_set_parameters of mod_usr.t:
6 !> To extrapolate a linear force free field from a observed magnetogram
7 !> prepared in a data file, e.g., 'hmiM720sxxxx.dat' replace
8 !> call init_bc_fff_data('hmiM720sxxxx.dat',unit_length,unit_magneticfield)
9 !> 'hmiM720sxxxx.dat' must be a binary file containing nx1,nx2,xc1,xc2,dxm1,
10 !> dxm2, Bz0(nx1,nx2). Integers nx1 and nx2 give the resolution of the
11 !> uniform-grid magentogram. Others are double-precision floats. xc1 and xc2
12 !> are coordinates of the central point of the magnetogram. dxm1 and dxm2
13 !> are the cell sizes for each direction, Bz0 is the vertical conponent
14 !> of magetic field on the solar surface from observations.
15 !>2 In the subroutine usr_init_one_grid of mod_usr.t,
16 !> add lines like:
17 !>
18 !> double precision :: Bf(ixG^S,1:ndir), alpha, zshift
19 !>
20 !> alpha=0.d0 ! potential field
21 !> !alpha=0.08d0 ! non-potential linear force-free field
22 !> zshift=0.05d0 ! lift your box zshift heigher to the bottom magnetogram
23 !> call calc_lin_fff(ixG^L,ix^L,Bf,x,alpha,zshift)
24 !>
25 !>3 Notice that the resolution of input magnetogram must be better than the best
26 !> resolution of your AMR grid to have a good behavior close to the bottom layer
27 module mod_lfff
28  implicit none
29 
30  integer, save :: nx1,nx2
31  double precision, save :: bzmax,darea
32  double precision, allocatable, save :: bz0(:,:)
33  double precision, allocatable, save :: xa1(:),xa2(:)
34 
35 contains
36 
37  subroutine init_b_fff_data(magnetogramname,qLunit,qBunit)
39  use mod_comm_lib, only: mpistop
40 
41  double precision, intent(in) :: qLunit,qBunit
42  double precision :: xc1,xc2,dxm1,dxm2
43  integer, dimension(MPI_STATUS_SIZE) :: statuss
44  integer :: file_handle,i
45  character(len=*), intent(in) :: magnetogramname
46  logical :: aexist
47  ! nx1,nx2 are numbers of cells for each direction
48  ! xc1,xc2 are coordinates of the central point of the magnetogram
49  ! dxm1,dxm2 are cell sizes for each direction
50  ! Bz0 is the 2D Bz magnetogram
51  inquire(file=magnetogramname,exist=aexist)
52  if(.not. aexist) then
53  if(mype==0) write(*,'(2a)') "can not find file:",magnetogramname
54  call mpistop("no input magnetogram----init_b_fff_data")
55  end if
56  call mpi_file_open(icomm,magnetogramname,mpi_mode_rdonly,mpi_info_null,&
57  file_handle,ierrmpi)
58  call mpi_file_read_all(file_handle,nx1,1,mpi_integer,statuss,ierrmpi)
59  call mpi_file_read_all(file_handle,nx2,1,mpi_integer,statuss,ierrmpi)
60  allocate(bz0(nx1,nx2))
61  call mpi_file_read_all(file_handle,xc1,1,mpi_double_precision,statuss,ierrmpi)
62  call mpi_file_read_all(file_handle,xc2,1,mpi_double_precision,statuss,ierrmpi)
63  call mpi_file_read_all(file_handle,dxm1,1,mpi_double_precision,statuss,ierrmpi)
64  call mpi_file_read_all(file_handle,dxm2,1,mpi_double_precision,statuss,ierrmpi)
65  call mpi_file_read_all(file_handle,bz0,nx1*nx2,mpi_double_precision,&
66  statuss,ierrmpi)
67  call mpi_file_close(file_handle,ierrmpi)
68  allocate(xa1(nx1))
69  allocate(xa2(nx2))
70  do i=1,nx1
71  xa1(i) = xc1 + (dble(i) - dble(nx1)/2.d0 - 0.5d0)*dxm1
72  enddo
73  do i=1,nx2
74  xa2(i) = xc2 + (dble(i) - dble(nx2)/2.d0 - 0.5d0)*dxm2
75  enddo
76  ! declare and define global variables Lunit and Bunit to be your length unit in
77  ! cm and magnetic strength unit in Gauss first
78  dxm1=dxm1/qlunit
79  dxm2=dxm2/qlunit
80  xa1=xa1/qlunit
81  xa2=xa2/qlunit
82  darea=dxm1*dxm2
83  bz0=bz0/qbunit
84  bzmax=maxval(dabs(bz0(:,:)))
85 
86  ! normalize b
87  bz0=bz0/bzmax
88  if(mype==0) then
89  print*,'magnetogram xrange:',minval(xa1),maxval(xa1)
90  print*,'magnetogram yrange:',minval(xa2),maxval(xa2)
91  end if
92 
93  if(mype==0) then
94  print*,'extrapolating 3D force-free field from an observed Bz '
95  print*,'magnetogram of',nx1,'by',nx2,'pixels. Bzmax=',bzmax
96  endif
97 
98  end subroutine init_b_fff_data
99 
100 {^ifthreed
101  subroutine calc_lin_fff(ixI^L,ixO^L,Bf,x,alpha,zshift,idir)
102  ! PURPOSE:
103  ! Calculation to determine linear FFF from the field on
104  ! the lower boundary (Chiu and Hilton 1977 ApJ 212,873).
105  ! NOTE: Only works for Cartesian coordinates
106  ! INPUT: Bf,x
107  ! OUTPUT: updated b in w
109 
110  integer, intent(in) :: ixI^L, ixO^L
111  integer, optional, intent(in) :: idir
112  double precision, intent(in) :: x(ixI^S,1:ndim),alpha,zshift
113  double precision, intent(inout) :: Bf(ixI^S,1:ndir)
114 
115  double precision, dimension(ixO^S) :: cos_az,sin_az,zk,bigr,r,r2,r3,cos_ar,sin_ar,g,dgdz
116  double precision :: gx(ixO^S,1:ndim),twopiinv
117  integer :: idim,ixp1,ixp2
118 
119  bf=0.d0
120  twopiinv = 0.5d0/dpi*bzmax*darea
121  ! get cos and sin arrays
122  zk(ixo^s)=x(ixo^s,3)-xprobmin3+zshift
123  cos_az(ixo^s)=dcos(alpha*zk(ixo^s))
124  sin_az(ixo^s)=dsin(alpha*zk(ixo^s))
125  ! looping Bz0 pixels
126  do ixp2=1,nx2
127  do ixp1=1,nx1
128  bigr(ixo^s)=dsqrt((x(ixo^s,1)-xa1(ixp1))**2+&
129  (x(ixo^s,2)-xa2(ixp2))**2)
130  r2=bigr**2+zk**2
131  r=dsqrt(r2)
132  r3=r**3
133  cos_ar=dcos(alpha*r)
134  sin_ar=dsin(alpha*r)
135  where(bigr/=0.d0)
136  bigr=1.d0/bigr
137  end where
138  where(r/=0.d0)
139  r=1.d0/r
140  end where
141  where(r2/=0.d0)
142  r2=1.d0/r2
143  end where
144  where(r3/=0.d0)
145  r3=1.d0/r3
146  end where
147  g=(zk*cos_ar*r-cos_az)*bigr
148  dgdz=(cos_ar*(r-zk**2*r3)-alpha*zk**2*sin_ar*r2+alpha*sin_az)*bigr
149  do idim=1,ndim
150  if(present(idir).and.idim/=idir) cycle
151  select case(idim)
152  case(1)
153  bf(ixo^s,1)=bf(ixo^s,1)+bz0(ixp1,ixp2)*((x(ixo^s,1)-xa1(ixp1))*dgdz(ixo^s)&
154  +alpha*g(ixo^s)*(x(ixo^s,2)-xa2(ixp2)))*bigr(ixo^s)
155  case(2)
156  bf(ixo^s,2)=bf(ixo^s,2)+bz0(ixp1,ixp2)*((x(ixo^s,2)-xa2(ixp2))*dgdz(ixo^s)&
157  -alpha*g(ixo^s)*(x(ixo^s,1)-xa1(ixp1)))*bigr(ixo^s)
158  case(3)
159  bf(ixo^s,3)=bf(ixo^s,3)+bz0(ixp1,ixp2)*(zk(ixo^s)*cos_ar(ixo^s)*r3(ixo^s)+alpha*&
160  zk(ixo^s)*sin_ar(ixo^s)*r2(ixo^s))
161  end select
162  end do
163  end do
164  end do
165  bf(ixo^s,:)=bf(ixo^s,:)*twopiinv
166 
167  end subroutine calc_lin_fff
168 
169  subroutine get_potential_field_potential(ixI^L,ixO^L,potential,x,zshift)
170  ! PURPOSE:
171  ! Calculation scalar potential of potential field given
172  ! Bz at photosphere (Schmidt 1964 NASSP).
173  ! NOTE: Only works for Cartesian coordinates
174  ! INPUT: x,zshift
175  ! OUTPUT: potential
177 
178  integer, intent(in) :: ixI^L, ixO^L
179  double precision, intent(in) :: x(ixI^S,1:ndim),zshift
180  double precision, intent(inout) :: potential(ixI^S)
181 
182  double precision, dimension(ixO^S) :: zk,bigr
183  integer :: ixp1,ixp2
184 
185  zk(ixo^s)=x(ixo^s,3)-xprobmin3+zshift
186  potential=0.d0
187  ! looping Bz0 pixels see equation (2)
188  !$OMP PARALLEL DO PRIVATE(bigr) REDUCTION(+:potential)
189  do ixp2=1,nx2
190  do ixp1=1,nx1
191  bigr(ixo^s)=dsqrt((x(ixo^s,1)-xa1(ixp1))**2+&
192  (x(ixo^s,2)-xa2(ixp2))**2+&
193  zk(ixo^s)**2)
194  potential(ixo^s)=potential(ixo^s)+0.5d0*bz0(ixp1,ixp2)/bigr*darea/dpi
195  end do
196  end do
197  !$OMP END PARALLEL DO
198  end subroutine get_potential_field_potential
199 
200  subroutine get_potential_field_potential_sphere(ixI^L,x,potential,nth,nph,magnetogram,theta,phi,r_sphere)
201  ! PURPOSE:
202  ! Calculation scalar potential of potential field given
203  ! Bz at photosphere (Schmidt 1964 NASSP).
204  ! NOTE: Only works for spherical coordinates
205  ! OUTPUT: potential
207  integer, intent(in) :: ixI^L,nth,nph
208  real*8, intent(in) :: x(ixi^s,1:ndim)
209  ! magnetogram Br on photosphere
210  real*8, intent(in) :: magnetogram(nth,nph)
211  ! theta and phi grid of the photospheric magnetogram
212  real*8, intent(in) :: theta(nth),phi(nph)
213  ! radius of photosphere
214  real*8, intent(in) :: r_sphere
215  real*8, intent(out) :: potential(ixi^s)
216 
217  real*8 :: area(nth),distance(ixi^s),dtheta_half,dphi,inv2pi
218  integer :: ix1,ix2
219 
220  potential=0.d0
221  ! assume uniformly discretized theta and phi
222  dtheta_half=0.5d0*(theta(2)-theta(1))
223  dphi=phi(2)-phi(1)
224  area(1:nth)=2.d0*r_sphere**2*sin(theta(1:nth))*sin(dtheta_half)*sin(dphi)
225  inv2pi=-1.d0/(2.d0*dpi)
226 
227  !$OMP PARALLEL DO PRIVATE(distance) REDUCTION(+:potential)
228  do ix2=1,nph,2
229  do ix1=1,nth,2
230  distance(ixi^s)=sqrt(x(ixi^s,1)**2+r_sphere**2-2.d0*x(ixi^s,1)*r_sphere*&
231  (sin(x(ixi^s,2))*sin(theta(ix1))*cos(phi(ix2)-x(ixi^s,3))+cos(x(ixi^s,2))*&
232  cos(theta(ix1))))
233  where(distance(ixi^s)/=0.d0)
234  distance(ixi^s)=1.d0/distance(ixi^s)
235  end where
236  potential(ixi^s)=potential(ixi^s)+inv2pi*magnetogram(ix1,ix2)*distance(ixi^s)*area(ix1)
237  end do
238  end do
239  !$OMP END PARALLEL DO
240 
242 
243  !> get potential magnetic field energy given normal B on all boundaries
244  subroutine potential_field_energy_mg(benergy)
245  ! NOTE: Solve Poisson equation of scalar potential using multigrid solver
246  ! Only works for 3D Cartesian coordinates
249  use mod_forest
250  use mod_geometry
251  real*8, intent(out) :: benergy
252  real*8 :: block_benergy(ixm^t), mype_benergy
253  real*8, allocatable :: tmp(:,:,:)
254  integer :: iigrid, igrid, idir
255  integer :: id,nc, lvl, ixI^L
256 
258  iximin^d=ixmlo^d-1;
259  iximax^d=ixmhi^d+1;
260  allocate(tmp(ixi^s))
261  mype_benergy=0.d0
262  do iigrid=1,igridstail; igrid=igrids(iigrid);
263  block=>ps(igrid)
264  ^d&dxlevel(^d)=rnode(rpdx^d_,igrid);
265  id = igrid_to_node(igrid, mype)%node%id
266  lvl= mg%boxes(id)%lvl
267  nc = mg%box_size_lvl(lvl)
268  block_benergy=0.d0
269  do idir=1,ndim
270  call gradient(mg%boxes(id)%cc({0:nc+1},mg_iphi),ixi^l,ixm^ll,idir,tmp)
271  block_benergy(ixm^t)=block_benergy(ixm^t)+tmp(ixm^t)**2
272  end do
273  mype_benergy=mype_benergy+sum(0.5d0*block_benergy(ixm^t)*block%dvolume(ixm^t))
274  end do
275 
276  call mpi_allreduce(mype_benergy,benergy,1,mpi_double_precision,mpi_sum,icomm,ierrmpi)
277 
278  end subroutine potential_field_energy_mg
279 
280  !> Solve Poisson equation of scalar potential using multigrid solver
283  use mod_comm_lib, only: mpistop
285  double precision :: max_residual, res
286  integer :: i, max_its
287 
288  mg%operator_type = mg_laplacian
289  max_its=50
290  max_residual=1.d-8
291 
292  ! Set boundary conditions
293  mg%bc(:, mg_iphi)%bc_type = mg_bc_neumann
294 
295  do i=1,2*ndim
296  mg%bc(i, mg_iphi)%boundary_cond => multigrid_bc
297  end do
298 
299  ! Solve laplacian(phi) = 0
300  do i=1,max_its
301  !call mg_fas_vcycle(mg, max_res=res)
302  call mg_fas_fmg(mg,.true.,max_res=res)
303  if(mype==0) write(*,*) 'mg iteration',i,'residual:', res
304  if(res < max_residual) exit
305  end do
306  if(res > max_residual) call mpistop("get potential multigrid: no convergence")
307 
308  end subroutine get_potential_field_potential_mg
309 
310  !> To set boundary condition on physical boundaries for mg Poisson solver
311  subroutine multigrid_bc(box, nc, iv, nb, bc_type, bc)
314  type(mg_box_t), intent(in) :: box
315  integer, intent(in) :: nc
316  integer, intent(in) :: iv !< Index of variable
317  integer, intent(in) :: nb !< number of boundary from 1 to 6 for 3D
318  integer, intent(out) :: bc_type !< Type of b.c.
319  ! mg boundary values
320  double precision, intent(out) :: bc(nc, nc)
321  ! mg coordinates of boundary layer
322  double precision :: rr(nc, nc, 3)
323  double precision :: rmina,rminb,rmaxa,rmaxb,xmina,xminb,xmaxa,xmaxb
324  ! store normal magnetic field on the boundary
325  double precision :: wbn(ixG^T)
326  double precision, allocatable :: xcoarse(:,:,:)
327  integer :: iigrid,igrid,ix^D,idir,ixbca,ixbcb,ixbcn,dlvl,wnc
328 
329  bc_type = mg_bc_neumann
330 
331  call mg_get_face_coords(box, nb, nc, rr)
332 
333  idir=(nb+1)/2
334  bc=0.d0
335  ! send normal B on boundaries from AMRVAC to mg Poisson solver
336  ! for Neumann boundary condition
337  select case(idir)
338  case(1)
339  if(mod(nb,2)==0) then
340  ixbcn=ixmhi1
341  else
342  ixbcn=ixmlo1-1
343  end if
344  rmina=rr(1,1,2)-0.5d0*box%dr(2)
345  rmaxa=rr(nc,1,2)+0.5d0*box%dr(2)
346  rminb=rr(1,1,3)-0.5d0*box%dr(3)
347  rmaxb=rr(1,nc,3)+0.5d0*box%dr(3)
348  do iigrid=1,igridstail; igrid=igrids(iigrid);
349  block=>ps(igrid)
350  if(.not.block%is_physical_boundary(nb)) cycle
351  if(stagger_grid) then
352  !wbn(ixbcn^%1ixG^T)=block%ws(ixbcn^%1ixG^T,idir)
353  wbn(ixbcn^%1ixg^t)=block%ws(ixbcn^%1ixg^t,idir)
354  else
355  wbn(ixbcn^%1ixg^t)=half*(block%w(ixbcn^%1ixg^t,iw_mag(idir))+block%w(ixbcn+1^%1ixg^t,iw_mag(idir)))
356  end if
357  xmina=block%x(1,1,1,2)-0.5d0*rnode(rpdx2_,igrid)
358  xmaxa=block%x(1,ixghi2,1,2)+0.5d0*rnode(rpdx2_,igrid)
359  xminb=block%x(1,1,1,3)-0.5d0*rnode(rpdx3_,igrid)
360  xmaxb=block%x(1,1,ixghi3,3)+0.5d0*rnode(rpdx3_,igrid)
361  if(xmina<rr(1,1,2) .and. xmaxa>rr(nc,1,2) .and. &
362  xminb<rr(1,1,3) .and. xmaxb>rr(1,nc,3)) then
363  do ix2=1,nc
364  do ix1=1,nc
365  ixbca=ceiling((rr(ix1,ix2,2)-xmina)/rnode(rpdx2_,igrid))
366  ixbcb=ceiling((rr(ix1,ix2,3)-xminb)/rnode(rpdx3_,igrid))
367  bc(ix1,ix2)=wbn(ixbcn,ixbca,ixbcb)
368  end do
369  end do
370  else if(block%x(1,ixmlo2,1,2)>rmina .and. block%x(1,ixmhi2,1,2)<rmaxa .and. &
371  block%x(1,1,ixmlo3,3)>rminb .and. block%x(1,1,ixmhi3,3)<rmaxb) then
372  dlvl=node(plevel_,igrid)-box%lvl
373  wnc=nc/2**dlvl
374  allocate(xcoarse(wnc,wnc,2))
375  do ix2=1,wnc
376  do ix1=1,wnc
377  xcoarse(ix1,ix2,1)=sum(block%x(1,(ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,1,2))/dble(2**dlvl)
378  xcoarse(ix1,ix2,2)=sum(block%x(1,1,(ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells,3))/dble(2**dlvl)
379  ixbca=ceiling((xcoarse(ix1,ix2,1)-rmina)/box%dr(2))
380  ixbcb=ceiling((xcoarse(ix1,ix2,2)-rminb)/box%dr(3))
381  bc(ixbca,ixbcb)=sum(wbn(ixbcn,(ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,&
382  (ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells))/dble(2**(2*dlvl))
383  end do
384  end do
385  deallocate(xcoarse)
386  end if
387  end do
388  case(2)
389  if(mod(nb,2)==0) then
390  ixbcn=ixmhi2
391  else
392  ixbcn=ixmlo2-1
393  end if
394  rmina=rr(1,1,1)-0.5d0*box%dr(1)
395  rmaxa=rr(nc,1,1)+0.5d0*box%dr(1)
396  rminb=rr(1,1,3)-0.5d0*box%dr(3)
397  rmaxb=rr(1,nc,3)+0.5d0*box%dr(3)
398  do iigrid=1,igridstail; igrid=igrids(iigrid);
399  block=>ps(igrid)
400  if(.not.block%is_physical_boundary(nb)) cycle
401  if(stagger_grid) then
402  wbn(ixbcn^%2ixg^t)=block%ws(ixbcn^%2ixg^t,idir)
403  else
404  wbn(ixbcn^%2ixg^t)=half*(block%w(ixbcn^%2ixg^t,iw_mag(idir))+block%w(ixbcn+1^%2ixg^t,iw_mag(idir)))
405  end if
406  xmina=block%x(1,1,1,1)-0.5d0*rnode(rpdx1_,igrid)
407  xmaxa=block%x(ixghi1,1,1,1)+0.5d0*rnode(rpdx1_,igrid)
408  xminb=block%x(1,1,1,3)-0.5d0*rnode(rpdx3_,igrid)
409  xmaxb=block%x(1,1,ixghi3,3)+0.5d0*rnode(rpdx3_,igrid)
410  if(xmina<rr(1,1,1) .and. xmaxa>rr(nc,1,1) .and. &
411  xminb<rr(1,1,3) .and. xmaxb>rr(1,nc,3)) then
412  do ix2=1,nc
413  do ix1=1,nc
414  ixbca=ceiling((rr(ix1,ix2,1)-xmina)/rnode(rpdx1_,igrid))
415  ixbcb=ceiling((rr(ix1,ix2,3)-xminb)/rnode(rpdx3_,igrid))
416  bc(ix1,ix2)=wbn(ixbca,ixbcn,ixbcb)
417  end do
418  end do
419  else if(block%x(ixmlo1,1,1,1)>rmina .and. block%x(ixmhi1,1,1,1)<rmaxa .and. &
420  block%x(1,1,ixmlo3,3)>rminb .and. block%x(1,1,ixmhi3,3)<rmaxb) then
421  dlvl=node(plevel_,igrid)-box%lvl
422  wnc=nc/2**dlvl
423  allocate(xcoarse(wnc,wnc,2))
424  do ix2=1,wnc
425  do ix1=1,wnc
426  xcoarse(ix1,ix2,1)=sum(block%x((ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,1,1,1))/dble(2**dlvl)
427  xcoarse(ix1,ix2,2)=sum(block%x(1,1,(ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells,3))/dble(2**dlvl)
428  ixbca=ceiling((xcoarse(ix1,ix2,1)-rmina)/box%dr(1))
429  ixbcb=ceiling((xcoarse(ix1,ix2,2)-rminb)/box%dr(3))
430  bc(ixbca,ixbcb)=sum(wbn((ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,ixbcn,&
431  (ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells))/dble(2**(2*dlvl))
432  end do
433  end do
434  deallocate(xcoarse)
435  end if
436  end do
437  case(3)
438  if(mod(nb,2)==0) then
439  ixbcn=ixmhi3
440  else
441  ixbcn=ixmlo3-1
442  end if
443  rmina=rr(1,1,1)-0.5d0*box%dr(1)
444  rmaxa=rr(nc,1,1)+0.5d0*box%dr(1)
445  rminb=rr(1,1,2)-0.5d0*box%dr(2)
446  rmaxb=rr(1,nc,2)+0.5d0*box%dr(2)
447  do iigrid=1,igridstail; igrid=igrids(iigrid);
448  block=>ps(igrid)
449  if(.not.block%is_physical_boundary(nb)) cycle
450  if(stagger_grid) then
451  wbn(ixbcn^%3ixg^t)=block%ws(ixbcn^%3ixg^t,idir)
452  else
453  wbn(ixbcn^%3ixg^t)=half*(block%w(ixbcn^%3ixg^t,iw_mag(idir))+block%w(ixbcn+1^%3ixg^t,iw_mag(idir)))
454  end if
455  xmina=block%x(1,1,1,1)-0.5d0*rnode(rpdx1_,igrid)
456  xmaxa=block%x(ixghi1,1,1,1)+0.5d0*rnode(rpdx1_,igrid)
457  xminb=block%x(1,1,1,2)-0.5d0*rnode(rpdx2_,igrid)
458  xmaxb=block%x(1,ixghi2,1,2)+0.5d0*rnode(rpdx2_,igrid)
459  if(xmina<rr(1,1,1) .and. xmaxa>rr(nc,1,1) .and. &
460  xminb<rr(1,1,2) .and. xmaxb>rr(1,nc,2)) then
461  do ix2=1,nc
462  do ix1=1,nc
463  ixbca=ceiling((rr(ix1,ix2,1)-xmina)/rnode(rpdx1_,igrid))
464  ixbcb=ceiling((rr(ix1,ix2,2)-xminb)/rnode(rpdx2_,igrid))
465  bc(ix1,ix2)=wbn(ixbca,ixbcb,ixbcn)
466  end do
467  end do
468  else if(block%x(ixmlo1,1,1,1)>rmina .and. block%x(ixmhi1,1,1,1)<rmaxa .and. &
469  block%x(1,ixmlo2,1,2)>rminb .and. block%x(1,ixmhi2,1,2)<rmaxb) then
470  dlvl=node(plevel_,igrid)-box%lvl
471  wnc=nc/2**dlvl
472  allocate(xcoarse(wnc,wnc,2))
473  do ix2=1,wnc
474  do ix1=1,wnc
475  xcoarse(ix1,ix2,1)=sum(block%x((ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,1,1,1))/dble(2**dlvl)
476  xcoarse(ix1,ix2,2)=sum(block%x(1,(ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells,1,2))/dble(2**dlvl)
477  ixbca=ceiling((xcoarse(ix1,ix2,1)-rmina)/box%dr(1))
478  ixbcb=ceiling((xcoarse(ix1,ix2,2)-rminb)/box%dr(2))
479  bc(ixbca,ixbcb)=sum(wbn((ix1-1)*2**dlvl+1+nghostcells:ix1*2**dlvl+nghostcells,&
480  (ix2-1)*2**dlvl+1+nghostcells:ix2*2**dlvl+nghostcells,ixbcn))/dble(2**(2*dlvl))
481  end do
482  end do
483  deallocate(xcoarse)
484  end if
485  end do
486  end select
487 
488  end subroutine multigrid_bc
489 }
490 end module mod_lfff
subroutine, public mpistop(message)
Exit MPI-AMRVAC with an error message.
Definition: mod_comm_lib.t:208
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
Module with geometry-related routines (e.g., divergence, curl)
Definition: mod_geometry.t:2
subroutine gradient(q, ixIL, ixOL, idir, gradq)
Calculate gradient of a scalar q within ixL in direction idir.
Definition: mod_geometry.t:321
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.
integer, parameter ndim
Number of spatial dimensions for grid variables.
logical stagger_grid
True for using stagger grid.
integer icomm
The MPI communicator.
integer mype
The rank of the current MPI task.
integer, parameter plevel_
integer, dimension(:), allocatable, parameter d
integer ixm
the mesh range of a physical block without ghost cells
integer ierrmpi
A global MPI error return code.
double precision, dimension(:,:), allocatable rnode
Corner coordinates.
integer nghostcells
Number of ghost cells surrounding a grid.
integer, dimension(:,:), allocatable node
double precision, dimension(ndim) dxlevel
Program to extrapolate linear force-free fields in 3D Cartesian coordinates, based on exact Green fun...
Definition: mod_lfff.t:27
double precision, save darea
Definition: mod_lfff.t:31
double precision, dimension(:,:), allocatable, save bz0
Definition: mod_lfff.t:32
subroutine get_potential_field_potential(ixIL, ixOL, potential, x, zshift)
Definition: mod_lfff.t:170
subroutine multigrid_bc(box, nc, iv, nb, bc_type, bc)
To set boundary condition on physical boundaries for mg Poisson solver.
Definition: mod_lfff.t:312
integer, save nx1
Definition: mod_lfff.t:30
subroutine get_potential_field_potential_mg()
Solve Poisson equation of scalar potential using multigrid solver.
Definition: mod_lfff.t:282
subroutine potential_field_energy_mg(benergy)
get potential magnetic field energy given normal B on all boundaries
Definition: mod_lfff.t:245
integer, save nx2
Definition: mod_lfff.t:30
subroutine init_b_fff_data(magnetogramname, qLunit, qBunit)
Definition: mod_lfff.t:38
double precision, dimension(:), allocatable, save xa2
Definition: mod_lfff.t:33
subroutine calc_lin_fff(ixIL, ixOL, Bf, x, alpha, zshift, idir)
Definition: mod_lfff.t:102
double precision, save bzmax
Definition: mod_lfff.t:31
subroutine get_potential_field_potential_sphere(ixIL, x, potential, nth, nph, magnetogram, theta, phi, r_sphere)
Definition: mod_lfff.t:201
double precision, dimension(:), allocatable, save xa1
Definition: mod_lfff.t:33
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.