! $Id$ ! !** AUTOMATIC CPARAM.INC GENERATION **************************** ! Declare (for generation of cparam.inc) the number of f array ! variables and auxiliary variables added by this module ! ! PENCILS PROVIDED x_mn; y_mn; z_mn; r_mn; r_mn1 ! PENCILS PROVIDED phix; phiy ! PENCILS PROVIDED pomx; pomy ! PENCILS PROVIDED rcyl_mn; rcyl_mn1; phi_mn ! PENCILS PROVIDED evr(3); rr(3); evth(3) ! !*************************************************************** module Grid ! use Cparam use Cdata use General, only: keep_compiler_quiet use Messages ! implicit none ! private ! public :: construct_grid public :: pencil_criteria_grid public :: pencil_interdep_grid public :: calc_pencils_grid public :: initialize_grid public :: get_grid_mn public :: box_vol public :: save_grid public :: coords_aux public :: real_to_index, inverse_grid public :: grid_bound_data public :: set_coorsys_dimmask public :: construct_serial_arrays public :: coarsegrid_interp ! public :: find_star public :: calc_bound_coeffs public :: grid_profile interface grid_profile module procedure grid_profile_0D module procedure grid_profile_1D endinterface ! interface calc_pencils_grid module procedure calc_pencils_grid_pencpar module procedure calc_pencils_grid_std endinterface calc_pencils_grid ! integer, parameter :: BOT=1, TOP=2 ! contains !*********************************************************************** !subroutine construct_grid(x,y,z,dx,dy,dz,x00,y00,z00) subroutine construct_grid(x,y,z,dx,dy,dz) ! ! Constructs a non-equidistant grid x(xi) of an equidistant grid xi with ! grid spacing dxi=1. For grid_func='linear' this is equivalent to an ! equidistant grid. ! ! dx_1 and dx_tilde are the coefficients that enter the formulae for the ! 1st and 2nd derivative: ! ! ``df/dx'' = ``df/dxi'' * dx_1 ! ``d2f/dx2'' = ``df2/dxi2'' * dx_1**2 + dx_tilde * ``df/dxi'' * dx_1 ! ! These coefficients are also very useful when adapting grid dependend stuff ! such as the timestep. A simple substitution ! 1./dx -> dx_1 ! should suffice in most cases. ! ! 25-jun-04/tobi+wolf: coded ! 3-mar-15/MR: calculation of d[xyz]2_bound added: contain twice the distances of ! three neighbouring points from the boundary point ! 9-mar-17/MR: removed unneeded use of optional parameters in calls to grid_profile ! 6-mar-18/MR: moved call construct_serial_arrays here from initialize_grid ! use Cdata, only: xprim, yprim, zprim, dx_1, dx_tilde, dy_1, dy_tilde, dz_1, dz_tilde real, dimension(mx), intent(out) :: x real, dimension(my), intent(out) :: y real, dimension(mz), intent(out) :: z real, intent(out) :: dx,dy,dz ! real :: x00,y00,z00 real :: xi1lo,xi1up,g1lo,g1up real :: xi2lo,xi2up,g2lo,g2up real :: xi3lo,xi3up,g3lo,g3up real, dimension(3,2) :: xi_step real, dimension(3,3) :: dxyz_step real :: xi1star,xi2star,xi3star,bound_prim1,bound_prim2 ! real, dimension(mx) :: g1,g1der1,g1der2,xi1,xprim2 real, dimension(my) :: g2,g2der1,g2der2,xi2,yprim2 real, dimension(mz) :: g3,g3der1,g3der2,xi3,zprim2 ! real, dimension(0:nprocx) :: xi1proc, g1proc real, dimension(0:nprocy) :: xi2proc, g2proc real, dimension(0:nprocz) :: xi3proc, g3proc ! real :: a,b,c integer :: i ! lequidist=(grid_func=='linear') ! ! Abbreviations ! x0 = xyz0(1) y0 = xyz0(2) z0 = xyz0(3) Lx = Lxyz(1) Ly = Lxyz(2) Lz = Lxyz(3) ! ! Set the lower boundary and the grid size. ! x00 = x0 y00 = y0 z00 = z0 ! dx = Lx / merge(nxgrid, max(nxgrid-1,1), lperi(1)) dy = Ly / merge(nygrid, max(nygrid-1,1), (lperi(2).or.lpole(2))) dz = Lz / merge(nzgrid, max(nzgrid-1,1), lperi(3)) ! ! Shift the lower boundary if requested, but only for periodic directions. ! if (lshift_origin(1)) x00 = x0 + 0.5 * dx if (lshift_origin(2)) y00 = y0 + 0.5 * dy if (lshift_origin(3)) z00 = z0 + 0.5 * dz ! ! Shift the lower boundary if requested, but only for periodic directions. ! Contrary to the upper case (lshift_origin) ! if (lshift_origin_lower(1)) x00 = x0 - 0.5 * dx if (lshift_origin_lower(2)) y00 = y0 - 0.5 * dy if (lshift_origin_lower(3)) z00 = z0 - 0.5 * dz ! ! Produce index arrays xi1, xi2, and xi3: ! xi = 0, 1, 2, ..., N-1 for non-periodic grid ! xi = 0.5, 1.5, 2.5, ..., N-0.5 for periodic grid ! do i=1,mx; xi1(i)=i-nghost-1+ipx*nx; enddo do i=1,my; xi2(i)=i-nghost-1+ipy*ny; enddo do i=1,mz; xi3(i)=i-nghost-1+ipz*nz; enddo ! if (lperi(1)) xi1 = xi1 + 0.5 if (lperi(2).or.lpole(2)) xi2 = xi2 + 0.5 if (lperi(3)) xi3 = xi3 + 0.5 ! ! Produce index arrays for processor boundaries, which are needed for ! particle migration (see redist_particles_bounds). The select cases ! should use these arrays to set g{2,3}proc using the grid function. ! xi1proc = (/ (real(i * nx), i = 0, nprocx) /) xi1p: if (lactive_dimension(1) .and. .not. lperi(1)) then xi1proc = xi1proc - 0.5 xi1proc(0) = 0.0 xi1proc(nprocx) = real(nxgrid - 1) endif xi1p ! xi2proc = (/ (real(i * ny), i = 0, nprocy) /) xi2p: if (lactive_dimension(2) .and. .not. (lperi(2) .or. lpole(2))) then xi2proc = xi2proc - 0.5 xi2proc(0) = 0.0 xi2proc(nprocy) = real(nygrid - 1) endif xi2p ! xi3proc = (/ (real(i * nz), i = 0, nprocz) /) xi3p: if (lactive_dimension(3) .and. .not. lperi(3)) then xi3proc = xi3proc - 0.5 xi3proc(0) = 0.0 xi3proc(nprocz) = real(nzgrid - 1) endif xi3p ! ! The following is correct for periodic and non-periodic case ! Periodic: x(xi=0) = x0 and x(xi=N) = x1 ! Non-periodic: x(xi=0) = x0 and x(xi=N-1) = x1 ! xi1lo=0.; xi1up=nxgrid-merge(0.,1.,lperi(1)) xi2lo=0.; xi2up=nygrid-merge(0.,1.,(lperi(2).or.lpole(2))) xi3lo=0.; xi3up=nzgrid-merge(0.,1.,lperi(3)) ! if (lpole(2)) xyz_star(2) = pi/2. ! ! Construct nonequidistant grid ! ! x coordinate ! if (nxgrid==1) then x = x00 + 0.5 * dx ! hopefully, we will only ever multiply by the following quantities: ! [PABourdin] We should find a way to have valid grid functions as: ! xprim = dx ; dx_1 = 1/dx ! together with a degenerated-case flag: ldegenerated_x = .true. xprim = 0. xprim2 = 0. dx_1 = 0. dx_tilde = 0. g1proc(0) = x00 g1proc(1) = x00 + Lx else ! select case (grid_func(1)) ! case ('linear','sinh') a=coeff_grid(1)*dx xi1star=find_star(a*xi1lo,a*xi1up,x00,x00+Lx,xyz_star(1),grid_func(1))/a call grid_profile(a*(xi1 -xi1star),grid_func(1),g1,g1der1,g1der2) call grid_profile(a*(xi1lo-xi1star),grid_func(1),g1lo) call grid_profile(a*(xi1up-xi1star),grid_func(1),g1up) call grid_profile(a * (xi1proc - xi1star), grid_func(1), g1proc) ! x =x00+Lx*(g1 - g1lo)/(g1up-g1lo) xprim = Lx*(g1der1*a )/(g1up-g1lo) xprim2= Lx*(g1der2*a**2)/(g1up-g1lo) g1proc = x00 + Lx * (g1proc - g1lo) / (g1up - g1lo) ! case ('cos','tanh') ! Approximately equidistant at the boundaries, linear in the middle a=pi*dx/trans_width(1) b=dxi_fact(1) xi1star = xi1lo + (xi1up-xi1lo) / (1.0 + (Lx/(xyz_star(1)-x00) - 1.0) / b) call grid_profile(a*(xi1 -xi1star),grid_func(1),g1,g1der1,g1der2,param=b) call grid_profile(a*(xi1lo-xi1star),grid_func(1),g1lo,param=b) call grid_profile(a*(xi1up-xi1star),grid_func(1),g1up,param=b) call grid_profile(a * (xi1proc - xi1star), grid_func(1), g1proc, param=b) ! x =x00+Lx*(g1 - g1lo)/(g1up-g1lo) xprim = Lx*(g1der1*a )/(g1up-g1lo) xprim2= Lx*(g1der2*a**2)/(g1up-g1lo) g1proc = x00 + Lx * (g1proc - g1lo) / (g1up - g1lo) ! case ('arsinh') ! Approximately linear in infinity, linear in the middle a=pi*dx/trans_width(1) b=dxi_fact(1) xi1star = xi1lo + (xi1up-xi1lo) / (1.0 + (Lx/(xyz_star(1)-x00) - 1.0) / b) call grid_profile(a*(xi1 -xi1star),grid_func(1),g1,g1der1,g1der2,param=b) call grid_profile(a*(xi1lo-xi1star),grid_func(1),g1lo,param=b) call grid_profile(a*(xi1up-xi1star),grid_func(1),g1up,param=b) call grid_profile(a * (xi1proc - xi1star), grid_func(1), g1proc, param=b) ! ! Slope should be 1 at the minimum: b = minval (g1der1) if (b < 1.0) then g1der1 = g1der1 + 1.0 - b g1 = g1 + (1.0 - b) * a * (xi1 - xi1star) g1lo = g1lo + (1.0 - b) * a * (xi1lo - xi1star) g1up = g1up + (1.0 - b) * a * (xi1up - xi1star) g1proc = g1proc + (1.0 - b) * a * (xi1proc - xi1star) endif ! x =x00+Lx*(g1 - g1lo)/(g1up-g1lo) xprim = Lx*(g1der1*a )/(g1up-g1lo) xprim2= Lx*(g1der2*a**2)/(g1up-g1lo) g1proc = x00 + Lx * (g1proc - g1lo) / (g1up - g1lo) ! case ('step-linear') ! xi_step(1,1)=xi_step_frac(1,1)*(nxgrid-1.0) xi_step(1,2)=xi_step_frac(1,2)*(nxgrid-1.0) dxyz_step(1,1)=(xyz_step(1,1)-x00)/(xi_step(1,1)-0.0) dxyz_step(1,2)=(xyz_step(1,2)-xyz_step(1,1))/ & (xi_step(1,2)-xi_step(1,1)) dxyz_step(1,3)=(x00+Lx-xyz_step(1,2))/(nxgrid-1.0-xi_step(1,2)) ! call grid_profile(xi1,grid_func(1),g1,g1der1,g1der2, & dxyz=dxyz_step(1,:),xistep=xi_step(1,:),delta=xi_step_width(1,:)) call grid_profile(xi1lo,grid_func(1),g1lo, & dxyz=dxyz_step(1,:),xistep=xi_step(1,:),delta=xi_step_width(1,:)) call grid_profile(xi1proc, grid_func(1), g1proc, dxyz=dxyz_step(1,:), xistep=xi_step(1,:), delta=xi_step_width(1,:)) ! x = x00 + g1-g1lo xprim = g1der1 xprim2= g1der2 g1proc = x00 + g1proc - g1lo ! case ('duct') a = pi/(max(nxgrid-1,1)) call grid_profile(a*xi1 -pi/2,grid_func(1),g1,g1der1,g1der2) call grid_profile(a*xi1lo-pi/2,grid_func(1),g1lo) call grid_profile(a*xi1up-pi/2,grid_func(1),g1up) call grid_profile(a * xi1proc - 0.5 * pi, grid_func(1), g1proc) ! x =x00+Lx*(g1-g1lo)/2 xprim = Lx*(g1der1*a )/2 xprim2= Lx*(g1der2*a**2)/2 g1proc = x00 + 0.5 * Lx * (g1proc - g1lo) ! if (lfirst_proc_x) then bound_prim1=x(l1+1)-x(l1) do i=1,nghost x(l1-i)=x(l1)-i*bound_prim1 xprim(1:l1)=bound_prim1 enddo endif if (llast_proc_x) then bound_prim2=x(l2)-x(l2-1) do i=1,nghost x(l2+i)=x(l2)+i*bound_prim2 xprim(l2:mx)=bound_prim2 enddo endif ! ! half-duct profile : like the duct above but the grid are closely spaced ! at the outer boundary. ! case ('half-duct') a =-pi/(2.*max(nxgrid-1,1)) call grid_profile(pi/2.+a*xi1,grid_func(1),g1,g1der1,g1der2) call grid_profile(pi/2.+a*xi1lo,grid_func(1),g1lo) call grid_profile(pi/2.+a*xi1up,grid_func(1),g1up) call grid_profile(0.5 * pi + a * xi1proc, grid_func(1), g1proc) ! x =x00+Lx*g1 xprim = Lx*(g1der1*a ) xprim2= Lx*(g1der2*a**2) g1proc = x00 + Lx * g1proc ! if (lfirst_proc_x) then bound_prim1=x(l1+1)-x(l1) do i=1,nghost x(l1-i)=x(l1)-i*bound_prim1 xprim(1:l1)=bound_prim1 enddo endif if (llast_proc_x) then bound_prim2=x(l2)-x(l2-1) do i=1,nghost x(l2+i)=x(l2)+i*bound_prim2 xprim(l2:mx)=bound_prim2 enddo endif ! case ('log') ! Grid distance increases logarithmically ! .i.e., grid spacing increases linearly ! d[log x] = const ==> dx = const*x a= log(xyz1(1)/xyz0(1))/(xi1up-xi1lo) b= .5*(xi1up+xi1lo-log(xyz1(1)*xyz0(1))/a) ! call grid_profile(a*(xi1-b) ,grid_func(1),g1,g1der1,g1der2) call grid_profile(a*(xi1lo-b),grid_func(1),g1lo) call grid_profile(a*(xi1up-b),grid_func(1),g1up) call grid_profile(a * (xi1proc - b), grid_func(1), g1proc) ! x =x00+Lx*(g1 - g1lo)/(g1up-g1lo) xprim = Lx*(g1der1*a )/(g1up-g1lo) xprim2= Lx*(g1der2*a**2)/(g1up-g1lo) g1proc = x00 + Lx * (g1proc - g1lo) / (g1up - g1lo) ! case ('power-law') ! Grid distance increases as a power-law set by c=coeff_grid(1) ! i.e., grid distance increases as an inverse power-law ! d[x^c] = const ==> dx = const/x^(c-1) if (lroot) then print*,"Constructing power-law grid. It will " print*,"generate a grid obeying the constrain " print*,"" print*," d[x^c] = const ==> dx = const/x^(c-1)" print*,"" print*,"where c is the value of coeff_grid." endif ! if (coeff_grid(1) == 0.) & call fatal_error('construct_grid:', 'Cannot create '//& 'a grid for a power-law of zero. Use "log" instead.') c= coeff_grid(1) a= (xyz1(1)**c-xyz0(1)**c)/(xi1up-xi1lo) b= .5*(xi1up+xi1lo-(xyz1(1)**c+xyz0(1)**c)/a) ! call grid_profile(a*(xi1-b) ,grid_func(1),g1,g1der1,g1der2,param=c) call grid_profile(a*(xi1lo-b),grid_func(1),g1lo,param=c) call grid_profile(a*(xi1up-b),grid_func(1),g1up,param=c) call grid_profile(a * (xi1proc - b), grid_func(1), g1proc, param=c) ! x =x00+Lx*(g1 - g1lo)/(g1up-g1lo) xprim = Lx*(g1der1*a )/(g1up-g1lo) xprim2= Lx*(g1der2*a**2)/(g1up-g1lo) g1proc = x00 + Lx * (g1proc - g1lo) / (g1up - g1lo) ! case ('squared') ! Grid distance increases linearily a = real(max(nxgrid, 1)) !b=-max(nxgrid,1)/10 ! 23-nov-20/ccyang: is it correct to use integer arithmetic like this? b = 0.0 call grid_profile(a*(xi1 -b),grid_func(1),g1,g1der1,g1der2) call grid_profile(a*(xi1lo-b),grid_func(1),g1lo) call grid_profile(a*(xi1up-b),grid_func(1),g1up) call grid_profile(a * (xi1proc - b), grid_func(1), g1proc) ! x =x00+Lx*(g1 - g1lo)/(g1up-g1lo) xprim = Lx*(g1der1*a )/(g1up-g1lo) xprim2= Lx*(g1der2*a**2)/(g1up-g1lo) g1proc = x00 + Lx * (g1proc - g1lo) / (g1up - g1lo) ! if (lfirst_proc_x) then bound_prim1=x(l1+1)-x(l1) do i=1,nghost x(l1-i)=x(l1)-i*bound_prim1 xprim(1:l1)=bound_prim1 enddo endif if (llast_proc_x) then bound_prim2=x(l2)-x(l2-1) do i=1,nghost x(l2+i)=x(l2)+i*bound_prim2 xprim(l2:mx)=bound_prim2 enddo endif ! case ('frozensphere') ! Just like sinh, except set dx constant below a certain radius, and ! constant for top ghost points. a=coeff_grid(1)*dx xi1star=find_star(a*xi1lo,a*xi1up,x00,x00+Lx,0.8*xyz_star(1),grid_func(1))/a call grid_profile(a*(xi1 -xi1star),grid_func(1),g1,g1der1,g1der2) call grid_profile(a*(xi1lo-xi1star),grid_func(1),g1lo) call grid_profile(a*(xi1up-xi1star),grid_func(1),g1up) call grid_profile(a * (xi1proc - xi1star), grid_func(1), g1proc) ! x =x00+Lx*(g1 - g1lo)/(g1up-g1lo) xprim = Lx*(g1der1*a )/(g1up-g1lo) xprim2= Lx*(g1der2*a**2)/(g1up-g1lo) g1proc = x00 + Lx * (g1proc - g1lo) / (g1up - g1lo) ! if (llast_proc_x) then bound_prim2=x(l2-2)-x(l2-3) do i=1,nghost+2 x(l2-2+i)=x(l2-2)+i*bound_prim2 xprim(l2-2:mx)=bound_prim2 enddo endif ! case default call fatal_error('construct_grid', & 'No such x grid function - '//trim(grid_func(1))) endselect ! ! Fitting a polynomial function to the grid function to set the xprim2 zero at the boundary, ! this makes sure that the second derivative is zero, if you choose a2 for example. ! Still in the testing phase, that why it is commented out. ! ! if ( .not.lequidist(1) ) then ! if (xprim2(l1) .ne. 0) then ! print*, 'use polynomia extension in x bottom' ! if (lfirst_proc_x) then ! a3=xprim2(l1+1)/6. ! a1=xprim(l1+1) - xprim2(l1+1)/2. ! a0=x(l1+1) - xprim(l1+1) + xprim2(l1+1)/3. ! x(l1)=a0 ! xprim(l1) = a1 ! xprim2(l1) = 0. ! endif ! endif ! if (xprim2(l2) .ne. 0) then ! print*, 'use polynomia extension in x top' ! if (llast_proc_x) then ! a3=-xprim2(l2-1)/6. ! a1=xprim(l2-1) + xprim2(l2-1)/2. ! a0=x(l2-1) +xprim(l2-1) + xprim2(l2-1)/3. ! x(l2)=a0 ! xprim(l2) = a1 ! xprim2(l2) = 0 ! endif ! endif ! endif ! dx2=xprim**2 dx_1=1./xprim dx_tilde=-xprim2/dx2 ! endif ! ! y coordinate ! if (nygrid==1) then y = y00 + 0.5 * dy ! hopefully, we will only ever multiply by the following quantities: ! [PABourdin] We should find a way to have valid grid functions as: ! xprim = dx ; dx_1 = 1/dx ! together with a degenerated-case flag: ldegenerated_x = .true. yprim = 0. yprim2 = 0. dy_1 = 0. dy_tilde = 0. g2proc(0) = y00 g2proc(1) = y00 + Ly else ! select case (grid_func(2)) ! case ('linear','sinh') ! a=coeff_grid(2)*dy xi2star=find_star(a*xi2lo,a*xi2up,y00,y00+Ly,xyz_star(2),grid_func(2))/a call grid_profile(a*(xi2 -xi2star),grid_func(2),g2,g2der1,g2der2) call grid_profile(a*(xi2lo-xi2star),grid_func(2),g2lo) call grid_profile(a*(xi2up-xi2star),grid_func(2),g2up) call grid_profile(a * (xi2proc - xi2star), grid_func(2), g2proc) ! y =y00+Ly*(g2 - g2lo)/(g2up-g2lo) yprim = Ly*(g2der1*a )/(g2up-g2lo) yprim2= Ly*(g2der2*a**2)/(g2up-g2lo) g2proc = y00 + Ly * (g2proc - g2lo) / (g2up - g2lo) ! case ('cos','tanh') ! Approximately equidistant at the boundaries, linear in the middle a=pi*dy/trans_width(2) b=dxi_fact(2) xi2star = xi2lo + (xi2up-xi2lo) / (1.0 + (Ly/(xyz_star(2)-y00) - 1.0) / b) call grid_profile(a*(xi2 -xi2star),grid_func(2),g2,g2der1,g2der2,param=b) call grid_profile(a*(xi2lo-xi2star),grid_func(2),g2lo,param=b) call grid_profile(a*(xi2up-xi2star),grid_func(2),g2up,param=b) call grid_profile(a * (xi2proc - xi2star), grid_func(2), g2proc, param=b) ! y =y00+Ly*(g2 - g2lo)/(g2up-g2lo) yprim = Ly*(g2der1*a )/(g2up-g2lo) yprim2= Ly*(g2der2*a**2)/(g2up-g2lo) g2proc = y00 + Ly * (g2proc - g2lo) / (g2up - g2lo) ! case ('arsinh') ! Approximately linear in infinity, linear in the middle a=pi*dy/trans_width(2) b=dxi_fact(2) xi2star = xi2lo + (xi2up-xi2lo) / (1.0 + (Ly/(xyz_star(2)-y00) - 1.0) / b) call grid_profile(a*(xi2 -xi2star),grid_func(2),g2,g2der1,g2der2,param=b) call grid_profile(a*(xi2lo-xi2star),grid_func(2),g2lo,param=b) call grid_profile(a*(xi2up-xi2star),grid_func(2),g2up,param=b) call grid_profile(a * (xi2proc - xi2star), grid_func(2), g2proc, param=b) ! ! Slope should be 1 at the minimum: b = minval (g2der1) if (b < 1.0) then g2der1 = g2der1 + 1.0 - b g2 = g2 + (1.0 - b) * a * (xi2 - xi2star) g2lo = g2lo + (1.0 - b) * a * (xi2lo - xi2star) g2up = g2up + (1.0 - b) * a * (xi2up - xi2star) g2proc = g2proc + (1.0 - b) * a * (xi2proc - xi2star) endif ! y =y00+Ly*(g2 - g2lo)/(g2up-g2lo) yprim = Ly*(g2der1*a )/(g2up-g2lo) yprim2= Ly*(g2der2*a**2)/(g2up-g2lo) g2proc = y00 + Ly * (g2proc - g2lo) / (g2up - g2lo) ! case ('duct') ! a = pi/max(nygrid-1, 1) call grid_profile(a*xi2 -pi/2,grid_func(2),g2,g2der1,g2der2) call grid_profile(a*xi2lo-pi/2,grid_func(2),g2lo) call grid_profile(a*xi2up-pi/2,grid_func(2),g2up) call grid_profile(a * xi2proc - 0.5 * pi, grid_func(2), g2proc) ! y =y00+Ly*(g2-g2lo)/2 yprim = Ly*(g2der1*a )/2 yprim2= Ly*(g2der2*a**2)/2 g2proc = y00 + 0.5 * Ly * (g2proc - g2lo) ! if (lfirst_proc_y) then bound_prim1=y(m1+1)-y(m1) do i=1,nghost y(m1-i)=y(m1)-i*bound_prim1 yprim(1:m1)=bound_prim1 enddo endif if (llast_proc_y) then bound_prim2=y(m2)-y(m2-1) do i=1,nghost y(m2+i)=y(m2)+i*bound_prim2 yprim(m2:my)=bound_prim2 enddo endif ! case ('step-linear') ! xi_step(2,1)=xi_step_frac(2,1)*merge(nygrid+0.,nygrid-1.,lpole(2)) xi_step(2,2)=xi_step_frac(2,2)*merge(nygrid+0.,nygrid-1.,lpole(2)) dxyz_step(2,1)=(xyz_step(2,1)-y00)/(xi_step(2,1)-0.0) dxyz_step(2,2)=(xyz_step(2,2)-xyz_step(2,1))/ & (xi_step(2,2)-xi_step(2,1)) dxyz_step(2,3)=(y00+Ly-xyz_step(2,2))/& (merge(nygrid+0.,nygrid-1.,lpole(2))-xi_step(2,2)) ! if (lpole(2)) then do i=1,my xi2(i)=0.5-nghost+0.5+(ipy*ny+i-1)*& (nygrid+2.*(nghost-0.5)-1)/(nygrid+2.*nghost-1) enddo xi2proc = (/ (real(1 - nghost) + (real(nghost + i * ny) - 0.5) * real(mygrid - 2) & / real(mygrid - 1), i = 0, nprocy) /) endif ! call grid_profile(xi2,grid_func(2),g2,g2der1,g2der2,dxyz= & dxyz_step(2,:),xistep=xi_step(2,:),delta=xi_step_width(2,:)) call grid_profile(xi2lo,grid_func(2),g2lo,dxyz= & dxyz_step(2,:),xistep=xi_step(2,:),delta=xi_step_width(2,:)) call grid_profile(xi2proc, grid_func(2), g2proc, dxyz=dxyz_step(2,:), xistep=xi_step(2,:), delta=xi_step_width(2,:)) ! if (lpole(2)) then call grid_profile(xi2up,grid_func(2),g2up,dxyz= & dxyz_step(2,:),xistep=xi_step(2,:),delta=xi_step_width(2,:)) y = y00 + g2 -0.5*(g2lo+g2up-pi) g2proc = y00 + g2proc - 0.5 * (g2lo + g2up - pi) else !if(iproc<2) print*, 'iproc, y00, g2-g2lo=', iproc, y00, g2-g2lo y = y00 + g2-g2lo g2proc = y00 + g2proc - g2lo endif ! yprim = g2der1 yprim2= g2der2 ! case ('trough') ! Grid distance is almost equidistant at boundaries and then decreases ! at the middle a=0.05 b=40 c=160-40 ! call grid_profile(xi2, grid_func(2),g2,g2der1,g2der2,param=a,xistep=(/b,c/),delta=(/0.5,0.1/)) call grid_profile(xi2lo,grid_func(2),g2lo,param=a,xistep=(/b,c/),delta=(/0.5,0.1/)) call grid_profile(xi2up,grid_func(2),g2up,param=a,xistep=(/b,c/),delta=(/0.5,0.1/)) call grid_profile(xi2proc, grid_func(2), g2proc, param=a, xistep=(/b,c/), delta=(/0.5,0.1/)) ! y =y00+Ly*(g2 - g2lo)/(g2up-g2lo) yprim = Ly*g2der1/(g2up-g2lo) yprim2= Ly*g2der2/(g2up-g2lo) g2proc = y00 + Ly * (g2proc - g2lo) / (g2up - g2lo) ! case default call fatal_error('construct_grid', & 'No such y grid function - '//trim(grid_func(2))) ! endselect ! ! Added parts for spherical coordinates and cylindrical coordinates. ! From now on dy = d\theta but dy_1 = 1/rd\theta and similarly for \phi. ! corresponding r and rsin\theta factors for equ.f90 (where CFL timesteps ! are estimated) are removed. ! if (lpole(2)) then !apply grid symmetry across the poles if (lfirst_proc_y) then y( 1:nghost) = -y( m1i:m1:-1) yprim( 1:nghost) = yprim( m1i:m1:-1) yprim2(1:nghost) = -yprim2(m1i:m1:-1) endif if (llast_proc_y) then y( m2+1:m2+nghost) = 2.*pi-y( m2:m2i:-1) yprim( m2+1:m2+nghost) = yprim( m2:m2i:-1) yprim2(m2+1:m2+nghost) = -yprim2(m2:m2i:-1) endif endif ! dy2=yprim**2 dy_1=1./yprim dy_tilde=-yprim2/dy2 ! endif ! ! z coordinate ! if (nzgrid==1) then z = z00 + 0.5 * dz ! hopefully, we will only ever multiply by the following quantities: ! [PABourdin] We should find a way to have valid grid functions as: ! xprim = dx ; dx_1 = 1/dx ! together with a degenerated-case flag: ldegenerated_x = .true. zprim = 0. zprim2 = 0. dz_1 = 0. dz_tilde = 0. g3proc(0) = z00 g3proc(1) = z00 + Lz else ! select case (grid_func(3)) ! case ('linear','sinh') ! a=coeff_grid(3)*dz xi3star=find_star(a*xi3lo,a*xi3up,z00,z00+Lz,xyz_star(3),grid_func(3))/a call grid_profile(a*(xi3 -xi3star),grid_func(3),g3,g3der1,g3der2) call grid_profile(a*(xi3lo-xi3star),grid_func(3),g3lo) call grid_profile(a*(xi3up-xi3star),grid_func(3),g3up) call grid_profile(a * (xi3proc - xi3star), grid_func(3), g3proc) ! z =z00+Lz*(g3 - g3lo)/(g3up-g3lo) zprim = Lz*(g3der1*a )/(g3up-g3lo) zprim2= Lz*(g3der2*a**2)/(g3up-g3lo) g3proc = z00 + Lz * (g3proc - g3lo) / (g3up - g3lo) ! case ('cos','tanh') ! Approximately equidistant at the boundaries, linear in the middle a=pi*dz/trans_width(3) b=dxi_fact(3) xi3star = xi3lo + (xi3up-xi3lo) / (1.0 + (Lz/(xyz_star(3)-z00) - 1.0) / b) call grid_profile(a*(xi3 -xi3star),grid_func(3),g3,g3der1,g3der2,param=b) call grid_profile(a*(xi3lo-xi3star),grid_func(3),g3lo,param=b) call grid_profile(a*(xi3up-xi3star),grid_func(3),g3up,param=b) call grid_profile(a * (xi3proc - xi3star), grid_func(3), g3proc, param=b) ! z =z00+Lz*(g3 - g3lo)/(g3up-g3lo) zprim = Lz*(g3der1*a )/(g3up-g3lo) zprim2= Lz*(g3der2*a**2)/(g3up-g3lo) g3proc = z00 + Lz * (g3proc - g3lo) / (g3up - g3lo) ! case ('arsinh') ! Approximately linear in infinity, linear in the middle a=pi*dz/trans_width(3) b=dxi_fact(3) xi3star = xi3lo + (xi3up-xi3lo) / (1.0 + (Lz/(xyz_star(3)-z00) - 1.0) / b) call grid_profile(a*(xi3 -xi3star),grid_func(3),g3,g3der1,g3der2,param=b) call grid_profile(a*(xi3lo-xi3star),grid_func(3),g3lo,param=b) call grid_profile(a*(xi3up-xi3star),grid_func(3),g3up,param=b) call grid_profile(a * (xi3proc - xi3star), grid_func(3), g3proc, param=b) ! ! Slope should be 1 at the minimum: b = minval (g3der1) if ((grid_func(3) == 'arsinh') .and. (b < 1.0)) then g3der1 = g3der1 + 1.0 - b g3 = g3 + (1.0 - b) * a * (xi3 - xi3star) g3lo = g3lo + (1.0 - b) * a * (xi3lo - xi3star) g3up = g3up + (1.0 - b) * a * (xi3up - xi3star) g3proc = g3proc + (1.0 - b) * a * (xi3proc - xi3star) endif ! z =z00+Lz*(g3 - g3lo)/(g3up-g3lo) zprim = Lz*(g3der1*a )/(g3up-g3lo) zprim2= Lz*(g3der2*a**2)/(g3up-g3lo) g3proc = z00 + Lz * (g3proc - g3lo) / (g3up - g3lo) ! case ('step-linear') ! xi_step(3,1)=xi_step_frac(3,1)*(nzgrid-1.0) xi_step(3,2)=xi_step_frac(3,2)*(nzgrid-1.0) dxyz_step(3,1)=(xyz_step(3,1)-z00)/(xi_step(3,1)-0.0) dxyz_step(3,2)=(xyz_step(3,2)-xyz_step(3,1))/ & (xi_step(3,2)-xi_step(3,1)) dxyz_step(3,3)=(z00+Lz-xyz_step(3,2))/(nzgrid-1.0-xi_step(3,2)) ! call grid_profile(xi3,grid_func(3),g3,g3der1,g3der2, & dxyz=dxyz_step(3,:),xistep=xi_step(3,:),delta=xi_step_width(3,:)) call grid_profile(xi3lo,grid_func(3),g3lo, & dxyz=dxyz_step(3,:),xistep=xi_step(3,:),delta=xi_step_width(3,:)) call grid_profile(xi3proc, grid_func(3), g3proc, dxyz=dxyz_step(3,:), xistep=xi_step(3,:), delta=xi_step_width(3,:)) ! z = z00 + g3-g3lo zprim = g3der1 zprim2= g3der2 g3proc = z00 + g3proc - g3lo ! case ('linear+log') ! Grid distance is equidistant first and then increases logarithmically ! .i.e., grid spacing increases linearly ! d[log z] = const ==> dz = const*z a=4.0/float(nzgrid) b=200 c=5 ! call grid_profile(a*(xi3-b) ,grid_func(3),g3,g3der1,g3der2,param=c) call grid_profile(a*(xi3lo-b),grid_func(3),g3lo,param=c) call grid_profile(a*(xi3up-b),grid_func(3),g3up,param=c) call grid_profile(a * (xi3proc - b), grid_func(3), g3proc, param=c) ! z =z00+Lz*(g3 - g3lo)/(g3up-g3lo) zprim = Lz*(g3der1*a )/(g3up-g3lo) zprim2= Lz*(g3der2*a**2)/(g3up-g3lo) g3proc = z00 + Lz * (g3proc - g3lo) / (g3up - g3lo) ! case ('squared') ! Grid distance increases linearily a=max(nzgrid,1) !b=-max(nzgrid,1)/10 ! 23-nov-20/ccyang: is it correct to use integer arithmetic like this? b = 0.0 call grid_profile(a*(xi3 -b),grid_func(3),g3,g3der1,g3der2) call grid_profile(a*(xi3lo-b),grid_func(3),g3lo) call grid_profile(a*(xi3up-b),grid_func(3),g3up) call grid_profile(a * (xi3proc - b), grid_func(3), g3proc) ! z =z00+Lz*(g3 - g3lo)/(g3up-g3lo) zprim = Lz*(g3der1*a )/(g3up-g3lo) zprim2= Lz*(g3der2*a**2)/(g3up-g3lo) g3proc = z00 + Lz * (g3proc - g3lo) / (g3up - g3lo) ! if (lfirst_proc_z) then bound_prim1=z(n1+1)-z(n1) do i=1,nghost z(n1-i)=z(n1)-i*bound_prim1 zprim(1:n1)=bound_prim1 enddo endif if (llast_proc_z) then bound_prim2=z(n2)-z(n2-1) do i=1,nghost z(n2+i)=z(n2)+i*bound_prim2 zprim(n2:mz)=bound_prim2 enddo endif ! case ('trough') ! Grid distance is almost equidistant at boundaries and then decreases ! at the middle a=0.05 b=47 c=192-47 ! call grid_profile(xi3 ,grid_func(3),g3,g3der1,g3der2,param=a,xistep=(/b,c/),delta=(/0.5,0.1/)) call grid_profile(xi3lo,grid_func(3),g3lo,param=a,xistep=(/b,c/),delta=(/0.5,0.1/)) call grid_profile(xi3up,grid_func(3),g3up,param=a,xistep=(/b,c/),delta=(/0.5,0.1/)) call grid_profile(xi3proc, grid_func(3), g3proc, param=a, xistep=(/b,c/), delta=(/0.5,0.1/)) ! z =z00+Lz*(g3 - g3lo)/(g3up-g3lo) zprim = Lz*g3der1/(g3up-g3lo) zprim2= Lz*g3der2/(g3up-g3lo) g3proc = z00 + Lz * (g3proc - g3lo) / (g3up - g3lo) ! case default call fatal_error('construct_grid', & 'No such z grid function - '//trim(grid_func(3))) endselect ! dz2=zprim**2 dz_1=1./zprim dz_tilde=-zprim2/dz2 ! endif ! ! Record the processor boundary. ! procx_bounds = g1proc procy_bounds = g2proc procz_bounds = g3proc ! call grid_bound_data ! ! Fargo (orbital advection acceleration) is implemented for polar coordinates only. ! Die otherwise. ! if (lfargo_advection.and.coord_system=='cartesian') then if (lroot) then print*,"" print*,"Fargo advection is only implemented for" print*,"polar coordinates. Switch" print*," coord_system='cylindric' or 'spherical'" print*,"in init_pars of start.in if you" print*,"want to use the fargo algorithm" print*,"" endif call fatal_error("construct_grid","") endif ! ! Set equator if possible. ! if (abs(xyz0(2)+xyz1(2))-pi<1e-3) lequatory=.true. ! ! Set the serial grid arrays, that contain the coordinate values ! from all processors. ! !!! call construct_serial_arrays !!! creates trouble ! endsubroutine construct_grid !*********************************************************************** subroutine set_coorsys_dimmask ! ! Sets switches for the different coordinate systems. ! ! 18-dec-15/MR: outsourced from initialize_grid ! 16-jan-17/MR: added initialization of dimensionality mask. ! lcartesian_coords=.false. lspherical_coords=.false. lcylindrical_coords=.false. lpipe_coords=.false. ! if (coord_system=='cartesian') then lcartesian_coords=.true. ! ! Introduce new names (spherical_coords), in addition to the old ones. ! elseif ( coord_system=='spherical' & .or.coord_system=='spherical_coords') then lspherical_coords=.true. ! ! Introduce new names (cylindrical_coords), in addition to the old ones. ! elseif ( coord_system=='cylindric' & .or.coord_system=='cylindrical_coords') then lcylindrical_coords=.true. else if (coord_system=='pipeflows') then lpipe_coords=.true. elseif (coord_system=='Lobachevskii') then call fatal_error('set_coorsys_dimmask', & 'Lobachevskii ccordinates not implemented') endif ! ! Initialize dimensionality mask. ! if (nxgrid==1) then if (nygrid==1) then dim_mask(1)=3 else dim_mask(1:2)=(/2,3/) endif else if (nygrid==1) dim_mask(1:2)=(/1,3/) endif endsubroutine set_coorsys_dimmask !*********************************************************************** subroutine initialize_grid ! ! Coordinate-related issues: nonuniform meshes, different coordinate systems ! ! 20-jul-10/wlad: moved here from register ! 3-mar-14/MR: outsourced calculation of box_volume into box_vol ! 29-sep-14/MR: outsourced calculation of auxiliary quantities for curvilinear ! coordinates into coords_aux; set coordinate switches at the ! beginning ! 9-jun-15/MR: calculation of Area_* added ! 18-dec-15/MR: outsourced setting of switches to set_coords_switches; added ! initialization of Yin-Yang grid; added dimensionality mask: ! lists the indices of the non-degenerate directions in the first ! dimensionality elements of dim_mask ! 10-oct-17/MR: avoided communication in calculation of r_int and r_ext ! 10-jan-17/MR: moved call construct_serial_arrays to beginning ! use Sub, only: remove_prof use Mpicomm use IO, only: lcollective_IO use General, only: indgen, itoa, find_proc ! real :: fact, dxmin_x, dxmin_y, dxmin_z, dxmax_x, dxmax_y, dxmax_z integer :: xj,yj,zj,itheta,nphi,na,ne,mm,nn,nni integer, dimension(3) :: idzl,idzr ! call construct_serial_arrays ! ! For curvilinear coordinate systems, calculate auxiliary quantities as, e.g., for spherical coordinates 1/r, cot(theta)/r, etc. ! call coords_aux(x,y,z) ! ! determine global minimum and maximum of grid spacing in any direction ! if (lequidist(1) .or. nxgrid <= 1) then dxmin_x = dx dxmax_x = dx else dxmin_x = minval(xprim(l1:l2)) dxmax_x = maxval(xprim(l1:l2)) endif ! if (lequidist(2) .or. nygrid <= 1) then if (lspherical_coords.or.lcylindrical_coords) then dxmin_y = dy*minval(x(l1:l2)) dxmax_y = dy*maxval(x(l1:l2)) else dxmin_y = dy ! possibly incorrect for pipe coordinates dxmax_y = dy endif else dxmin_y = minval(yprim(m1:m2)) dxmax_y = maxval(yprim(m1:m2)) endif ! if (lequidist(3) .or. nzgrid <= 1) then if (lspherical_coords) then dxmin_z = dz*minval(x(l1:l2))*minval(sinth(m1:m2)) dxmax_z = dz*maxval(x(l1:l2))*maxval(sinth(m1:m2)) else dxmin_z = dz ! possibly incorrect for pipe coordinates dxmax_z = dz endif else dxmin_z = minval(zprim(n1:n2)) dxmax_z = maxval(zprim(n1:n2)) endif ! ! Find minimum/maximum grid spacing. Note that ! minval( (/dxmin_x,dxmin_y,dxmin_z/), MASK=((/nxgrid,nygrid,nzgrid/) > 1) ) ! will be undefined if all n[xyz]grid==1, so we have to add the fourth ! component with a test that is always true ! dxmin = minval( (/dxmin_x, dxmin_y, dxmin_z, huge(dx)/), & MASK=((/nxgrid, nygrid, nzgrid, 2/) > 1) ) call mpiallreduce_min(dxmin,dxmin_x) dxmin=dxmin_x ! if (dxmin == 0) & call fatal_error ("initialize_grid", "check Lx,Ly,Lz: is one of them 0?", .true.) ! dxmax = maxval( (/dxmax_x, dxmax_y, dxmax_z, epsilon(dx)/), & MASK=((/nxgrid, nygrid, nzgrid, 2/) > 1) ) call mpiallreduce_max(dxmax,dxmax_x) dxmax=dxmax_x ! ! Grid spacing. For non-equidistant grid or non-Cartesian coordinates ! the grid spacing is calculated in the (m,n) loop. ! if (lcartesian_coords .and. all(lequidist)) then ! ! FAG replaced old_cdtv flag with more general coordinate independent lmaximal ! if (old_cdtv) then ! dxyz_2 = max(dx_1(l1:l2)**2,dy_1(m1)**2,dz_1(n1)**2) ! else dline_1(:,1)=dx_1(l1:l2) dline_1(:,2)=dy_1(m1) dline_1(:,3)=dz_1(n1) ! if (lmaximal_cdtv) then dxyz_2 = max(dline_1(:,1)**2, dline_1(:,2)**2, dline_1(:,3)**2) dxyz_4 = max(dline_1(:,1)**4, dline_1(:,2)**4, dline_1(:,3)**4) dxyz_6 = max(dline_1(:,1)**6, dline_1(:,2)**6, dline_1(:,3)**6) else dxyz_2 = dline_1(:,1)**2 + dline_1(:,2)**2 + dline_1(:,3)**2 dxyz_4 = dline_1(:,1)**4 + dline_1(:,2)**4 + dline_1(:,3)**4 dxyz_6 = dline_1(:,1)**6 + dline_1(:,2)**6 + dline_1(:,3)**6 endif ! dxyz_2 = dline_1(:,1)**2+dline_1(:,2)**2+dline_1(:,3)**2 ! dxyz_4 = dline_1(:,1)**4+dline_1(:,2)**4+dline_1(:,3)**4 ! dxyz_6 = dline_1(:,1)**6+dline_1(:,2)**6+dline_1(:,3)**6 ! endif ! ! Fill pencil with maximum gridspacing. Will be overwritten ! during the mn loop in the non-equidistant or non-Cartesian case. ! dxmax_pencil = dxmax dxmin_pencil = dxmin ! endif ! ! Box volume ! call box_vol ! ! Initialize Yin-Yang grid. ! if (lyinyang) call yyinit ! ! Volume element and area of coordinate surfaces. ! Note that in the area a factor depending only on the coordinate x_i which defines the surface by x_i=const. is dropped. ! Area_xy=1.; Area_yz=1.; Area_xz=1. ! if (lcartesian_coords) then ! ! x-extent ! if (nxgrid/=1) then dVol_x=xprim Area_xy=Area_xy*Lxyz(1) Area_xz=Area_xz*Lxyz(1) else dVol_x=1. endif ! ! y-extent ! if (nygrid/=1) then dVol_y=yprim Area_xy=Area_xy*Lxyz(2) Area_yz=Area_yz*Lxyz(2) else dVol_y=1. endif ! ! z-extent ! if (nzgrid/=1) then dVol_z=zprim Area_yz=Area_yz*Lxyz(3) Area_xz=Area_xz*Lxyz(3) else dVol_z=1. endif ! ! single value volume element dVol applicable only to Cartesian equidistant grid ! if (all(lequidist)) dVol=dVol_x(l1:l2)*dVol_y(m1)*dVol_z(n1) ! ! Spherical coordinate system ! elseif (lspherical_coords) then ! ! Volume element - it is wrong for spherical, since ! sinth also changes with y-position. ! ! WL: Is this comment above still relevant? ! ! Split up volume differential as (dr) * (r*dtheta) * (r*sinth*dphi) ! and assume that sinth=1 if there is no theta extent. ! This should always give a volume of 4pi/3*(r2^3-r1^3) for constant integrand ! r extent: ! if (nxgrid/=1) then dVol_x=x**2*xprim Area_xy=Area_xy*1./2.*(xyz1(1)**2-xyz0(1)**2) Area_xz=Area_xz*1./2.*(xyz1(1)**2-xyz0(1)**2) else ! dVol_x=1./3.*(xyz1(1)**3-xyz0(1)**3) !??? ! Area_xy=Area_xy*1./2.*(xyz1(1)**2-xyz0(1)**2) !??? ! Area_xz=Area_xz*1./2.*(xyz1(1)**2-xyz0(1)**2) dVol_x=1./3. Area_xy=Area_xy*1./2. Area_xz=Area_xz*1./2. endif ! ! Theta extent (if non-radially symmetric) ! if (nygrid/=1) then dVol_y=sinth*yprim Area_xy=Area_xy*Lxyz(2) Area_yz=Area_yz*(cos(xyz0(2))-cos(xyz1(2))) else dVol_y=2. Area_xy=Area_xy*pi Area_yz=Area_yz*2. endif ! ! phi extent (if non-axisymmetry) ! if (nzgrid/=1) then dVol_z=zprim Area_xz=Area_xz*Lxyz(3) Area_yz=Area_yz*Lxyz(3) else dVol_z=2.*pi Area_xz=Area_xz*2.*pi Area_yz=Area_yz*2.*pi endif ! ! weighted coordinates for integration purposes ! Need to modify for 2-D and 1-D cases! !AB: for now, allow only if nxgrid>1. Dhruba, please check ! r2_weight=x(l1:l2)**2 sinth_weight=sinth if (nxgrid>1) then do itheta=1,nygrid sinth_weight_across_proc(itheta)=sin(xyz0(2)+dy*itheta) enddo else sinth_weight_across_proc=1. endif ! ! Calculate the volume of the box, for spherical coordinates. ! MR: Why needed? Box_volume should be the same (but more accurate). ! nVol, nVol1 never used! ! nVol=0. do xj=l1,l2 do yj=m1,m2 do zj=n1,n2 nVol=nVol+x(xj)*x(xj)*sinth(yj) enddo enddo enddo nVol1=1./nVol ! ! Trapezoidal rule ! if (lfirst_proc_x) r2_weight( 1)=.5*r2_weight( 1) if (llast_proc_x ) r2_weight(nx)=.5*r2_weight(nx) ! if (lfirst_proc_y) sinth_weight(m1)=.5*sinth_weight(m1) if (llast_proc_y ) sinth_weight(m2)=.5*sinth_weight(m2) sinth_weight_across_proc(1 )=0.5*sinth_weight_across_proc(1) sinth_weight_across_proc(nygrid)=0.5*sinth_weight_across_proc(nygrid) ! ! End of coord_system=='spherical_coords' query. ! elseif (lcylindrical_coords) then ! ! Volume element. ! if (nxgrid/=1) then dVol_x=x*xprim Area_xy=Area_xy*1./2.*(xyz1(1)**2-xyz0(1)**2) Area_xz=Area_xz*Lxyz(1) else ! dVol_x=1./2.*(xyz1(1)**2-xyz0(1)**2) !??? dVol_x=1./2. Area_xy=Area_xy*dVol_x(1) Area_xz=Area_xz*Lxyz(1) endif ! ! theta extent (non-cylindrically symmetric) ! if (nygrid/=1) then dVol_y=yprim Area_xy=Area_xy*Lxyz(2) Area_yz=Area_yz*Lxyz(2) else dVol_y=2.*pi Area_xy=Area_xy*2.*pi Area_yz=Area_yz*2.*pi endif ! ! z extent (vertically extended) ! if (nzgrid/=1) then dVol_z=zprim Area_xz=Area_xz*Lxyz(3) Area_yz=Area_yz*Lxyz(3) else dVol_z=1. endif ! ! Trapezoidal rule ! rcyl_weight=rcyl_mn if (lfirst_proc_x) rcyl_weight( 1)=.5*rcyl_weight( 1) if (llast_proc_x ) rcyl_weight(nx)=.5*rcyl_weight(nx) ! ! Pipe coordinates (for hydraulic applications) ! elseif (lpipe_coords) then ! ! Compute profile function. ! select case (pipe_func) ! case ('error_function') call warning('initialize_grid','calculation of Area_* not implemented') fact=.5/CrossSec_w**2 glnCrossSec=glnCrossSec0*(-exp(-fact*(x(l1:l2)-CrossSec_x1)**2) & +exp(-fact*(x(l1:l2)-CrossSec_x2)**2)) ! case default call fatal_error('initialize_grid', & 'no such pipe function - '//trim(pipe_func)) endselect ! ! Lobachevskii space ! ! else ! ! Stop if no existing coordinate system is specified ! else call fatal_error('initialize_grid', & 'no such coordinate system: '//coord_system) endif ! ! Inverse volume elements ! dVol1_x = dVol_x dVol1_y = dVol_y dVol1_z = dVol_z ! ! Avoid 1/0 at axis or in origin. ! if (lspherical_coords) then if (dVol1_x(l1) == 0.) dVol1_x(l1)=.5*x(l1+1)**2*xprim(l1+1) if (dVol1_y(m1) == 0.) dVol1_y(m1)=.5*sinth(m1+1)*yprim(m1+1) if (dVol1_y(m2) == 0.) dVol1_y(m2)=.5*sinth(m2-1)*yprim(m2-1) elseif (lcylindrical_coords) then if (dVol1_x(l1) == 0.) dVol1_x(l1)=.5*x(l1+1)*xprim(l1+1) endif dVol1_x = 1./dVol1_x dVol1_y = 1./dVol1_y dVol1_z = 1./dVol1_z ! if (lspherical_coords.or.lcylindrical_coords) then ! ! Define inner and outer radii for non-cartesian coords. ! If the user did not specify them yet (in start.in), ! these are the first point of the first x-processor, ! and the last point of the last x-processor. ! r_int and r_ext can be taken from the domain's x-extent as periodic BC can't ! appear for r. ! if (r_int == 0) r_int=xyz0(1) if (r_ext ==impossible) r_ext=xyz1(1) if (lroot.and.ip<14) print*,'initialize_grid, r_int,r_ext=',r_int,r_ext endif ! ! For a non-periodic mesh, multiply boundary points by 1/2. ! Do it for each direction in turn. ! If a direction has no extent, it is automatically periodic ! and the corresponding step is therefore not called. ! if (.not.lperi(1)) then if (lfirst_proc_x) dVol_x(l1)=.5*dVol_x(l1) if (llast_proc_x ) dVol_x(l2)=.5*dVol_x(l2) endif ! if (.not.lperi(2)) then if (lfirst_proc_y) dVol_y(m1)=.5*dVol_y(m1) if (llast_proc_y ) dVol_y(m2)=.5*dVol_y(m2) endif ! if (.not.lperi(3)) then if (lfirst_proc_z) dVol_z(n1)=.5*dVol_z(n1) if (llast_proc_z ) dVol_z(n2)=.5*dVol_z(n2) endif ! ! Check max possible ncoarse. ! if (nzgrid==1) then ! switch off coarsening if no z-extent ncoarse=0 elseif (ncoarse>nz/nghost) then ncoarse=nz/nghost call warning('initialize_grid','there are jumped-over processors due to grid coarsening'// & ' -> ncoarse reduced to floor(nz/nghost)='//trim(itoa(ncoarse))) endif ! ! Set lcoarse for coarsening grid near poles. ! lcoarse=lspherical_coords .and. lpole(2) .and. ncoarse>1 .and. nprocy>1 if (.not.lcoarse) ncoarse=1 if (lcoarse.and.(lfirst_proc_y.or.llast_proc_y)) then !MR: TB generalized to more than two procs, in which coarsening happens. if (lfirst_proc_y) then mexts=(/m1,m1+ncoarse-2/) elseif (llast_proc_y) then mexts=(/m2-ncoarse+2,m2/) endif ! !MR: missing - nprocy=1 case: two m intervals in one proc! ! allocate(nphis(mexts(1):mexts(2)), nphis1(mexts(1):mexts(2)), & nphis2(mexts(1):mexts(2))) allocate(nexts(mexts(1):mexts(2),2)) allocate(ninds(-nghost:nghost,mexts(1):mexts(2),n1:n2)) ninds=0 do mm=mexts(1),mexts(2) if (lfirst_proc_y) then nphi = ceiling(ncoarse/(mm-m1+1.)) else nphi = ceiling(ncoarse/(m2-mm+1.)) endif if (mod(nzgrid,nphi)/=0) & call warning('initialize_grid','coarsened grid non-periodic: nphi(m='// & trim(itoa(mm))//')='//trim(itoa(nphi))) nphis(mm)=nphi nphis1(mm)=1./nphi nphis2(mm)=1./(nphi*nphi) na=n1 + mod((nphi-mod(nz,nphi))*ipz,nphi) do nn=na,n2,nphi ninds(0,mm,nn)=nn ninds(-1:-nghost:-1,mm,nn)=nn-indgen(nghost)*nphi where(ninds(-nghost:-1,mm,nn)n2) & ninds(+1:+nghost,mm,nn)=(ninds(+1:+nghost,mm,nn)-n2-1)/nphi+1+n2 ne=nn do nni=nn+1,min(n2,nn+nphi-1) idzl=nni-nn+(indgen(nghost)-1)*nphi idzr=nn+nphi-nni+(indgen(nghost)-1)*nphi ninds(:,mm,nni)=(/idzl(nghost:1:-1),-nn,idzr/) enddo enddo do nni=n1,na-1 idzl=nphis(mm)-na+nni+(indgen(nghost)-1)*nphis(mm) idzr=na-nni+(indgen(nghost)-1)*nphis(mm) ninds(:,mm,nni)=(/idzl(nghost:1:-1),0,idzr/) enddo nexts(mm,:)=(/na,ne/) enddo if (lfirst_proc_x.and.lfirst_proc_z) then print*, 'On processors '//trim(itoa(iproc))//' ...(iprocy='//trim(itoa(ipy))// & ')... '//trim(itoa(find_proc(nprocx-1,ipy,nprocz-1)))// & ' the grid is coarsened for m = '// & trim(itoa(mexts(1)))//' ... '//trim(itoa(mexts(2))) print*, 'with coarsening factors:' print'(30(1x,i2))', nphis endif !write(iproc+30,*) 'nexts=', nexts(:,:) !write(iproc+30,*) 'nphis=', nphis(:) !write(iproc+30,'(7(i4,1x))') (ninds(:,mm,:), mm=mexts(1),mexts(2)) else lcoarse=.false. ! now processor-dependent! endif ! ! Print the value for which output is being produced. ! (Have so far only bothered about single processor output.) ! This is now done on all the processors. ! lpoint=min(max(l1,lpoint),l2) mpoint=min(max(m1,mpoint),m2) npoint=min(max(n1,npoint),n2) lpoint2=min(max(l1,lpoint2),l2) mpoint2=min(max(m1,mpoint2),m2) npoint2=min(max(n1,npoint2),n2) ! ! Clean up profile files. ! if (lroot.or..not.lcollective_IO) call remove_prof('z') lwrite_prof=.true. if (lslope_limit_diff) then if (lroot) print*,'initialize_grid: Set up half grid x12, y12, z12' call generate_halfgrid(x12,y12,z12) endif ! endsubroutine initialize_grid !*********************************************************************** subroutine coarsegrid_interp(f,ivar1,ivar2) use General, only: ioptest real, dimension(:,:,:,:) :: f integer, optional :: ivar1,ivar2 integer :: mm,ll,nn,coarse_neigh,iv,iv1,iv2 integer, dimension(3), parameter :: left_inds=(/-3,-2,-1/), right_inds=(/1,2,3/) integer, dimension(6) :: neighs,dists real, dimension(6) :: ws real :: fac iv1=ioptest(ivar1,1) iv2=ioptest(ivar2,size(f,4)) do mm=mexts(1),mexts(2) do nn=n1,n2 if (ninds(0,mm,nn)<=0) then ! ! Interpolation neighbors and distances: ! if (ninds(0,mm,nn)==0) then coarse_neigh=nexts(mm,1) ! right coarsegrid neighbor neighs=(/ninds(left_inds,mm,coarse_neigh),coarse_neigh,ninds(right_inds(1:2),mm,coarse_neigh)/) else coarse_neigh=-ninds(0,mm,nn) ! left coarsegrid neighbor neighs=(/ninds(left_inds(2:3),mm,coarse_neigh),coarse_neigh,ninds(right_inds,mm,coarse_neigh)/) endif dists=(/ninds(left_inds,mm,nn),-ninds(right_inds,mm,nn)/) ! should be precalculated fac=1./nphis(mm)**5 ws(1)=-0.0083333333333333*(dists(2)*dists(3)*dists(4)*dists(5)*dists(6))*fac ws(2)= 0.041666666666667 *(dists(1)*dists(3)*dists(4)*dists(5)*dists(6))*fac ws(3)=-0.083333333333333 *(dists(1)*dists(2)*dists(4)*dists(5)*dists(6))*fac ws(4)= 0.083333333333333 *(dists(1)*dists(2)*dists(3)*dists(5)*dists(6))*fac ws(5)=-0.041666666666667 *(dists(1)*dists(2)*dists(3)*dists(4)*dists(6))*fac ws(6)= 0.0083333333333333*(dists(1)*dists(2)*dists(3)*dists(4)*dists(5))*fac if (abs(sum(ws)-1.)>1e-7) write(iproc+40,'(6(e12.5,1x), e12.5)') ws, sum(ws) do ll=l1,l2 do iv=iv1,iv2 f(ll,mm,nn,iv) = sum(ws*f(ll,mm,neighs,iv)) !if (ll==0.and.iv==5) then ! write(iproc+50,'(i3,7(e14.7,1x))') nn, f(ll,mm,nn,iv), sum(ws*f(ll,mm,neighs,iv)), f(ll,mm,neighs,iv) !endif enddo enddo endif enddo enddo endsubroutine coarsegrid_interp !*********************************************************************** subroutine quintic_interp(nn,a,ninds,dc) integer :: nn real, dimension(:,:) :: a integer, dimension(:) :: ninds real :: dc ! interpolation weights: 1/([-120,24,-12,12,-24,120]*dx^5) = ! [-0.00833333,0.0416667,-0.0833333,0.0833333,-0.0416667,0.00833333]*dx^-5 integer :: iv,izu real, dimension(6), parameter :: coefs = & (/-0.00833333,0.0416667,-0.0833333,0.0833333,-0.0416667,0.00833333/) izu=ipz*nz+nghost ! dels=zgrid(inds)-zgrid(nn) do iv=1,mvar a(nn,iv) = sum(a(ninds,iv)*coefs)/dc**(-5) enddo endsubroutine quintic_interp !*********************************************************************** subroutine save_grid(lrestore) ! ! Saves grid into local statics (needed for downsampled output) ! ! 6-mar-14/MR: coded ! use General, only : loptest ! logical, optional :: lrestore ! real, dimension(mx), save :: xs,dx_1s,dx_tildes real, dimension(my), save :: ys,dy_1s,dy_tildes real, dimension(mz), save :: zs,dz_1s,dz_tildes real, dimension(nx), save :: r_mns,r1_mns,r2_mns real, dimension(my), save :: sinths,sin1ths,sin2ths,cosths, & cotths,cos1ths,tanths real, dimension(mz), save :: sinphs, cosphs real, dimension(nx), save :: rcyl_mns,rcyl_mn1s,rcyl_mn2s real, save :: dxs, dys, dzs logical, save :: lfirst=.true. if (loptest(lrestore)) then if (lfirst) then call fatal_error('save_grid','first call must have lrestore=F') else dx=dxs; dy=dys; dz=dzs x=xs; y=ys; z=zs dx_1=dx_1s; dy_1=dy_1s; dz_1=dz_1s dx_tilde=dx_tildes; dy_tilde=dy_tildes; dz_tilde=dz_tildes if (lspherical_coords) then r_mn=r_mns; r1_mn=r1_mns; r2_mn=r2_mns sinth=sinths; sin1th=sin1ths; sin2th=sin2ths; costh=cosths cotth=cotths; cos1th=cos1ths; tanth=tanths sinphs=sinph; cosphs=cosph elseif (lcylindrical_coords) then rcyl_mn=rcyl_mns; rcyl_mn1=rcyl_mn1s; rcyl_mn2=rcyl_mn2s endif endif elseif (lfirst) then lfirst=.false. dxs=dx; dys=dy; dzs=dz xs=x; ys=y; zs=z dx_1s=dx_1; dy_1s=dy_1; dz_1s=dz_1 dx_tildes=dx_tilde; dy_tildes=dy_tilde; dz_tildes=dz_tilde if (lspherical_coords) then r_mns=r_mn; r1_mns=r1_mn; r2_mns=r2_mn sinths=sinth; sin1ths=sin1th; sin2ths=sin2th; cosths=costh cotths=cotth; cos1ths=cos1th; tanths=tanth sinph=sinphs; cosph=cosphs elseif (lcylindrical_coords) then rcyl_mns=rcyl_mn; rcyl_mn1s=rcyl_mn1; rcyl_mn2s=rcyl_mn2 endif endif endsubroutine save_grid !*********************************************************************** subroutine coords_aux(x,y,z) ! ! 6-mar-14/MR: outsourced from initialize_grid ! real, dimension(:) :: x,y,z ! real, dimension(size(y)) :: lat real, parameter :: sinth_min=1e-5,costh_min=1e-5 ! if (lspherical_coords) then ! ! For spherical coordinates ! r_mn=x(l1:l2) if (x(l1)==0.) then r1_mn(2:)=1./x(l1+1:l2) r1_mn(1)=0. else r1_mn=1./x(l1:l2) endif r2_mn=r1_mn**2 ! ! Calculate sin(theta). Make sure that sinth=1 if there is no y extent, ! regardless of the value of y. This is needed for correct integrations. ! if (ny==1) then sinth=1. else sinth=sin(y) endif ! ! Calculate cos(theta) via latitude, which allows us to ensure ! that sin(lat(midpoint)) = 0 exactly. ! if (luse_latitude) then lat=pi/2-y costh=sin(lat) else costh=cos(y) endif ! ! Calculate 1/sin(theta). To avoid the axis we check that sinth ! is always larger than a minmal value, sinth_min. The problem occurs ! on theta=pi, because the theta range is normally only specified ! with no more than 6 digits, e.g. theta = 0., 3.14159. ! where(abs(sinth)>sinth_min) sin1th=1./sinth elsewhere sin1th=0. endwhere sin2th=sin1th**2 ! ! Calculate cot(theta). ! cotth=costh*sin1th ! ! Calculate 1/cos(theta). To avoid the equator we check that costh ! is always larger than a minmal value, costh_min. The problem occurs ! on theta=pi/2., because the theta range is normally only specified ! with no more than 6 digits, e.g. theta = 0., 3.14159. ! where(abs(costh)>costh_min) cos1th=1./costh elsewhere cos1th=0. endwhere ! ! Calculate tan(theta). ! tanth=sinth*cos1th ! ! also calculate sin(phi) and cos(phi) useful to calculate spherical harmonic ! decomposition. ! if (nzgrid.gt.1) then sinph=sin(z) cosph=cos(z) else sinph=0.;cosph=1 endif ! elseif (lcylindrical_coords) then ! ! Note: for consistency with spherical, 1/rcyl should really be rcyl1_mn, ! not rcyl_mn1. ! rcyl_mn=x(l1:l2) if (x(l1)==0.) then rcyl_mn1(2:)=1./x(l1+1:l2) rcyl_mn1(1)=0. else rcyl_mn1=1./x(l1:l2) endif rcyl_mn2=rcyl_mn1**2 ! endif ! endsubroutine coords_aux !*********************************************************************** subroutine box_vol ! ! calculates box volume ! ! 3-mar-14/MR: outsourced from initialize_grid ! 6-mar-14/MR: changed into subroutine setting global variable box_volume ! box_volume=1. if (lcartesian_coords) then ! ! x,y,z-extent ! if (nxgrid/=1) box_volume = box_volume*Lxyz(1) if (nygrid/=1) box_volume = box_volume*Lxyz(2) if (nzgrid/=1) box_volume = box_volume*Lxyz(3) ! ! Spherical coordinate system ! elseif (lspherical_coords) then ! ! Assume that sinth=1 if there is no theta extent. ! This should always give a volume of 4pi/3*(r2^3-r1^3) for constant integrand ! ! r extent ! if (nxgrid/=1) & box_volume = box_volume*1./3.*(xyz1(1)**3-xyz0(1)**3) ! ! theta extent (if non-radially symmetric) ! if (nygrid/=1) then box_volume = box_volume*(-(cos(xyz1(2)) -cos(xyz0(2)))) else box_volume = box_volume*2. endif ! ! phi extent (if non-axisymmetry) ! if (nzgrid/=1) then box_volume = box_volume*Lxyz(3) else box_volume = box_volume*2.*pi endif ! elseif (lcylindrical_coords) then ! ! Box volume and volume element. ! if (nxgrid/=1) & box_volume = box_volume*.5*(xyz1(1)**2-xyz0(1)**2) ! ! theta extent (non-cylindrically symmetric) ! if (nygrid/=1) then box_volume = box_volume*Lxyz(2) else box_volume = box_volume*2.*pi endif ! ! z extent (vertically extended) ! if (nzgrid/=1) box_volume = box_volume*Lxyz(3) ! endif ! ! Volume calculation for pipe coordinates missing! ! endsubroutine box_vol !*********************************************************************** subroutine pencil_criteria_grid ! ! All pencils that this special module depends on are specified here. ! ! 15-nov-06/tony: coded ! if (any(lfreeze_varext).or.any(lfreeze_varint)) then if (lcylinder_in_a_box.or.lcylindrical_coords) then lpenc_requested(i_rcyl_mn)=.true. else lpenc_requested(i_r_mn)=.true. endif endif ! endsubroutine pencil_criteria_grid !*********************************************************************** subroutine pencil_interdep_grid(lpencil_in) ! ! Interdependency among pencils provided by this module are specified here. ! ! 15-nov-06/tony: coded ! logical, dimension(npencils) :: lpencil_in ! if (lpencil_in(i_rcyl_mn1)) lpencil_in(i_rcyl_mn)=.true. if (lpencil_in(i_r_mn1)) lpencil_in(i_r_mn)=.true. if (lpencil_in(i_evr)) lpencil_in(i_r_mn)=.true. if (lpencil_in(i_evth).or.lpencil_in(i_evr)) then lpencil_in(i_pomx)=.true. lpencil_in(i_pomy)=.true. lpencil_in(i_rcyl_mn)=.true. lpencil_in(i_r_mn1)=.true. endif if ( lpencil_in(i_pomx) & .or. lpencil_in(i_pomy) & .or. lpencil_in(i_phix) & .or. lpencil_in(i_phiy)) then if (lcartesian_coords) then lpencil_in(i_rcyl_mn1)=.true. endif endif if (lspherical_coords.and.lpencil_in(i_phi_mn)) then lpencil_in(i_x_mn)=.true. lpencil_in(i_y_mn)=.true. endif ! if (lpencil_in(i_rr)) then lpencil_in(i_x_mn)=.true. lpencil_in(i_y_mn)=.true. lpencil_in(i_z_mn)=.true. endif ! endsubroutine pencil_interdep_grid !*********************************************************************** subroutine calc_pencils_grid_std(f,p) ! ! Envelope adjusting calc_pencils_hydro_pencpar to the standard use with ! lpenc_loc=lpencil ! ! 10-oct-17/MR: coded ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p ! intent(in) :: f intent(inout) :: p ! call calc_pencils_grid_pencpar(f,p,lpencil) ! endsubroutine calc_pencils_grid_std !*********************************************************************** subroutine calc_pencils_grid_pencpar(f,p,lpenc_loc) ! ! Calculate Grid/geometry related pencils. Uses arbitrary pencil mask lpenc_loc. ! Most basic pencils should come first, as others may depend on them. ! ! 15-nov-06/tony: coded ! 27-aug-07/wlad: generalized for cyl. and sph. coordinates ! real, dimension (mx,my,mz,mfarray) :: f type (pencil_case) :: p logical, dimension(npencils) :: lpenc_loc ! intent(in) :: f intent(inout) :: p intent(in) :: lpenc_loc ! if (lcartesian_coords) then !coordinates vectors if (lpenc_loc(i_x_mn)) p%x_mn = x(l1:l2) if (lpenc_loc(i_y_mn)) p%y_mn = spread(y(m),1,nx) if (lpenc_loc(i_z_mn)) p%z_mn = spread(z(n),1,nx) !spherical distance if (lpenc_loc(i_r_mn)) p%r_mn = sqrt(x(l1:l2)**2+y(m)**2+z(n)**2) !cylindrical distance (pomega) if (lpenc_loc(i_rcyl_mn)) p%rcyl_mn = sqrt(x(l1:l2)**2+y(m)**2) !azimuthal angle (phi) if (lpenc_loc(i_phi_mn)) p%phi_mn = atan2(y(m),x(l1:l2)) !inverse cylindrical distance 1/pomega if (lpenc_loc(i_rcyl_mn1)) p%rcyl_mn1=1./max(p%rcyl_mn,tini) !inverse spherical distance 1/r if (lpenc_loc(i_r_mn1)) p%r_mn1 =1./max(p%r_mn,tini) !pomega unit vectors: pomx=cos(phi) and pomy=sin(phi) where phi=azimuthal angle if (lpenc_loc(i_pomx)) p%pomx = x(l1:l2)*p%rcyl_mn1 if (lpenc_loc(i_pomy)) p%pomy = y( m )*p%rcyl_mn1 !phi unit vectors if (lpenc_loc(i_phix)) p%phix =-y( m )*p%rcyl_mn1 if (lpenc_loc(i_phiy)) p%phiy = x(l1:l2)*p%rcyl_mn1 ! elseif (lcylindrical_coords) then if (lpenc_loc(i_x_mn)) p%x_mn = x(l1:l2)*cos(y(m)) if (lpenc_loc(i_y_mn)) p%y_mn = x(l1:l2)*sin(y(m)) if (lpenc_loc(i_z_mn)) p%z_mn = spread(z(n),1,nx) if (lpenc_loc(i_r_mn)) p%r_mn = sqrt(x(l1:l2)**2+z(n)**2) if (lpenc_loc(i_rcyl_mn)) p%rcyl_mn = x(l1:l2) if (lpenc_loc(i_phi_mn)) p%phi_mn = spread(y(m),1,nx) if (lpenc_loc(i_rcyl_mn1)) p%rcyl_mn1=1./max(p%rcyl_mn,tini) if (lpenc_loc(i_r_mn1)) p%r_mn1 =1./max(p%r_mn,tini) if (lpenc_loc(i_pomx)) p%pomx = 1. if (lpenc_loc(i_pomy)) p%pomy = 0. if (lpenc_loc(i_phix)) p%phix = 0. if (lpenc_loc(i_phiy)) p%phiy = 1. elseif (lspherical_coords) then if (lpenc_loc(i_x_mn)) p%x_mn = x(l1:l2)*sin(y(m))*cos(z(n)) if (lpenc_loc(i_y_mn)) p%y_mn = x(l1:l2)*sin(y(m))*sin(z(n)) if (lpenc_loc(i_z_mn)) p%z_mn = x(l1:l2)*cos(y(m)) if (lpenc_loc(i_r_mn)) p%r_mn = x(l1:l2) if (lpenc_loc(i_rcyl_mn)) p%rcyl_mn = x(l1:l2)*sin(y(m)) if (lpenc_loc(i_phi_mn)) p%phi_mn = spread(z(n),1,nx) if (lpenc_loc(i_rcyl_mn1)) p%rcyl_mn1=1./max(p%rcyl_mn,tini) if (lpenc_loc(i_r_mn1)) p%r_mn1 =1./max(p%r_mn,tini) if (lpenc_loc(i_pomx).or.lpenc_loc(i_pomy).or.& lpenc_loc(i_phix).or.lpenc_loc(i_phiy)) & call fatal_error('calc_pencils_grid', & 'pomx, pomy, phix and phix not implemented for '// & 'spherical polars') endif ! ! set position vector ! if (lpenc_loc(i_rr)) then if (lcartesian_coords) then p%rr(:,1)=p%x_mn p%rr(:,2)=p%y_mn p%rr(:,3)=p%z_mn else call fatal_error('calc_pencils_grid', & 'position vector not implemented for '//& 'non-cartesian coordinates') endif endif ! ! evr is the radial unit vector ! if (lpenc_loc(i_evr)) then if (lcartesian_coords) then p%evr(:,1) = p%rcyl_mn*p%r_mn1*p%pomx p%evr(:,2) = p%rcyl_mn*p%r_mn1*p%pomy p%evr(:,3) = z(n)*p%r_mn1 else call fatal_error('calc_pencils_grid', & 'radial unit vector not implemented for '//& 'non-cartesian coordinates') endif endif ! ! evth is the co-latitudinal unit vector ! if (lpenc_loc(i_evth)) then if (lcartesian_coords) then p%evth(:,1) = z(n)*p%r_mn1*p%pomx p%evth(:,2) = z(n)*p%r_mn1*p%pomy p%evth(:,3) = -p%rcyl_mn*p%r_mn1 else call fatal_error('calc_pencils_grid', & 'co-latitudinal unit vector not implemented for '//& 'non-cartesian coordinates') endif endif ! call keep_compiler_quiet(f) ! endsubroutine calc_pencils_grid_pencpar !*********************************************************************** subroutine grid_profile_0D(xi,grid_func,g,gder1,gder2,param,dxyz,xistep,delta) ! ! Scalar wrapper for the "elemental" subroutine 'grid_profile_1D'. ! ! 14-dec-10/Bourdin.KIS: coded ! real :: xi character(len=*) :: grid_func real :: g real, optional :: gder1, gder2 real, optional :: param real, optional, dimension(3) :: dxyz real, optional, dimension(2) :: xistep, delta ! intent(in) :: xi, grid_func, param, dxyz, xistep, delta intent(out) :: g, gder1, gder2 ! real, dimension(1) :: tmp_g, tmp_gder1, tmp_gder2 ! call grid_profile ((/xi/), grid_func, tmp_g, tmp_gder1, tmp_gder2, param, dxyz, xistep, delta) ! g = tmp_g(1) if (present (gder1)) gder1 = tmp_gder1(1) if (present (gder2)) gder2 = tmp_gder2(1) ! endsubroutine grid_profile_0D !*********************************************************************** subroutine grid_profile_1D(xi,grid_func,g,gder1,gder2,param,dxyz,xistep,delta) ! ! Specify the functional form of the grid profile function g ! and calculate g,g',g''. ! ! 25-jun-04/tobi+wolf: coded ! 9-mar-17/MR: blocked used of unpresent optionals ! real, dimension(:) :: xi character(len=*) :: grid_func real, dimension(size(xi,1)) :: g real, dimension(size(xi,1)), optional :: gder1,gder2 real, optional :: param real, optional, dimension(3) :: dxyz real, optional, dimension(2) :: xistep,delta real :: m ! intent(in) :: xi,grid_func,param,dxyz,xistep,delta intent(out) :: g,gder1,gder2 ! select case (grid_func) ! case ('linear') g=xi if (present(gder1)) gder1=1.0 if (present(gder2)) gder2=0.0 ! case ('sinh', 'sinh2') g=sinh(xi) if (present(gder1)) gder1=cosh(xi) if (present(gder2)) gder2=sinh(xi) ! case ('cos') ! Cos grid: ! Approximately equidistant at the boundaries, linear in the middle if (.not. present (param)) & call fatal_error ('grid_profile', "'cos' needs its parameter.") ! ! Transition goes from slope 1 (xi < pi/2) to slope param (xi > pi/2). m = 0.5 * (param - 1) where (xi <= -pi/2.0) g = xi + pi/2.0 + 1 elsewhere (xi >= pi/2.0) g = param * (xi - pi/2.0) + pi/2.0 + m * pi/2.0 + 1 + pi/2.0 * (1 + m) elsewhere g = xi + m * (xi - cos (xi)) + 1 + pi/2.0 * (1 + m) endwhere if (present (gder1)) then where (xi <= -pi/2.0) gder1 = 1.0 elsewhere (xi >= pi/2.0) gder1 = param elsewhere gder1 = 1 + 0.5 * (m - 1) * (1 + sin (xi)) endwhere endif if (present (gder2)) then where ((xi <= -pi/2.0) .or. (xi >= pi/2.0)) gder2 = 0.0 elsewhere gder2 = 0.5 * (m - 1) * cos (xi) endwhere endif ! case ('arsinh') ! Area sinh grid: ! Approximately equidistant at the boundaries, linear in the middle if (.not. present (param)) & call fatal_error ('grid_profile', "'arsinh' needs its parameter.") ! ! ('asinh' is not available in F95, therefore it is replaced by 'ln'.) g = xi + param * (xi * log (xi + sqrt (xi**2 + 1)) - sqrt (xi**2 + 1)) if (present (gder1)) gder1 = 1.0 + param * log (xi + sqrt (xi**2 + 1)) if (present (gder2)) gder2 = param / (sqrt (xi**2 + 1)) ! case ('tanh') ! Tanh grid: ! Approximately equidistant at the boundaries, linear in the middle if (.not. present (param)) & call fatal_error ('grid_profile', "'tanh' needs its parameter.") ! ! Transition goes from slope 1 (xi -> -oo) to slope param (xi -> +oo). m = 0.5 * (param - 1) g = xi * (m + 1) + m * log (cosh (xi)) if (present (gder1)) gder1 = m * (1 + tanh (xi)) + 1 if (present (gder2)) gder2 = m * (1 - tanh (xi)**2) ! case ('duct') g=sin(xi) if (present(gder1)) gder1= cos(xi) if (present(gder2)) gder2=-sin(xi) ! case ('half-duct') ! duct, but only on one boundary: ! Points are much denser near the boundaries than in the middle g=cos(xi) if (present(gder1)) gder1=-sin(xi) if (present(gder2)) gder2=-cos(xi) ! case ('squared') ! Grid distance increases linearily g=0.5*xi**2 if (present(gder1)) gder1= xi if (present(gder2)) gder2= 0. ! case ('log') ! Grid distance increases logarithmically g=exp(xi) if (present(gder1)) gder1= g if (present(gder2)) gder2= g ! case ('linear+log') ! Grid distance is equidistant near lower boundary and then ! increases logarithmically if (present(param)) then g=(1+xi)*0.5*(1.-tanh(param*xi))+exp(xi)*0.5*(1.+tanh(param*xi)) if (present(gder1)) then gder1= 0.5*((1.-tanh(param*xi))+ & (1+xi)*param*(tanh(param*xi)**2-1.0)+ & exp(xi)*(1.+param+tanh(param*xi)-param*tanh(param*xi)**2)) endif if (present(gder2)) then gder2 = param*(tanh(param*xi)**2-1.) + & param**2*(1+xi)*(1.-tanh(param*xi)**2)*tanh(param*xi) + & 0.5*exp(xi)*(1.+2*param+(1.-2*param**2)*tanh(param*xi)- & 2*param*tanh(param*xi)**2+2*param**2*tanh(param*xi)**3) endif else g=exp(xi) if (present(gder1)) gder1= g if (present(gder2)) gder2= g endif ! case ('power-law') ! Grid distance increases according to a power-law if (.not. present (param)) & call fatal_error ('grid_profile', "'power-law' needs its parameter.") g=xi**(1./param) if (present(gder1)) gder1=1./param* xi**(1/param-1) if (present(gder2)) gder2=1./param*(1./param-1.)*xi**(1/param-2) ! case ('trough') ! Grid distance larger and linear near boundaries and smaller and linear in ! middle. Needs xistep(2), delta(2) and param parameters specified in ! start.in. E.g.,for solar atmosphere can use delta=0.4,0.5, xistar=100,300, ! param=0.02 on a grid of mzgrid=646 if (.not. present (param)) & call fatal_error ('grid_profile', "'trough' needs its parameter.") g=(-alog(exp(param*(xi-xistep(1)))+exp(-param*(xi-xistep(1))))/param+ & alog(exp(param*(xi-xistep(2)))+exp(-param*(xi-xistep(2))))/param+ & (2.+delta(1))*xi)/(2.+delta(1)) if (present(gder1)) then gder1=(2.+delta(1)-tanh(param*(xi-xistep(1)))+ & tanh(param*(xi-xistep(2))))/(2.+delta(1)) endif if (present(gder2)) then gder2=(-param*(1.-(tanh(param*(xi-xistep(1))))**2)+ & param*(1.-(tanh(param*(xi-xistep(2))))**2))/ & (2.+delta(1)) endif ! case ('frozensphere') ! Just like sinh, except set dx constant below a certain radius. m = 4. where (xi<0) g = m*xi elsewhere g=sinh(xi) endwhere if (present(gder1)) then where (xi<0) gder1 = m elsewhere gder1=cosh(xi) endwhere endif if (present(gder2)) then where (xi<0) gder2 = 0. elsewhere gder2=sinh(xi) endwhere endif ! case ('step-linear') if (.not. (present(dxyz) .and. present(xistep) .and. present(delta))) & call fatal_error('grid_profile',"'step-linear' needs its parameters.") if (xistep(1)/=0.0) then g= & dxyz(1)*0.5*(xi-delta(1)*log(cosh(dble((xi-xistep(1))/delta(1))))) + & dxyz(2)*0.5*( delta(1)*log(cosh(dble((xi-xistep(1))/delta(1)))) - & delta(2)*log(cosh(dble((xi-xistep(2))/delta(2)))) )+ & dxyz(3)*0.5*(xi+delta(2)*log(cosh(dble((xi-xistep(2))/delta(2)))) ) ! if (present(gder1)) then gder1= & dxyz(1)*0.5*(1.0-tanh(dble((xi-xistep(1))/delta(1))) ) + & dxyz(2)*0.5*( tanh(dble((xi-xistep(1))/delta(1))) - & tanh(dble((xi-xistep(2))/delta(2))) ) + & dxyz(3)*0.5*(1.0+tanh(dble((xi-xistep(2))/delta(2))) ) ! endif if (present(gder2)) then gder2= & + 0.5/delta(1)*(dxyz(2)-dxyz(1))/cosh(dble((xi-xistep(1))/delta(1)))**2 & + 0.5/delta(2)*(dxyz(3)-dxyz(2))/cosh(dble((xi-xistep(2))/delta(2)))**2 endif else ! if xistep(1)=0.0, only a single step is needed g= & dxyz(2)*0.5*(xi-delta(2)*log(cosh(dble((xi-xistep(2))/delta(2)))) )+ & dxyz(3)*0.5*(xi+delta(2)*log(cosh(dble((xi-xistep(2))/delta(2)))) ) ! if (present(gder1)) then gder1= & dxyz(2)*0.5*(1.0-tanh(dble((xi-xistep(2))/delta(2))) ) + & dxyz(3)*0.5*(1.0+tanh(dble((xi-xistep(2))/delta(2))) ) ! endif if (present(gder2)) then gder2= & + 0.5/delta(2)*(dxyz(3)-dxyz(2))/cosh(dble((xi-xistep(2))/delta(2)))**2 endif endif ! case default call fatal_error('grid_profile','grid function not implemented: '//trim(grid_func)) ! endselect ! endsubroutine grid_profile_1D !*********************************************************************** function find_star(xi_lo,xi_up,x_lo,x_up,x_star,grid_func) result (xi_star) ! ! Finds the xi that corresponds to the inflection point of the grid-function ! by means of a newton-raphson root-finding algorithm. ! ! 25-jun-04/tobi+wolf: coded ! real, intent(in) :: xi_lo,xi_up,x_lo,x_up,x_star character(len=*), intent(in) :: grid_func ! real :: xi_star,dxi,tol real :: g_lo,gder_lo real :: g_up,gder_up real :: f ,fder integer, parameter :: maxit=1000 logical :: lreturn integer :: it ! if (xi_lo>=xi_up) & call fatal_error('find_star','xi1 >= xi2 -- this should not happen') ! tol=epsi*(xi_up-xi_lo) xi_star= (xi_up+xi_lo)/2 ! lreturn=.false. ! do it=1,maxit ! call grid_profile(xi_lo-xi_star,grid_func,g_lo,gder_lo) call grid_profile(xi_up-xi_star,grid_func,g_up,gder_up) ! f =-(x_up-x_star)*g_lo +(x_lo-x_star)*g_up fder= (x_up-x_star)*gder_lo-(x_lo-x_star)*gder_up ! dxi=f/fder xi_star=xi_star-dxi ! if (lreturn) return ! if (abs(dxi) 0) then dir: do i = 1, 3 if (lactive_dimension(i)) then call inverse_grid(i, x(:,i), xi(:,i), local=.true.) else xi(:,i) = ngp1 endif enddo dir endif nonzero ! endsubroutine real_to_index !*********************************************************************** subroutine inverse_grid(dir, x, xi, local) ! ! Transform the x coordinates in real space to the xi coordinates in ! index space in dir direction, where dir = 1, 2, or, 3. ! If local is present and .true., the index space is with respect to ! the local grid. ! ! 24-dec-14/ccyang: coded. ! use General, only: arcsinh ! integer, intent(in) :: dir real, dimension(:), intent(in) :: x real, dimension(:), intent(out) :: xi logical, intent(in), optional :: local ! integer, dimension(size(xi)) :: inear character(len=linelen) :: msg logical :: loc integer :: shift, inear_max real :: h, a, b, c ! ! Sanity check. ! if (any(lshift_origin) .or. any(lshift_origin_lower)) & call fatal_error('inverse_grid', 'lshift_origin and lshift_origin_lower are not supported. ') ! ! Global or local index space? ! loc = .false. if (present(local)) loc = local ! ! Check the direction. ! h = 0.0 shift = 0 ckdir: select case (dir) case (1) ckdir h = dx if (loc) shift = nx * ipx if (loc) inear_max = nghost + nx case (2) ckdir h = dy if (loc) shift = ny * ipy if (loc) inear_max = nghost + ny case (3) ckdir h = dz if (loc) shift = nz * ipz if (loc) inear_max = nghost + nz case default ckdir write(msg,*) 'unknown direction dir = ', dir call fatal_error('inverse_grid', trim(msg)) endselect ckdir ! ! Make the inversion according to the grid function. ! func: select case (grid_func(dir)) ! case ('linear') func xi = (x - xyz0(dir)) / h ! case ('sinh') func a = coeff_grid(dir) * Lxyz(dir) b = sinh(a) c = cosh(a) - 1.0 a = (xyz_star(dir) - xyz0(dir)) / Lxyz(dir) a = a * b / sqrt(1.0 + 2.0 * a * (1.0 - a) * c) b = (sqrt(1.0 + a * a) * b - a * c) / Lxyz(dir) xi = (arcsinh(a) + arcsinh(b * (x - xyz0(dir)) - a)) / (coeff_grid(dir) * h) ! case default func call fatal_error('inverse_grid in grid.f90', & 'unknown grid function ' // trim(grid_func(dir))) ! endselect func ! ! Shift to match the global index space. ! if (lperi(dir)) then xi = xi + real(nghost) + 0.5 else xi = xi + real(nghost + 1) endif ! ! Convert to the local index space if requested. ! getloc: if (loc) then xi = xi - real(shift) inear = nint(xi) where ((x < xyz1(dir) .and. inear > inear_max) .or. & (x < xyz0(dir) .and. inear > nghost)) xi = nearest(xi, -1.0) endif getloc ! endsubroutine inverse_grid !*********************************************************************** subroutine symmetrize_grid(grid,ngrid,idim) integer :: ngrid,idim real, dimension(ngrid) :: grid if ( lequidist(idim) .and. xyz0(idim)+xyz1(idim)