submodule (inplane_operators_m) inplane_operators_s implicit none contains module subroutine init_inplane_operators(self, equi, mesh, equi_on_mesh, norm_scale, norm_flutter) class(inplane_operators_t), intent(inout) :: self class(equilibrium_t), target, intent(inout) :: equi type(mesh_cart_t), target, intent(inout) :: mesh class(equilibrium_storage_t), target, intent(inout) :: equi_on_mesh real(GP), intent(in) :: norm_scale real(GP), intent(in) :: norm_flutter real(GP) :: hf hf = mesh%get_spacing_f()/norm_scale self%hf_inv = 1.0_GP / hf self%hf_inv_sqrd = 1.0_GP / hf**2 self%delta = 1.0_GP / norm_scale self%deltabeta = norm_flutter self%equi => equi self%mesh => mesh self%equi_on_mesh => equi_on_mesh self%n_points = mesh%get_n_points() self%n_points_inner = mesh%get_n_points_inner() allocate(self%wrk1(self%n_points)) allocate(self%wrk2(self%n_points)) ! Select curvature operators for different geometries select type(equi) type is(dommaschk_t) ! btor based curvature operators self%curvature => curvature_btor_based self%curvature_sqrd => curvature_sqrd_btor_based class default ! ddy based curvature operators, ! can only bu used under the assumption Btor \propto 1/x, self%curvature => curvature_ddy_based self%curvature_sqrd => curvature_sqrd_ddy_based end select end subroutine pure real(GP) module function ddx(self, u, l) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: u integer, intent(in) :: l integer :: lwest, least lwest = self%mesh%get_index_neighbor(-1, 0, l) least = self%mesh%get_index_neighbor( 1, 0, l) ddx = (u(least) - u(lwest)) * 0.5_GP * self%hf_inv end function pure real(GP) module function ddy(self, u, l) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: u integer, intent(in) :: l integer :: lsouth, lnorth lsouth = self%mesh%get_index_neighbor( 0, -1, l) lnorth = self%mesh%get_index_neighbor( 0, 1, l) ddy = (u(lnorth) - u(lsouth)) * 0.5_GP * self%hf_inv end function real(GP) module function nabla_pol_sqrd(self, u, l, co) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: u integer, intent(in) :: l real(GP), dimension(self%n_points), optional, intent(in) :: co integer :: l_lft, l_rgt, l_btm, l_top real(GP) :: x, y, phi, x_lft, x_rgt real(GP) :: co_mx_l, co_mx_r, co_my_b, co_my_t, jac, jac_lft, jac_rgt l_lft = self%mesh%get_index_neighbor(-1, 0, l) l_rgt = self%mesh%get_index_neighbor(1, 0, l) l_btm = self%mesh%get_index_neighbor(0, -1, l) l_top = self%mesh%get_index_neighbor(0, 1, l) ! Jacobian (in x) phi = self%mesh%get_phi() x = self%mesh%get_x(l) y = self%mesh%get_y(l) x_lft = x - self%mesh%get_spacing_f() x_rgt = x + self%mesh%get_spacing_f() jac = self%equi%jacobian(x, y, phi) jac_lft = self%equi%jacobian(x_lft, y, phi) jac_rgt = self%equi%jacobian(x_rgt, y, phi) ! Map J*co to x and y staggered mesh if (present(co)) then co_mx_l = 0.5_GP * (co(l_lft) * jac_lft + co(l) * jac) co_mx_r = 0.5_GP * (co(l) * jac + co(l_rgt) * jac_rgt) co_my_b = 0.5_GP * jac * (co(l_btm) + co(l)) co_my_t = 0.5_GP * jac * (co(l) + co(l_top)) else co_mx_l = 0.5_GP * (jac_lft + jac) co_mx_r = 0.5_GP * (jac + jac_rgt) co_my_b = jac co_my_t = jac endif nabla_pol_sqrd = ( co_mx_r * (u(l_rgt) - u(l)) & - co_mx_l * (u(l) - u(l_lft)) & + co_my_t * (u(l_top) - u(l)) & - co_my_b * (u(l) - u(l_btm))) & * self%hf_inv_sqrd / jac end function pure real(GP) module function arakawa(self, u, v, l) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: u real(GP), dimension(self%n_points), intent(in) :: v integer, intent(in) :: l integer :: in, jn, ln real(GP), dimension(-1:1,-1:1) :: zeta, psi do jn = -1, 1 do in = -1, 1 ln = self%mesh%get_index_neighbor(in, jn, l) zeta(in,jn) = u(ln) psi(in,jn) = v(ln) enddo enddo arakawa = self%arakawa_loc(zeta, psi) end function pure real(GP) module function doubled_arakawa(self, u, v, l) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: u real(GP), dimension(self%n_points), intent(in) :: v integer, intent(in) :: l integer :: in, jn, ln real(GP), dimension(-2:2, -2:2) :: u_loc, v_loc ! Values on local Cartesian grid around grid point l (stencil) ! Map values to local 5x5 grid do jn = -2, 2 do in = -2, 2 ln = self%mesh%get_index_neighbor(in, jn, l) u_loc(in,jn) = u(ln) v_loc(in,jn) = v(ln) enddo enddo doubled_arakawa = self%doubled_arakawa_loc(u_loc, v_loc) end function real(GP) module function curvature_btor_based(self, u, l) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: u integer, intent(in) :: l integer :: in, jn, ln real(GP) :: xn, yn, phi, x, y real(GP), dimension(-1:1,-1:1) :: zeta, psi phi = self%mesh%get_phi() do jn = -1, 1 do in = -1, 1 ln = self%mesh%get_index_neighbor(in, jn, l) xn = self%mesh%get_x(ln) yn = self%mesh%get_y(ln) zeta(in,jn) = self%equi%jacobian(xn, yn, phi) / self%equi_on_mesh%btor(ln) psi(in,jn) = u(ln) enddo enddo x = self%mesh%get_x(l) y = self%mesh%get_y(l) curvature_btor_based = -self%delta * self%arakawa_loc(zeta, psi) & / self%equi%jacobian(x, y, phi) end function pure real(GP) module function curvature_ddy_based(self, u, l) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: u integer, intent(in) :: l integer :: lnorth, lsouth lnorth = self%mesh%get_index_neighbor(0, 1, l) lsouth = self%mesh%get_index_neighbor(0, -1, l) curvature_ddy_based = - (u(lnorth) - u(lsouth)) * self%hf_inv end function real(GP) module function curvature_sqrd_btor_based(self, u, l) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: u integer, intent(in) :: l integer :: in, jn, ln real(GP) :: x, y, phi real(GP), dimension(-2:2, -2:2) :: jac_loc, a_loc, u_loc phi = self%mesh%get_phi() ! Map values to local 5x5 grid do jn = -2, 2 do in = -2, 2 ln = self%mesh%get_index_neighbor(in, jn, l) x = self%mesh%get_x(ln) y = self%mesh%get_y(ln) jac_loc(in, jn) = self%equi%jacobian(x, y, phi) a_loc(in, jn) = -1.0_GP / self%equi_on_mesh%btor(ln) u_loc(in, jn) = u(ln) enddo enddo curvature_sqrd_btor_based = self%delta**2 * self%doubled_arakawa_loc(a_loc, u_loc, jac_loc) end function pure real(GP) module function curvature_sqrd_ddy_based(self, u, l) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: u integer, intent(in) :: l integer :: lnorth, lsouth lnorth = self%mesh%get_index_neighbor(0, 1, l) lsouth = self%mesh%get_index_neighbor(0, -1, l) curvature_sqrd_ddy_based = (u(lnorth) + u(lsouth) - 2*u(l)) * 4.0_GP * self%hf_inv_sqrd end function real(GP) module function polarisation_advection(self, dia_fac, co, adv_pot, potv, tiv, nev, l) class(inplane_operators_t), intent(in) :: self real(GP), intent(in) :: dia_fac real(GP), dimension(self%n_points), intent(in) :: co real(GP), dimension(self%n_points), intent(in) :: adv_pot real(GP), dimension(self%n_points), intent(in) :: potv real(GP), dimension(self%n_points), intent(in) :: tiv real(GP), dimension(self%n_points), intent(in) :: nev integer, intent(in) :: l integer :: in, jn, ln real(GP), dimension(-2:2, -2:2) :: co_loc, adv_pot_loc, pot_loc, pi_loc, ne_loc, btor_loc !! values on local Cartesian grid around grid point l (stencil) real(GP), dimension(-2:1, -1:1) :: dvdx !! x-derivative (electric field) on x-staggered grid real(GP), dimension(-2:1, -1:1) :: co_mx, adv_pot_mx !! Mapped to x-staggered grid real(GP), dimension(-1:1, -2:1) :: dvdy !! y-derivative (electric field) on y-staggered grid real(GP), dimension(-1:1, -2:1) :: co_my, adv_pot_my !! Mapped to y-staggered grid real(GP) :: ne_mx, dpidx, dpotdx real(GP) :: ne_my, dpidy, dpotdy real(GP) :: arxl, arxr, aryt, aryb real(GP) :: btorl, btorr, btort, btorb ! Map values to local 5x5 grid do jn = -2, 2 do in = -2, 2 ln = self%mesh%get_index_neighbor(in, jn, l) co_loc(in,jn) = co(ln) adv_pot_loc(in,jn) = adv_pot(ln) pot_loc(in,jn) = potv(ln) pi_loc(in,jn) = tiv(ln)*nev(ln) ne_loc(in,jn) = nev(ln) btor_loc(in,jn) = self%equi_on_mesh%btor(ln) enddo enddo ! Map and compute necessary x-derivatives on x-staggered grid do jn = -1,1 do in = -2,1 co_mx(in, jn) = 0.5_GP * ( co_loc(in+1, jn) + co_loc(in, jn) ) adv_pot_mx(in, jn) = 0.5_GP * ( adv_pot_loc(in+1, jn) + adv_pot_loc(in, jn) ) ne_mx = 0.5_GP * ( ne_loc(in+1, jn) + ne_loc(in, jn) ) dpidx = ( pi_loc(in+1, jn) - pi_loc(in, jn) ) * self%hf_inv dpotdx = ( pot_loc(in+1, jn) - pot_loc(in, jn) ) * self%hf_inv dvdx(in, jn) = dpotdx + dia_fac * dpidx / ne_mx enddo enddo ! Map and compute necessary y-derivatives on y-staggered grid do jn = -2,1 do in = -1,1 co_my(in, jn) = 0.5_GP * ( co_loc(in, jn+1) + co_loc(in, jn) ) adv_pot_my(in, jn) = 0.5_GP * ( adv_pot_loc(in, jn+1) + adv_pot_loc(in, jn) ) ne_my = 0.5_GP * ( ne_loc(in, jn+1) + ne_loc(in, jn) ) dpidy = ( pi_loc(in, jn+1) - pi_loc(in, jn) ) * self%hf_inv dpotdy = ( pot_loc(in, jn+1) - pot_loc(in, jn) ) * self%hf_inv dvdy(in, jn) = dpotdy + dia_fac * dpidy / ne_my enddo enddo ! Compute toroidal field at staggered points btorl = 0.5_GP * (btor_loc(-1,0) + btor_loc(0,0)) btorr = 0.5_GP * (btor_loc(1,0) + btor_loc(0,0)) btorb = 0.5_GP * (btor_loc(0,-1) + btor_loc(0,0)) btort = 0.5_GP * (btor_loc(0,1) + btor_loc(0,0)) ! Compute Arakawa brackets (x-staggered ) for points (in=-1, jn=0) and (in=0, jn=0) arxl = self%arakawa_loc(adv_pot_mx(-2:0,-1:1), dvdx(-2:0,-1:1)) / btorl arxr = self%arakawa_loc(adv_pot_mx(-1:1,-1:1), dvdx(-1:1,-1:1)) / btorr ! Compute Arakawa brackets (y-staggered ) for points (in=0, jn=-1) and (in=0, jn=0) aryb = self%arakawa_loc(adv_pot_my(-1:1,-2:0), dvdy(-1:1,-2:0)) / btorb aryt = self%arakawa_loc(adv_pot_my(-1:1,-1:1), dvdy(-1:1,-1:1)) / btort polarisation_advection = ( co_mx(0,0) * arxr - co_mx(-1,0) * arxl & + co_my(0,0) * aryt - co_my(0,-1) * aryb ) & * self%hf_inv end function module subroutine apply_dissipation_inplane(self, u, ndiss, du, co) class(inplane_operators_t), intent(inout) :: self real(GP), dimension(self%n_points), intent(in) :: u integer, intent(in) :: ndiss real(GP), dimension(self%n_points_inner), intent(out) :: du real(GP), dimension(self%n_points), optional, intent(in) :: co integer:: ki, l, idiss , sgn if (mod(ndiss,2) == 1) then sgn = 1 else sgn = -1 endif !$omp do private(l) do l = 1, self%n_points self%wrk1(l) = u(l) self%wrk2(l) = 0.0_GP enddo !$omp end do do idiss = 1, ndiss if (idiss < ndiss) then !$omp do private(ki, l) do ki = 1, self%n_points_inner l = self%mesh%inner_indices(ki) self%wrk2(l) = self%nabla_pol_sqrd(self%wrk1, l) enddo !$omp end do !$omp do private(l) do l = 1, self%n_points self%wrk1(l) = self%wrk2(l) enddo !$omp end do else !$omp do private(ki, l) do ki = 1, self%n_points_inner l = self%mesh%inner_indices(ki) du(ki) = sgn * self%nabla_pol_sqrd(self%wrk1, l, co) enddo !$omp end do endif enddo end subroutine pure real(GP) module function arakawa_loc(self, zeta, psi) class(inplane_operators_t), intent(in) :: self real(GP), dimension(-1:1,-1:1), intent(in) :: zeta real(GP), dimension(-1:1,-1:1), intent(in) :: psi real(GP) :: jplusplus, jpluscross, jcrossplus real(GP), parameter :: one_over_twelve = 1.0_GP / 12.0_GP jplusplus = ( zeta(1,0) - zeta(-1,0) ) * ( psi(0,1) - psi(0,-1) ) & - ( zeta(0,1) - zeta(0,-1) ) * ( psi(1,0) - psi(-1,0) ) jpluscross = zeta(1,0) * ( psi(1,1) - psi(1,-1) ) & - zeta(-1,0) * ( psi(-1,1) - psi(-1,-1) ) & - zeta(0,1) * ( psi(1,1) - psi(-1,1) ) & + zeta(0,-1) * ( psi(1,-1) - psi(-1,-1) ) jcrossplus = zeta(1,1) * ( psi(0,1) - psi(1,0) ) & - zeta(-1,-1) * ( psi(-1,0) - psi(0,-1) ) & - zeta(-1,1) * ( psi(0,1) - psi(-1,0) ) & + zeta(1,-1) * ( psi(1,0) - psi(0,-1) ) arakawa_loc = (jplusplus + jpluscross + jcrossplus ) * self%hf_inv_sqrd * one_over_twelve end function pure real(GP) module function doubled_arakawa_loc(self, u_loc, v_loc, jac_loc) class(inplane_operators_t), intent(in) :: self real(GP), dimension(-2:2,-2:2), intent(in) :: u_loc real(GP), dimension(-2:2,-2:2), intent(in) :: v_loc real(GP), dimension(-2:2,-2:2), intent(in), optional :: jac_loc real(GP), dimension(-2:2,-2:2) :: uj_loc integer :: in, jn real(GP), dimension(-2:1, -1:1) :: u_mx, v_mx ! Mapped to x-staggered grid real(GP), dimension(-1:1, -2:1):: u_my, v_my ! Mapped to y-staggered grid real(GP) :: arx_l, arx_r, ary_t, ary_b, & dudy_mx_l, dudy_mx_r, dudx_my_b, dudx_my_t, & jac_c, jac_mx_l, jac_mx_r, jac_my_b, jac_my_t if (present(jac_loc)) then do jn = -2,2 do in = -2,2 uj_loc(in, jn) = u_loc(in, jn) * jac_loc(in, jn) enddo enddo else uj_loc = u_loc endif ! Map to x-staggered grid do jn = -1,1 do in = -2,1 u_mx(in, jn) = 0.5_GP * (uj_loc(in+1, jn) + uj_loc(in, jn)) v_mx(in, jn) = 0.5_GP * (v_loc(in+1, jn) + v_loc(in, jn)) enddo enddo ! Map to y-staggered grid do jn = -2,1 do in = -1,1 u_my(in, jn) = 0.5_GP * (uj_loc(in, jn+1) + uj_loc(in, jn)) v_my(in, jn) = 0.5_GP * (v_loc(in, jn+1) + v_loc(in, jn)) enddo enddo ! Compute Arakawa brackets (x-staggered ) for points (in=-1, jn=0) and (in=0, jn=0) arx_l = self%arakawa_loc(u_mx(-2:0,-1:1), v_mx(-2:0,-1:1)) arx_r = self%arakawa_loc(u_mx(-1:1,-1:1), v_mx(-1:1,-1:1)) ! Compute Arakawa brackets (y-staggered ) for points (in=0, jn=-1) and (in=0, jn=0) ary_b = self%arakawa_loc(u_my(-1:1,-2:0), v_my(-1:1,-2:0)) ary_t = self%arakawa_loc(u_my(-1:1,-1:1), v_my(-1:1,-1:1)) ! Compute dudy on x-staggered grid dudy_mx_l = 0.25_GP * self%hf_inv * (uj_loc(0,1) - uj_loc(0,-1) + uj_loc(-1,1) - uj_loc(-1,-1)) dudy_mx_r = 0.25_GP * self%hf_inv * (uj_loc(1,1) - uj_loc(1,-1) + uj_loc(0,1) - uj_loc(0,-1)) ! Compute dudx on y-staggered grid dudx_my_b = 0.25_GP * self%hf_inv * (uj_loc(1,0) - uj_loc(-1,0) + uj_loc(1,-1) - uj_loc(-1,-1)) dudx_my_t = 0.25_GP * self%hf_inv * (uj_loc(1,1) - uj_loc(-1,1) + uj_loc(1,0) - uj_loc(-1,0)) ! Map jacobian to x and y staggered grid if (present(jac_loc)) then jac_c = jac_loc(0, 0) jac_mx_l = 0.5_GP * (jac_loc(-1, 0) + jac_loc(0, 0)) jac_mx_r = 0.5_GP * (jac_loc(0, 0) + jac_loc(1, 0)) jac_my_b = 0.5_GP * (jac_loc(0, -1) + jac_loc(0, 0)) jac_my_t = 0.5_GP * (jac_loc(0, 0) + jac_loc(0, 1)) else jac_c = 1.0_GP jac_mx_l = 1.0_GP jac_mx_r = 1.0_GP jac_my_b = 1.0_GP jac_my_t = 1.0_GP endif doubled_arakawa_loc = self%hf_inv / jac_c * ( & - arx_r * dudy_mx_r / jac_mx_r + arx_l * dudy_mx_l / jac_mx_l & + ary_t * dudx_my_t / jac_my_t - ary_b * dudx_my_b / jac_my_b ) end function pure real(GP) module function flutter_grad(self, apar, u, l) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: apar real(GP), dimension(self%n_points), intent(in) :: u integer, intent(in) :: l integer :: in, jn, ln real(GP), dimension(-1:1,-1:1) :: zeta, psi do jn = -1, 1 do in = -1, 1 ln = self%mesh%get_index_neighbor(in, jn, l) zeta(in,jn) = apar(ln) / self%equi_on_mesh%btor(ln) psi(in,jn) = u(ln) enddo enddo flutter_grad = self%deltabeta * self%arakawa_loc(zeta, psi) end function pure real(GP) module function flutter_div(self, apar, u, l) class(inplane_operators_t), intent(in) :: self real(GP), dimension(self%n_points), intent(in) :: apar real(GP), dimension(self%n_points), intent(in) :: u integer, intent(in) :: l integer :: in, jn, ln real(GP), dimension(-1:1,-1:1) :: zeta, psi do jn = -1, 1 do in = -1, 1 ln = self%mesh%get_index_neighbor(in, jn, l) zeta(in,jn) = apar(ln) / self%equi_on_mesh%btor(ln) psi(in,jn) = u(ln) / self%equi_on_mesh%btor(ln) enddo enddo flutter_div = self%deltabeta * self%equi_on_mesh%btor(l) * self%arakawa_loc(zeta, psi) end function module subroutine destructor(self) type(inplane_operators_t), intent(inout) :: self self%equi => null() self%mesh => null() self%equi_on_mesh => null() self%curvature => null() self%curvature_sqrd => null() if (allocated(self%wrk2)) deallocate(self%wrk2) if (allocated(self%wrk1)) deallocate(self%wrk1) end subroutine end submodule