From b0289ffe8d3931903d056cec406c4b43470e5b93 Mon Sep 17 00:00:00 2001 From: DeIonizedPlasma Date: Tue, 15 Jul 2025 22:33:19 -0400 Subject: [PATCH 1/7] fixed crash due to missing empty dictionary in blank call to set_coil_currents --- src/examples/TokaMaker/ITER/ITER_disruption_forces.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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", From edf40520dc71afd5ac4a2cf28038daba63808e9d Mon Sep 17 00:00:00 2001 From: DeIonizedPlasma Date: Thu, 24 Jul 2025 23:20:36 -0400 Subject: [PATCH 2/7] Added vcoil plotting to plot_eddy --- .../OpenFUSIONToolkit/TokaMaker/_core.py | 26 +++++++++++++++++++ 1 file changed, 26 insertions(+) diff --git a/src/python/OpenFUSIONToolkit/TokaMaker/_core.py b/src/python/OpenFUSIONToolkit/TokaMaker/_core.py index 3e2f8b51..99441486 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 @@ -180,6 +182,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 @@ -634,6 +637,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 @@ -1615,6 +1620,27 @@ 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) + + if len(self._vcoils.keys())>0: + coil_currents, _ = self.get_coil_currents() + # 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'] + mask = numpy.logical_or(mask,mask_tmp) + J = coil_currents[coil_name] + if (coil_id:=coil_obj["id"]) in self.dist_coils.keys(): + face_currents = numpy.mean(self.dist_coils[coil_id][self.lc],axis=1) + if cell_centered: + mesh_currents[mask_tmp] = face_currents*J + else: + curr[self.lc][mask_tmp] = self.dist_coils[coil_id][self.lc][mask_tmp]*J + else: + if cell_centered: + mesh_currents[mask_tmp] = J + curr[self.lc][mask_tmp] = J + if cell_centered: return mask, mesh_currents else: From b4a40979c05a131cc46207b897db044b99dfb44b Mon Sep 17 00:00:00 2001 From: DeIonizedPlasma Date: Fri, 25 Jul 2025 17:29:41 -0400 Subject: [PATCH 3/7] Fixed overcomplicated implementation of eddy vcoil plotting --- src/python/OpenFUSIONToolkit/TokaMaker/_core.py | 15 ++------------- 1 file changed, 2 insertions(+), 13 deletions(-) diff --git a/src/python/OpenFUSIONToolkit/TokaMaker/_core.py b/src/python/OpenFUSIONToolkit/TokaMaker/_core.py index 99441486..00a886d4 100644 --- a/src/python/OpenFUSIONToolkit/TokaMaker/_core.py +++ b/src/python/OpenFUSIONToolkit/TokaMaker/_core.py @@ -1621,25 +1621,14 @@ def get_conductor_currents(self,psi,cell_centered=False): mesh_currents[mask_tmp] = numpy.sum(curr[self.lc[mask_tmp,:]],axis=1)/3.0 mask = numpy.logical_or(mask,mask_tmp) - if len(self._vcoils.keys())>0: - coil_currents, _ = self.get_coil_currents() # 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) - J = coil_currents[coil_name] - if (coil_id:=coil_obj["id"]) in self.dist_coils.keys(): - face_currents = numpy.mean(self.dist_coils[coil_id][self.lc],axis=1) - if cell_centered: - mesh_currents[mask_tmp] = face_currents*J - else: - curr[self.lc][mask_tmp] = self.dist_coils[coil_id][self.lc][mask_tmp]*J - else: - if cell_centered: - mesh_currents[mask_tmp] = J - curr[self.lc][mask_tmp] = J if cell_centered: return mask, mesh_currents From f61aa100b0e7d469652d604dd4edea2c23e5b477 Mon Sep 17 00:00:00 2001 From: DeIonizedPlasma Date: Fri, 25 Jul 2025 19:13:03 -0400 Subject: [PATCH 4/7] Added dpsi support for vcoil eddy plots --- src/python/OpenFUSIONToolkit/TokaMaker/_core.py | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/python/OpenFUSIONToolkit/TokaMaker/_core.py b/src/python/OpenFUSIONToolkit/TokaMaker/_core.py index 00a886d4..1f04d1a1 100644 --- a/src/python/OpenFUSIONToolkit/TokaMaker/_core.py +++ b/src/python/OpenFUSIONToolkit/TokaMaker/_core.py @@ -1668,6 +1668,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): From eb8dfd0ddad8e34c0a71137494d6aa65ffa80604 Mon Sep 17 00:00:00 2001 From: DeIonizedPlasma Date: Mon, 15 Sep 2025 11:42:20 -0400 Subject: [PATCH 5/7] Updated units on mygs.get_jtor_plasma() to be SI by undoing division by mu0 --- src/python/OpenFUSIONToolkit/TokaMaker/_core.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/python/OpenFUSIONToolkit/TokaMaker/_core.py b/src/python/OpenFUSIONToolkit/TokaMaker/_core.py index 9f240c90..d9b29bf7 100644 --- a/src/python/OpenFUSIONToolkit/TokaMaker/_core.py +++ b/src/python/OpenFUSIONToolkit/TokaMaker/_core.py @@ -1103,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 From d59ce6a1cc6508b91a03ac8199228213283300fa Mon Sep 17 00:00:00 2001 From: DeIonizedPlasma Date: Thu, 25 Sep 2025 19:50:12 -0400 Subject: [PATCH 6/7] added new variable to avoid reuse for clarity --- src/physics/grad_shaf_td.F90 | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) 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 From 999ef951b746dfb5ec1647dc4fa814b76a1b966e Mon Sep 17 00:00:00 2001 From: DeIonizedPlasma Date: Thu, 25 Sep 2025 19:51:02 -0400 Subject: [PATCH 7/7] Fixed bug when using linear solver with distribution specified on voltage coils --- src/physics/grad_shaf.F90 | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) 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