diff --git a/src/examples/TokaMaker/ITER/ITER_disruption_forces.ipynb b/src/examples/TokaMaker/ITER/ITER_disruption_forces.ipynb index 0a857034..d7a7b5f1 100644 --- a/src/examples/TokaMaker/ITER/ITER_disruption_forces.ipynb +++ b/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", diff --git a/src/physics/grad_shaf.F90 b/src/physics/grad_shaf.F90 index 1fa95877..2c51d34a 100644 --- a/src/physics/grad_shaf.F90 +++ b/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 diff --git a/src/physics/grad_shaf_td.F90 b/src/physics/grad_shaf_td.F90 index 74c9e6e3..cdf72aa3 100644 --- a/src/physics/grad_shaf_td.F90 +++ b/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 \ No newline at end of file +end module oft_gs_td diff --git a/src/python/OpenFUSIONToolkit/TokaMaker/_core.py b/src/python/OpenFUSIONToolkit/TokaMaker/_core.py index 2ac9b3e0..d9b29bf7 100644 --- a/src/python/OpenFUSIONToolkit/TokaMaker/_core.py +++ b/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):