Handle Vcoils in plot_eddy

 - Fixed bug when using linear solver with distribution specified on voltage coils
 - Correct units in `get_jtor_plasma()`

Author: hansec

src/examples/TokaMaker/ITER/ITER_disruption_forces.ipynb

@@ -402,7 +402,7 @@
     "sim_time = [0.0]\n",
     "mygs.settings.pm=False\n",
     "mygs.update_settings()\n",
-    "mygs.set_coil_currents()\n",
+    "mygs.set_coil_currents({})\n",
     "for i in range(100):\n",
     "    if i == 40:\n",
     "        dt = CQ_time/5.0\n",

src/physics/grad_shaf.F90

@@ -1806,12 +1806,14 @@ subroutine gs_wall_source(self,dpsi_dt,b)
 CLASS(oft_vector), intent(inout) :: dpsi_dt !< \f$ \psi \f$ at start of step
 CLASS(oft_vector), intent(inout) :: b !< Resulting source field
 real(r8), pointer, dimension(:) :: btmp,btmp_coils
-real(8) :: psitmp,goptmp(3,3),det,pt(3),v,ffp(3),t1,nturns
+real(8) :: psitmp,goptmp(3,3),det,pt(3),v,ffp(3),t1,nturns,cond_norm
 real(8), allocatable :: rhs_loc(:),cond_fac(:),rop(:),vcache(:),eta_reg(:),reg_source(:)
 real(8), pointer :: psi_vals(:),coil_vals(:)
-integer(4) :: j,m,l,k,kk,iCond
+integer(4) :: j,m,l,k,kk,iCond,jc
 integer(4), allocatable :: j_lag(:)
 logical :: curved
+CLASS(oft_scalar_bfem), POINTER :: lag_rep
+lag_rep=>self%fe_rep
 ! t1=omp_get_wtime()
 !---
 NULLIFY(btmp,btmp_coils,psi_vals,coil_vals)
@@ -1839,7 +1841,7 @@ subroutine gs_wall_source(self,dpsi_dt,b)
   END DO
 END IF
 !---
-!$omp parallel private(rhs_loc,j_lag,curved,goptmp,v,m,det,pt,psitmp,l,rop,nturns)
+!$omp parallel private(rhs_loc,j_lag,curved,goptmp,v,m,det,pt,psitmp,l,rop,nturns,cond_norm,jc)
 allocate(rhs_loc(self%fe_rep%nce+self%ncoils))
 allocate(rop(self%fe_rep%nce))
 allocate(j_lag(self%fe_rep%nce))
@@ -1863,7 +1865,11 @@ subroutine gs_wall_source(self,dpsi_dt,b)
       IF(self%Rcoils(l)<=0.d0)CYCLE
       nturns = self%coil_nturns(self%fe_rep%mesh%reg(j),l)
       IF(ABS(nturns)>1.d-8)THEN
-        rhs_loc(self%fe_rep%nce+l)=rhs_loc(self%fe_rep%nce+l)+mu0*2.d0*pi*psitmp*nturns*det
+        cond_norm=0.d0
+        do jc=1,lag_rep%nce ! Loop over degrees of freedom
+          cond_norm = cond_norm + self%dist_coil(j_lag(jc),l)*rop(jc)
+        end do
+        rhs_loc(self%fe_rep%nce+l)=rhs_loc(self%fe_rep%nce+l)+mu0*2.d0*pi*psitmp*nturns*cond_norm*det
       END IF
     END DO
     IF(eta_reg(self%fe_rep%mesh%reg(j))<0.d0)CYCLE

src/physics/grad_shaf_td.F90

@@ -482,7 +482,7 @@ subroutine apply_rhs(self,a,b)
 class(oft_vector), intent(inout) :: b !< Result of metric function
 integer(i4) :: i,m,jr,jc
 integer(i4), allocatable :: j(:)
-real(r8) :: vol,det,goptmp(3,3),elapsed_time,pt(3),eta_tmp,psi_tmp,eta_source,psi_lim,psi_max
+real(r8) :: vol,det,goptmp(3,3),elapsed_time,pt(3),eta_tmp,psi_tmp,eta_source,psi_lim,psi_max,cond_norm
 real(8) :: max_tmp,lim_tmp,vcont_val,nturns
 real(r8), allocatable :: rop(:),gop(:,:),lop(:,:),vals_loc(:),reg_source(:)
 real(r8), pointer, dimension(:) :: pol_vals,rhs_vals
@@ -528,7 +528,7 @@ subroutine apply_rhs(self,a,b)
 !------------------------------------------------------------------------------
 ! Operator integration
 !------------------------------------------------------------------------------
-!$omp parallel private(j,vals_loc,rop,gop,det,curved,goptmp,m,vol,jr,jc,pt,eta_tmp,psi_tmp,eta_source,nturns)
+!$omp parallel private(j,vals_loc,rop,gop,det,curved,goptmp,m,vol,jr,jc,pt,eta_tmp,psi_tmp,eta_source,nturns,cond_norm)
 allocate(j(lag_rep%nce),vals_loc(lag_rep%nce+self%gs_eq%ncoils)) ! Local DOF and matrix indices
 allocate(rop(lag_rep%nce),gop(3,lag_rep%nce)) ! Reconstructed gradient operator
 !$omp do schedule(static,1)
