TokaMaker: Add subroutines for computing local shear and NT H-mode access proxy (OpenFUSIONToolkit#338)

- Add function to compute simple NT H-mode access proxy from [Nelson et al. (2022)]
 - Add support for computing local shear (`TokaMaker_equilibrium.calc_local_shear()`)
 - Add support for interpolating general nodal fields in TokaMaker through `TokaMaker.get_field_eval()`
 - Add support for evaluating multiple points at once in interpolation objects returned by `TokaMaker.get_field_eval()`
    * Add start of generic FE interpolation support (just surface Lagrange for now)
 - Add the ability to retrieve fields on node points using FE projection (`TokaMaker_equilibrium.get_nodal_field()`)
 - Add support for uniform resampling in `TokaMaker.trace_surf()`
 - Fix bug in Ip and centroid calculations in `gs_comp_globals` when anisotropic pressure is used
    * This change impacts `TokaMaker_equilibrium.get_stats()` and all downstream consumers (eg. `print_info()`)

Author: hansec

src/base/oft_local.F90

@@ -15,6 +15,7 @@
 !--------------------------------------------------------------------------------
 MODULE oft_local
 USE, INTRINSIC :: iso_c_binding, only: c_int, c_ptr, c_long
+USE, INTRINSIC :: ieee_arithmetic, only: ieee_value, ieee_quiet_nan
 #ifdef __INTEL_COMPILER
 USE ifport ! Intel fortran portability library
 #endif

src/examples/TokaMaker/DIIID/DIIID_baseline_ex.ipynb

@@ -76,6 +76,9 @@ MODULE oft_blag_operators
 !------------------------------------------------------------------------------
 type, extends(bfem_interp) :: oft_lag_bvrinterp
   class(oft_vector), pointer :: u => NULL() !< Field for interpolation
+  class(oft_vector), pointer :: ux => NULL() !< Field for interpolation
+  class(oft_vector), pointer :: uy => NULL() !< Field for interpolation
+  class(oft_vector), pointer :: uz => NULL() !< Field for interpolation
   real(r8), pointer, dimension(:,:) :: vals => NULL() !< Local values
   class(oft_scalar_bfem), pointer :: lag_rep => NULL() !< Lagrange FE representation
 contains
@@ -87,6 +90,15 @@ MODULE oft_blag_operators
   procedure :: delete => lag_bvrinterp_delete
 end type oft_lag_bvrinterp
 !------------------------------------------------------------------------------
+!> Interpolate \f$ \nabla \f$ of a Lagrange field
+!------------------------------------------------------------------------------
+type, extends(oft_lag_bvrinterp) :: oft_lag_bvcinterp
+  logical :: cylindrical = .FALSE. !< Curl in cylindrical coordinates?
+contains
+  !> Reconstruct field
+  procedure :: interp => lag_bvcinterp
+end type oft_lag_bvcinterp
+!------------------------------------------------------------------------------
 !> Needs docs
 !------------------------------------------------------------------------------
 type, extends(oft_solver_bc) :: oft_blag_zerob
@@ -256,12 +268,23 @@ subroutine lag_bvrinterp_setup(self,lag_rep)
 self%mesh=>self%lag_rep%mesh
 !---Get local slice
 IF(.NOT.ASSOCIATED(self%vals))ALLOCATE(self%vals(3,self%lag_rep%ne))
-vtmp=>self%vals(1,:)
-CALL self%u%get_local(vtmp,1)
-vtmp=>self%vals(2,:)
-CALL self%u%get_local(vtmp,2)
-vtmp=>self%vals(3,:)
-CALL self%u%get_local(vtmp,3)
+IF(ASSOCIATED(self%u))THEN
+  vtmp=>self%vals(1,:)
+  CALL self%u%get_local(vtmp,1)
+  vtmp=>self%vals(2,:)
+  CALL self%u%get_local(vtmp,2)
+  vtmp=>self%vals(3,:)
+  CALL self%u%get_local(vtmp,3)
+ELSE IF(ASSOCIATED(self%ux))THEN
+  vtmp=>self%vals(1,:)
+  CALL self%ux%get_local(vtmp)
+  vtmp=>self%vals(2,:)
+  CALL self%uy%get_local(vtmp)
+  vtmp=>self%vals(3,:)
+  CALL self%uz%get_local(vtmp)
+ELSE
+  CALL oft_abort("No vector associated","lag_bvrinterp_setup",__FILE__)
+END IF
 end subroutine lag_bvrinterp_setup
 !------------------------------------------------------------------------------
 !> Destroy temporary internal storage
@@ -301,6 +324,49 @@ subroutine lag_bvrinterp(self,cell,f,gop,val)
 DEBUG_STACK_POP
 end subroutine lag_bvrinterp
 !------------------------------------------------------------------------------
+!> Reconstruct the curl of a boundary Lagrange vector field
+!------------------------------------------------------------------------------
+subroutine lag_bvcinterp(self,cell,f,gop,val)
+class(oft_lag_bvcinterp), intent(inout) :: self
+integer(i4), intent(in) :: cell !< Cell for interpolation
+real(r8), intent(in) :: f(:) !< Position in cell in logical coord [3]
+real(r8), intent(in) :: gop(3,3) !< Logical gradient vectors at f [3,3]
+real(r8), intent(out) :: val(:) !< Reconstructed field at f [3]
+integer(i4), allocatable :: j(:)
+integer(i4) :: jc
+real(r8) :: rop(3),fop,pt(3)
+DEBUG_STACK_PUSH
+!---
+IF(.NOT.ASSOCIATED(self%vals))CALL oft_abort('Setup has not been called!', &
+  'lag_bvcinterp',__FILE__)
+!---Get dofs
+allocate(j(self%lag_rep%nce))
+call self%lag_rep%ncdofs(cell,j) ! get DOFs
+!---Reconstruct field
+IF(self%cylindrical)THEN
+  val=0.d0
+  pt=self%mesh%log2phys(cell,f)
+  do jc=1,self%lag_rep%nce
+    call oft_blag_eval(self%lag_rep,cell,jc,f,fop)
+    call oft_blag_geval(self%lag_rep,cell,jc,f,rop,gop)
+    rop(3) = rop(2); rop(2) = 0.d0 ! Shuffle 2D gradient to 3D cylindrical
+    val(1) = val(1) + (self%vals(3,j(jc))*rop(2)/pt(1) - self%vals(2,j(jc))*rop(3))
+    val(2) = val(2) + (self%vals(1,j(jc))*rop(3) - self%vals(3,j(jc))*rop(1))
+    val(3) = val(3) + (pt(1)*self%vals(2,j(jc))*rop(1) + self%vals(2,j(jc))*fop - self%vals(1,j(jc))*rop(2))/pt(1)
+  end do
+ELSE
+  val=0.d0
+  do jc=1,self%lag_rep%nce
+    call oft_blag_geval(self%lag_rep,cell,jc,f,rop,gop)
+    val(1) = val(1) + (self%vals(3,j(jc))*rop(2) - self%vals(2,j(jc))*rop(3))
+    val(2) = val(2) + (self%vals(1,j(jc))*rop(3) - self%vals(3,j(jc))*rop(1))
+    val(3) = val(3) + (self%vals(2,j(jc))*rop(1) - self%vals(1,j(jc))*rop(2))
+  end do
+END IF
+deallocate(j)
+DEBUG_STACK_POP
+end subroutine lag_bvcinterp
+!------------------------------------------------------------------------------
 !> Zero a surface Lagrange scalar field at all boundary nodes
 !------------------------------------------------------------------------------
 subroutine zerob_apply(self,a)

src/fem/blag_operators.F90

@@ -237,7 +237,7 @@ MODULE oft_gs
   CLASS(oft_matrix), POINTER :: mop => NULL() !< Lagrange FE mass matrix
   CLASS(oft_matrix), POINTER :: mop_axis => NULL() !< Lagrange FE mass matrix with Dirichlet BCs on axis
   CLASS(oft_bmesh), POINTER :: mesh => NULL() !< Mesh
-  CLASS(oft_scalar_bfem), POINTER :: fe_rep => NULL() !< Lagrange FE representation
+  TYPE(oft_scalar_bfem), POINTER :: fe_rep => NULL() !< Lagrange FE representation
   TYPE(oft_ml_fem_type), POINTER :: ML_fe_rep => NULL() !< Multi-level Lagrange FE representation (only top level used)
   TYPE(oft_blag_zerob), POINTER :: zerob_bc => NULL() !< BC object for zeroing boundary nodes
   TYPE(oft_blag_zerogrnd), POINTER :: zerogrnd_bc => NULL() !< BC object for zeroing grounding node(s)
@@ -4720,7 +4720,7 @@ subroutine gs_prof_interp_apply(self,cell,f,gop,val)
     ELSE
       val(1)=self%equil%psiscale*self%equil%I%f_offset
     END IF
-  CASE(3)
+  CASE(3) ! Isotropic pressure or perpendicular component for anisotropic pressure
     CALL self%psi_eval%interp(cell,f,gop,psitmp)
     IF(in_plasma.AND.(psitmp(1)>self%equil%plasma_bounds(1)))THEN
       ! Handle anisotropic pressure
@@ -4741,6 +4741,19 @@ subroutine gs_prof_interp_apply(self,cell,f,gop,val)
     ELSE
       val(1)=0.d0
     END IF
+  CASE(5) ! Isotropic pressure or parallel component for anisotropic pressure
+    CALL self%psi_eval%interp(cell,f,gop,psitmp)
+    IF(in_plasma.AND.(psitmp(1)>self%equil%plasma_bounds(1)))THEN
+      ! Handle anisotropic pressure
+      IF(ASSOCIATED(self%equil%P_ani))THEN
+        CALL self%equil%P_ani%interp(cell,f,gop,pani)
+        val(1) = (self%equil%psiscale**2)*self%equil%p_scale*self%equil%P%F(psitmp(1))*pani(1)/mu0
+      ELSE
+        val(1) = (self%equil%psiscale**2)*self%equil%p_scale*self%equil%P%F(psitmp(1))/mu0
+      END IF
+    ELSE
+      val(1)=0.d0
+    END IF
   CASE DEFAULT
     CALL oft_abort('Unknown field mode','gs_prof_interp_apply',__FILE__)
 END SELECT

src/physics/grad_shaf.F90

@@ -27,7 +27,7 @@ MODULE oft_gs_util
 USE tracing_2d, ONLY: active_tracer, tracinginv_fs, set_tracer
 USE mhd_utils, ONLY: mu0
 USE oft_gs, ONLY: gs_factory, flux_func, gs_dflux, gs_itor_nl, gs_test_bounds, gs_b_interp, &
-  gsinv_interp, gs_psi2r, gs_psi2pt, gs_epsilon, gs_update_bounds
+  gsinv_interp, gs_psi2r, gs_psi2pt, gs_epsilon, gs_update_bounds, gs_bcrosskappa
 USE oft_gs_profiles
 USE grad_shaf_prof_phys, ONLY: create_jphi_ff, jphi_flux_func
 IMPLICIT NONE
@@ -197,10 +197,10 @@ subroutine gs_comp_globals(self,itor,centroid,vol,pvol,dflux,tflux,bp_vol)
 real(8), intent(out) :: dflux !< Diamagnetic flux
 real(8), intent(out) :: tflux !< Contained toroidal flux
 real(8), intent(out) :: bp_vol !< \f$ \int B_p^2 dV \f$
-type(oft_lag_brinterp) :: psi_eval
+type(oft_lag_brinterp) :: psi_eval,bcross_kappa_fun
 type(oft_lag_bginterp) :: psi_geval
 real(8) :: itor_loc,goptmp(3,3),v,psitmp(1),gpsitmp(3)
-real(8) :: pt(3),curr_cent(2),Btor,Bpol(2)
+real(8) :: pt(3),curr_cent(2),Btor,Bpol(2),pani(2),bcross_kappa(1)
 integer(4) :: i,m
 class(oft_bmesh), pointer :: smesh
 type(gs_factory), pointer :: device
@@ -210,6 +210,13 @@ subroutine gs_comp_globals(self,itor,centroid,vol,pvol,dflux,tflux,bp_vol)
 psi_eval%u=>self%psi
 CALL psi_eval%setup(device%fe_rep)
 CALL psi_geval%shared_setup(psi_eval)
+IF(ASSOCIATED(self%P_ani))THEN
+  CALL self%psi%new(bcross_kappa_fun%u)
+  CALL gs_bcrosskappa(self,bcross_kappa_fun%u)
+  CALL bcross_kappa_fun%setup(device%fe_rep)
+ELSE
+  NULLIFY(bcross_kappa_fun%u)
+END IF
 !---
 itor = 0.d0
 centroid = 0.d0
@@ -218,7 +225,7 @@ subroutine gs_comp_globals(self,itor,centroid,vol,pvol,dflux,tflux,bp_vol)
 dflux = 0.d0
 tflux = 0.d0
 bp_vol = 0.d0
-!$omp parallel do private(m,goptmp,v,psitmp,gpsitmp,pt,itor_loc,Btor,Bpol) &
+!$omp parallel do private(m,goptmp,v,psitmp,gpsitmp,pt,itor_loc,Btor,Bpol,bcross_kappa,pani) &
 !$omp reduction(+:itor) reduction(+:centroid) reduction(+:pvol) reduction(+:vol) reduction(+:dflux) &
 !$omp reduction(+:tflux) reduction(+:bp_vol)
 do i=1,smesh%nc
@@ -231,11 +238,17 @@ subroutine gs_comp_globals(self,itor,centroid,vol,pvol,dflux,tflux,bp_vol)
     !---Compute Magnetic Field
     IF(gs_test_bounds(self,pt))THEN
       IF(self%mode==0)THEN
-        itor_loc = (self%p_scale*pt(1)*self%P%Fp(psitmp(1)) &
-        + self%I%Fp(psitmp(1))*((self%ffp_scale**2)*self%I%f(psitmp(1))+self%ffp_scale*self%I%f_offset)/(pt(1)+gs_epsilon))
+        itor_loc = self%I%Fp(psitmp(1))*((self%ffp_scale**2)*self%I%f(psitmp(1))+self%ffp_scale*self%I%f_offset)/(pt(1)+gs_epsilon)
       ELSE
-        itor_loc = (self%p_scale*pt(1)*self%P%Fp(psitmp(1)) &
-        + .5d0*self%ffp_scale*self%I%Fp(psitmp(1))/(pt(1)+gs_epsilon))
+        itor_loc = 0.5d0*self%ffp_scale*self%I%Fp(psitmp(1))/(pt(1)+gs_epsilon)
+      END IF
+      ! Handle anisotropic pressure
+      IF(ASSOCIATED(self%P_ani))THEN
+        CALL self%P_ani%interp(i,device%fe_rep%quad%pts(:,m),goptmp,pani)
+        CALL bcross_kappa_fun%interp(i,device%fe_rep%quad%pts(:,m),goptmp,bcross_kappa)
+        itor_loc = itor_loc + self%p_scale*pt(1)*(self%P%fp(psitmp(1))*pani(2)+self%P%f(psitmp(1))*(pani(1)-pani(2))*bcross_kappa(1))
+      ELSE
+        itor_loc = itor_loc + self%p_scale*pt(1)*self%P%Fp(psitmp(1))
       END IF
       itor = itor + itor_loc*v*device%fe_rep%quad%wts(m)
       centroid = centroid + itor_loc*pt(1:2)*v*device%fe_rep%quad%wts(m)
@@ -271,6 +284,11 @@ subroutine gs_comp_globals(self,itor,centroid,vol,pvol,dflux,tflux,bp_vol)
 bp_vol=bp_vol*self%psiscale*self%psiscale
 dflux=dflux*self%psiscale
 tflux=tflux*self%psiscale
+IF(ASSOCIATED(bcross_kappa_fun%u))THEN
+  CALL bcross_kappa_fun%u%delete
+  DEALLOCATE(bcross_kappa_fun%u)
+  CALL bcross_kappa_fun%delete
+END IF
 CALL psi_eval%delete
 CALL psi_geval%delete
 end subroutine gs_comp_globals