@@ -561,11 +561,11 @@ subroutine apply_rhs(self,a,b)
             IF(self%gs_eq%Rcoils(jr)<=0.d0)CYCLE
             nturns = self%gs_eq%coil_nturns(mesh%reg(i),jr)
             IF(ABS(nturns)>1.d-8)THEN
-                eta_source=0.d0
+                cond_norm=0.d0
                 do jc=1,lag_rep%nce ! Loop over degrees of freedom
-                    eta_source = eta_source + self%gs_eq%dist_coil(j(jc),jr)*rop(jc)
+                    cond_norm = cond_norm + self%gs_eq%dist_coil(j(jc),jr)*rop(jc)
                 end do
-                vals_loc(lag_rep%nce+jr)=vals_loc(lag_rep%nce+jr)+mu0*2.d0*pi*psi_tmp*nturns*eta_source*det
+                vals_loc(lag_rep%nce+jr)=vals_loc(lag_rep%nce+jr)+mu0*2.d0*pi*psi_tmp*nturns*cond_norm*det
             END IF
         END DO
     end do
@@ -1438,4 +1438,4 @@ subroutine apply_wrap(self,a,b)
 DEALLOCATE(tmp_vec)
 DEBUG_STACK_POP
 end subroutine apply_wrap
-end module oft_gs_td
+end module oft_gs_td

src/python/OpenFUSIONToolkit/TokaMaker/_core.py

@@ -108,6 +108,8 @@ def __init__(self,OFT_env):
         self.coil_sets = {}
         ## Virtual coils, if present (currently only `'#VSC'`)
         self._virtual_coils = {'#VSC': {'id': -1 ,'facs': {}}}
+        ## Voltage coils dictionary. Currently only used for plotting on python side.
+        self._vcoils = {}
         ## Coil set names in order of id number
         self.coil_set_names = []
         ## Distribution coils, only (currently) saved for plotting utility
@@ -182,6 +184,7 @@ def reset(self):
         self._cond_dict = {}
         self._vac_dict = {}
         self._coil_dict = {}
+        self._vcoils = {}
         self.coil_sets = {}
         self._virtual_coils = {'#VSC': {'id': -1 ,'facs': {}}}
         self._F0 = 0.0
@@ -644,6 +647,8 @@ def set_vcoils(self,coil_resistivities):
         tokamaker_set_vcoil(self._tMaker_ptr,res_array,error_string)
         if error_string.value != b'':
             raise Exception(error_string.value)
+        #Merge dicts and overwrite with new values where necessary. Only used for plotting.
+        self._vcoils = self._vcoils | coil_resistivities
 
     def init_psi(self, r0=-1.0, z0=0.0, a=0.0, kappa=0.0, delta=0.0, curr_source=None):
         r'''! Initialize \f$\psi\f$ using uniform current distributions
@@ -1098,7 +1103,7 @@ def get_jtor_plasma(self):
         tokamaker_get_jtor(self._tMaker_ptr,curr,error_string)
         if error_string.value != b'':
             raise Exception(error_string.value)
-        return curr/mu0
+        return curr
 
     def get_psi(self,normalized=True):
         r'''! Get poloidal flux values on node points
@@ -1625,6 +1630,16 @@ def get_conductor_currents(self,psi,cell_centered=False):
                 if cell_centered:
                     mesh_currents[mask_tmp] = numpy.sum(curr[self.lc[mask_tmp,:]],axis=1)/3.0
                 mask = numpy.logical_or(mask,mask_tmp)
+                
+        # Treat vcoils as conductors when looking at induced currents
+        for coil_name, coil_obj in self.coil_sets.items():
+            if coil_name in self._vcoils.keys():
+                for sub_coil in coil_obj["sub_coils"]:
+                    mask_tmp = self.reg == sub_coil['reg_id']
+                    if cell_centered:
+                        mesh_currents[mask_tmp] = numpy.mean(curr[self.lc[mask_tmp]],axis=1)
+                    mask = numpy.logical_or(mask,mask_tmp)
+        
         if cell_centered:
             return mask, mesh_currents
         else:
@@ -1663,6 +1678,16 @@ def get_conductor_source(self,dpsi_dt):
                 if cond_reg.get('noncontinuous',False):
                     mesh_currents[mask_tmp] -= (mesh_currents[mask_tmp]*area[mask_tmp]).sum()/area[mask_tmp].sum()
                 mask = numpy.logical_or(mask,mask_tmp)
+
+        # Treat vcoils as conductors when looking at induced currents
+        for coil_name, coil_obj in self.coil_sets.items():
+            if coil_name in self._vcoils.keys():
+                for sub_coil in coil_obj["sub_coils"]:
+                    mask_tmp = self.reg == sub_coil['reg_id']
+                    field_tmp = -dpsi_dt/self._vcoils[coil_name]
+                    mesh_currents[mask_tmp] = numpy.mean(field_tmp[self.lc[mask_tmp]],axis=1)
+                    mask = numpy.logical_or(mask,mask_tmp)
+
         return mask, mesh_currents
 
     def plot_eddy(self,fig,ax,psi=None,dpsi_dt=None,nlevels=40,colormap='jet',clabel=r'$J_w$ [$A/m^2$]',symmap=False):