diff --git a/.gitignore b/.gitignore index 8ecd1dacb..77ea7de1d 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,7 @@ # Python artifacts *.pyc __pycache__ +.ipynb_checkpoints # Editor Items .settings diff --git a/src/physics/grad_shaf.F90 b/src/physics/grad_shaf.F90 index 98cf5ced6..45d024c53 100644 --- a/src/physics/grad_shaf.F90 +++ b/src/physics/grad_shaf.F90 @@ -264,6 +264,7 @@ MODULE oft_gs TYPE :: gs_equil LOGICAL :: diverted = .FALSE. !< Equilibrium is diverted? LOGICAL :: has_plasma = .TRUE. !< Solve with plasma? (otherwise vacuum) + LOGICAL :: skip_targets = .FALSE. !< Skip toroidal current target in non-linear solve? INTEGER(i4) :: mode = 0 !< RHS source mode (0 -> F*F', 1 -> F') INTEGER(i4) :: nx_points = 0 !< Number of X-points in current solution INTEGER(i4) :: nregularize = 0 !< Number of regularization terms @@ -279,7 +280,7 @@ MODULE oft_gs REAL(r8) :: mirror_n = -1.d0 !< Anisotropy exponent for mirror pressure profiles REAL(r8) :: mirror_bturn = 0.d0 !< Turning point for mirror pressure profiles REAL(r8) :: mirror_zthroat = 0.d0 !< Mirror peak field point - REAL(r8) :: Itor_target = -1.d0 !< Toroidal current target + REAL(r8) :: Ip_target = -1.d0 !< Toroidal current target REAL(r8) :: estore_target = -1.d0 !< Stored energy target REAL(r8) :: dflux_target = -1.d99 !< Diamagnetic flux target REAL(r8) :: pax_target = -1.d0 !< On-axis pressure target @@ -1052,7 +1053,8 @@ subroutine copy_eq(self,source) self%mirror_n=source%mirror_n self%mirror_bturn=source%mirror_bturn self%mirror_zthroat=source%mirror_zthroat -self%Itor_target=source%Itor_target +self%skip_targets=source%skip_targets +self%Ip_target=source%Ip_target self%estore_target=source%estore_target self%dflux_target=source%dflux_target self%pax_target=source%pax_target @@ -1121,13 +1123,13 @@ subroutine gs_init_psi(self,equil,ierr,r0,a,kappa,delta,curr_source) CALL equil%psi%set(0.d0) CALL gs_vacuum_solve(self,equil%psi,tmp_vec) itor=equil%itor() - IF(equil%Itor_target>0.d0)CALL equil%psi%scale(equil%Itor_target/itor) + IF(equil%Ip_target>0.d0)CALL equil%psi%scale(equil%Ip_target/itor) ELSE IF(PRESENT(curr_source))THEN CALL equil%psi%set(0.d0) CALL tmp_vec%restore_local(curr_source) CALL gs_vacuum_solve(self,equil%psi,tmp_vec) itor=equil%itor() - IF(equil%Itor_target>0.d0)CALL equil%psi%scale(equil%Itor_target/itor) + IF(equil%Ip_target>0.d0)CALL equil%psi%scale(equil%Ip_target/itor) ELSE !---Setup Solver eigsolver%A=>self%dels @@ -1153,9 +1155,9 @@ subroutine gs_init_psi(self,equil,ierr,r0,a,kappa,delta,curr_source) call equil%psi%scale(1.d0/psi_vals(mind)) equil%psimax=1.d0 END IF - IF(equil%Itor_target>0.d0)THEN + IF(equil%Ip_target>0.d0)THEN itor=equil%itor() - CALL equil%psi%scale(equil%Itor_target/itor) + CALL equil%psi%scale(equil%Ip_target/itor) END IF !---Cleanup CALL eigsolver%pre%delete @@ -2123,6 +2125,7 @@ subroutine gs_solve(self,equil,ierr) !--- error_flag=0 self%nl_its=0 +equil%skip_targets=.FALSE. IF(TRIM(self%lu_solver%package)=='none')THEN CALL oft_abort("LU solver required for GS solve","gs_solve",__FILE__) ELSE @@ -2299,12 +2302,17 @@ subroutine gs_solve(self,equil,ierr) param_mat(1,1)=1.d0 param_rhs(1)=0.d0 ELSE - IF(equil%Itor_target>0.d0)THEN - param_mat(1,:)=[itor_ffp,itor_press,0.d0] - param_rhs(1)=equil%Itor_target - ELSE + IF(equil%skip_targets)THEN param_mat(1,1)=1.d0 param_rhs(1)=equil%ffp_scale + ELSE + IF(equil%Ip_target>0.d0)THEN + param_mat(1,:)=[itor_ffp,itor_press,0.d0] + param_rhs(1)=equil%Ip_target + ELSE + param_mat(1,1)=1.d0 + param_rhs(1)=equil%ffp_scale + END IF END IF END IF @@ -2334,42 +2342,49 @@ subroutine gs_solve(self,equil,ierr) param_rhs(2)=equil%p_scale END IF ELSE - IF(equil%R0_target>0.d0)THEN - IF(ALL(self%target_weights>0.d0))THEN - equil%saddle_targets(1:2,equil%saddle_ntargets)=pt - ELSE - ! - psi_geval%u=>psi_vac - CALL psi_geval%setup(self%fe_rep) - CALL psi_geval%interp(cell,f,goptmp,gpsi0) - param_rhs(2)=-gpsi0(1) - ! - psi_geval%u=>psi_vcont - CALL psi_geval%setup(self%fe_rep) - CALL psi_geval%interp(cell,f,goptmp,gpsi0) - psi_geval%u=>psi_ffp - CALL psi_geval%setup(self%fe_rep) - CALL psi_geval%interp(cell,f,goptmp,gpsi1) - psi_geval%u=>psi_press - CALL psi_geval%setup(self%fe_rep) - CALL psi_geval%interp(cell,f,goptmp,gpsi2) - param_mat(2,:)=[gpsi1(1),gpsi2(1),gpsi0(1)] - END IF - ELSE IF(equil%dflux_target>-1.d98)THEN - param_rhs(2)=SIGN(equil%dflux_target**2,equil%dflux_target) - param_mat(2,:)=[SIGN(dflux_ffp**2,dflux_ffp)/equil%ffp_scale,0.d0,0.d0] - ELSE IF(equil%estore_target>0.d0)THEN - param_rhs(2)=equil%estore_target - param_mat(2,2)=estored*3.d0/2.d0 - ELSE IF(equil%pax_target>0.d0)THEN - param_mat(2,2)=equil%P%f(equil%plasma_bounds(2)) - param_rhs(2)=equil%pax_target - ELSE IF(equil%Ip_ratio_target>-1.d98)THEN - param_rhs(2)=0.d0 - param_mat(2,:)=[itor_ffp,-itor_press*equil%Ip_ratio_target,0.d0] - ELSE + IF(equil%skip_targets)THEN + IF(ALL(self%target_weights>0.d0))CALL oft_abort("Soft targets are not compatible with externally-handled targets", & + "gs_solve",__FILE__) param_mat(2,2)=1.d0 param_rhs(2)=equil%p_scale + ELSE + IF(equil%R0_target>0.d0)THEN + IF(ALL(self%target_weights>0.d0))THEN + equil%saddle_targets(1:2,equil%saddle_ntargets)=pt + ELSE + ! + psi_geval%u=>psi_vac + CALL psi_geval%setup(self%fe_rep) + CALL psi_geval%interp(cell,f,goptmp,gpsi0) + param_rhs(2)=-gpsi0(1) + ! + psi_geval%u=>psi_vcont + CALL psi_geval%setup(self%fe_rep) + CALL psi_geval%interp(cell,f,goptmp,gpsi0) + psi_geval%u=>psi_ffp + CALL psi_geval%setup(self%fe_rep) + CALL psi_geval%interp(cell,f,goptmp,gpsi1) + psi_geval%u=>psi_press + CALL psi_geval%setup(self%fe_rep) + CALL psi_geval%interp(cell,f,goptmp,gpsi2) + param_mat(2,:)=[gpsi1(1),gpsi2(1),gpsi0(1)] + END IF + ELSE IF(equil%dflux_target>-1.d98)THEN + param_rhs(2)=SIGN(equil%dflux_target**2,equil%dflux_target) + param_mat(2,:)=[SIGN(dflux_ffp**2,dflux_ffp)/equil%ffp_scale,0.d0,0.d0] + ELSE IF(equil%estore_target>0.d0)THEN + param_rhs(2)=equil%estore_target + param_mat(2,2)=estored*3.d0/2.d0 + ELSE IF(equil%pax_target>0.d0)THEN + param_mat(2,2)=equil%P%f(equil%plasma_bounds(2)) + param_rhs(2)=equil%pax_target + ELSE IF(equil%Ip_ratio_target>-1.d98)THEN + param_rhs(2)=0.d0 + param_mat(2,:)=[itor_ffp,-itor_press*equil%Ip_ratio_target,0.d0] + ELSE + param_mat(2,2)=1.d0 + param_rhs(2)=equil%p_scale + END IF END IF END IF @@ -2563,13 +2578,13 @@ subroutine gs_solve(self,equil,ierr) CALL equil%psi%scale(1.d0/equil%psimax) equil%ffp_scale=SQRT((equil%ffp_scale**2)/equil%psimax) equil%psimax=1.d0 - ! ELSE IF(self%Itor_target>0.d0)THEN + ! ELSE IF(self%Ip_target>0.d0)THEN ! itor=self%itor() ! IF(itor<=0.d0)THEN ! error_flag=-5 ! EXIT ! END IF - ! ffp_scale_prev=itor/self%Itor_target + ! ffp_scale_prev=itor/self%Ip_target ! CALL self%psi%scale(1.d0/ffp_scale_prev) ! self%ffp_scale=SQRT((self%ffp_scale**2)/ffp_scale_prev) END IF @@ -2760,11 +2775,11 @@ subroutine gs_lin_solve(self,equil,adjust_r0,ierr) param_mat=0.d0 param_rhs=0.d0 -IF(equil%Itor_target>0.d0)THEN +IF(equil%Ip_target>0.d0)THEN itor_ffp=equil%itor(psi_ffp) itor_press=equil%itor(psi_press) param_mat(1,:)=[itor_ffp,itor_press,0.d0] - param_rhs(1)=equil%Itor_target + param_rhs(1)=equil%Ip_target ELSE param_mat(1,1)=1.d0 param_rhs(1)=equil%ffp_scale diff --git a/src/physics/grad_shaf_fit.F90 b/src/physics/grad_shaf_fit.F90 index 56ac807af..df775d779 100644 --- a/src/physics/grad_shaf_fit.F90 +++ b/src/physics/grad_shaf_fit.F90 @@ -324,7 +324,7 @@ SUBROUTINE fit_gs(gs,inpath,outpath,fitI,fitP,fit_p_scale,fit_ffp_scale,fitR0,fi !---Count coefficients ncofs=0 IF(gs%device%free)THEN - IF(fit_FFPscale.OR.(gs_active%Itor_target>0.d0))ncofs = ncofs+1 + IF(fit_FFPscale.OR.(gs_active%Ip_target>0.d0))ncofs = ncofs+1 ELSE ncofs=ncofs+1 IF(fit_FFPscale)CALL oft_abort('Lambda cannot be fit in fixed boundary mode.', & @@ -349,14 +349,14 @@ SUBROUTINE fit_gs(gs,inpath,outpath,fitI,fitP,fit_p_scale,fit_ffp_scale,fitR0,fi IF(fit_FFPscale)THEN offset=1 cofs(1)=gs_active%ffp_scale - ELSE IF(gs_active%Itor_target>0.d0)THEN + ELSE IF(gs_active%Ip_target>0.d0)THEN offset=1 - cofs(1)=gs_active%Itor_target + cofs(1)=gs_active%Ip_target END IF ELSE offset=1 - IF(gs_active%Itor_target>0.d0)THEN - cofs(1)=gs_active%Itor_target + IF(gs_active%Ip_target>0.d0)THEN + cofs(1)=gs_active%Ip_target ELSE cofs(1)=gs_active%psiscale cofs_scale(1)=1.d0/ABS(gs_active%psiscale) @@ -466,7 +466,7 @@ SUBROUTINE fit_gs(gs,inpath,outpath,fitI,fitP,fit_p_scale,fit_ffp_scale,fitR0,fi geval_count=0 !---Initialize CALL fit_error(ncons,ncofs,cofs,error,info) - ! gs_active%Itor_target=-1.d0 + ! gs_active%Ip_target=-1.d0 ! IF(fit_FFPscale)cofs(1)=gs_active%ffp_scale IF(gs_active%device%ierr<0)CALL oft_abort('Initial equilibrium solve failed to converge','fit_gs',__FILE__) !--- @@ -635,7 +635,7 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) ALLOCATE(cofs_in(n)) cofs_in=cofs ffp_scale_in=gs_active%ffp_scale - ip_target_in=gs_active%Itor_target + ip_target_in=gs_active%Ip_target p_scale_in=gs_active%p_scale estore_target_in=gs_active%estore_target ! bounds_in=gs_active%spatial_bounds @@ -669,14 +669,14 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) IF(fit_FFPscale)THEN offset=1 gs_active%ffp_scale=cofs(1) - ELSE IF(gs_active%Itor_target>0.d0)THEN + ELSE IF(gs_active%Ip_target>0.d0)THEN offset=1 - gs_active%Itor_target=cofs(1) + gs_active%Ip_target=cofs(1) END IF ELSE offset=1 - IF(gs_active%Itor_target>0.d0)THEN - gs_active%Itor_target=cofs(1) + IF(gs_active%Ip_target>0.d0)THEN + gs_active%Ip_target=cofs(1) ELSE gs_active%psiscale=cofs(1) END IF @@ -741,7 +741,7 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) IF(fit_F0)gs_active%I%f_offset = cofs(offset+1) ! !---Centering ! ffp_scale_in=gs_active%ffp_scale - ! ip_target_in=gs_active%Itor_target + ! ip_target_in=gs_active%Ip_target ! p_scale_in=gs_active%p_scale ! estore_target_in=gs_active%estore_target ! bounds_in=gs_active%spatial_bounds @@ -777,7 +777,7 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) CALL psi_best%add(0.d0,1.d0,gs_active%psi) END IF ffp_scale_in=gs_active%ffp_scale - ip_target_in=gs_active%Itor_target + ip_target_in=gs_active%Ip_target p_scale_in=gs_active%p_scale estore_target_in=gs_active%estore_target ! bounds_in=gs_active%spatial_bounds @@ -801,13 +801,13 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) IF(gs_active%device%free)THEN IF(fit_FFPscale)THEN offset=offset+1 - ELSE IF(gs_active%Itor_target>0.d0)THEN - WRITE(*,'(2A,ES11.3)')oft_indent,'Itor_target =',gs_active%Itor_target/mu0 + ELSE IF(gs_active%Ip_target>0.d0)THEN + WRITE(*,'(2A,ES11.3)')oft_indent,'Ip_target =',gs_active%Ip_target/mu0 offset=offset+1 END IF ELSE - IF(gs_active%Itor_target>0.d0)THEN - WRITE(*,'(2A,ES11.3)')oft_indent,'Itor_target =',gs_active%Itor_target/mu0 + IF(gs_active%Ip_target>0.d0)THEN + WRITE(*,'(2A,ES11.3)')oft_indent,'Ip_target =',gs_active%Ip_target/mu0 ELSE WRITE(*,'(2A,ES11.3)')oft_indent,'Psi_scale =',gs_active%psiscale END IF @@ -890,18 +890,18 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) jac_mat(:,offset+1)=(jac_mat(:,offset+1)-err)/dx gs_active%ffp_scale=cofs(offset+1) offset=1 - ELSE IF(gs_active%Itor_target>0.d0)THEN + ELSE IF(gs_active%Ip_target>0.d0)THEN CALL reset_eq dx = dxi/cofs_scale(offset+1) - gs_active%Itor_target=cofs(offset+1) + dx + gs_active%Ip_target=cofs(offset+1) + dx CALL run_err(.FALSE.,jac_mat(:,offset+1),m,ierr) jac_mat(:,offset+1)=(jac_mat(:,offset+1)-err)/dx - gs_active%Itor_target=cofs(offset+1) + gs_active%Ip_target=cofs(offset+1) offset=1 END IF ELSE - IF(gs_active%Itor_target>0.d0)THEN - gs_active%Itor_target=cofs(1) + IF(gs_active%Ip_target>0.d0)THEN + gs_active%Ip_target=cofs(1) ELSE CALL reset_eq dx = dxi/cofs_scale(offset+1) @@ -1039,7 +1039,7 @@ SUBROUTINE fit_error_grad(m,n,cofs,err,jac_mat,ldjac_mat,iflag) ! SUBROUTINE reset_eq gs_active%ffp_scale=ffp_scale_in -gs_active%Itor_target=ip_target_in +gs_active%Ip_target=ip_target_in gs_active%p_scale=p_scale_in gs_active%estore_target=estore_target_in ! gs_active%spatial_bounds=bounds_in diff --git a/src/physics/grad_shaf_prof_phys.F90 b/src/physics/grad_shaf_prof_phys.F90 index 3dbbbf520..cc98d0a43 100644 --- a/src/physics/grad_shaf_prof_phys.F90 +++ b/src/physics/grad_shaf_prof_phys.F90 @@ -58,6 +58,7 @@ module grad_shaf_prof_phys type, extends(linterp_flux_func) :: jphi_flux_func integer(4) :: ngeom = 50 !< Number of points in psi for , <1/R> evaluation real(8) :: j0 = 0.d0 !< LCFS Jphi value + real(8) :: norm_last = 1.d0 !< Last Jphi normalization factor (for Ip target) real(8), pointer, dimension(:) :: jphi => NULL() !< Jphi(psi) profile values contains !> Needs docs @@ -500,6 +501,7 @@ subroutine jphi_copy(self,new) new%f_offset=self%f_offset new%ngeom = self%ngeom new%j0 = self%j0 + new%norm_last = self%norm_last ALLOCATE(new%jphi, SOURCE=self%jphi) END SELECT end subroutine jphi_copy @@ -521,12 +523,12 @@ subroutine jphi_update(self,gseq) class(jphi_flux_func), intent(inout) :: self class(gs_equil), intent(inout) :: gseq INTEGER(i4) :: i -REAL(r8) :: jphi_norm,pscale,pprime +REAL(r8) :: jphi_norm,pscale,pprime,dnorm REAL(r8), ALLOCATABLE :: ravgs(:,:),qtmp(:),psi_q(:) type(spline_type) :: R_spline self%plasma_bounds=gseq%plasma_bounds IF(gseq%mode/=1)CALL oft_abort("Jphi profile requires (F^2)' formulation","jphi_update",__FILE__) -! IF(gseq%Itor_target<0.d0)CALL oft_abort("Jphi profile requires Ip target","jphi_update",__FILE__) +IF(gseq%Ip_target<0.d0)CALL oft_abort("Jphi profile requires Ip target","jphi_update",__FILE__) IF(gseq%pax_target<0.d0)CALL oft_abort("Jphi profile requires Pax target","jphi_update",__FILE__) !---Get updated flux surface geometry for Jphi -> F*F' mapping ALLOCATE(ravgs(self%ngeom+1,3),psi_q(self%ngeom+1),qtmp(self%ngeom+1)) @@ -550,16 +552,24 @@ subroutine jphi_update(self,gseq) CALL spline_fit(R_spline,"extrap") DEALLOCATE(ravgs,psi_q,qtmp) !---Update jphi normalization to match Ip target -ALLOCATE(qtmp(self%npsi)) -DO i=1,self%npsi - CALL spline_eval(R_spline,self%x(i),0) - qtmp(i) = R_spline%f(1)*R_spline%f(2) -END DO -CALL gs_flux_int(gseq,self%x,self%jphi/qtmp,self%npsi,jphi_norm) -jphi_norm=ABS(gseq%Itor_target)/jphi_norm -DEALLOCATE(qtmp) -!---Get pressure profile CALL gseq%P%update(gseq) ! Make sure pressure profile is up to date with EQ +IF(gseq%skip_targets)THEN + CALL gs_itor_nl(gseq,jphi_norm) + dnorm=gseq%Ip_target/jphi_norm + jphi_norm=(1.d0+dnorm)*self%norm_last/2.d0 + self%norm_last=jphi_norm +ELSE + ALLOCATE(qtmp(self%npsi)) + DO i=1,self%npsi + CALL spline_eval(R_spline,self%x(i),0) + qtmp(i) = R_spline%f(1)*R_spline%f(2) + END DO + CALL gs_flux_int(gseq,self%x,self%jphi/qtmp,self%npsi,jphi_norm) + DEALLOCATE(qtmp) + jphi_norm=ABS(gseq%Ip_target)/jphi_norm + self%norm_last=jphi_norm +END IF +!---Get pressure profile IF(ASSOCIATED(gseq%P_ani))CALL oft_abort('Jphi profiles do not support anistopic pressure','jphi_update',__FILE__) !CALL gseq%P_ani%update(gseq) pscale=gseq%P%f(gseq%plasma_bounds(2)) pscale=gseq%pax_target/pscale @@ -573,8 +583,10 @@ subroutine jphi_update(self,gseq) self%yp(i) = 2.d0*(self%jphi(i)*jphi_norm - R_spline%f(1)*pprime*pscale)/R_spline%f(2) END DO ! Disable Ip matching and fix F*F' scale (matching is done here instead) -IF(gseq%Itor_target>0.d0)gseq%Itor_target=-gseq%Itor_target +gseq%skip_targets=.TRUE. +! IF(gseq%Ip_target>0.d0)gseq%Ip_target=-gseq%Ip_target gseq%ffp_scale=1.d0 +gseq%p_scale=pscale !---Clean up CALL spline_dealloc(R_spline) i=self%set_cofs(self%yp) diff --git a/src/physics/grad_shaf_td.F90 b/src/physics/grad_shaf_td.F90 index 1031ca9b8..d21224148 100644 --- a/src/physics/grad_shaf_td.F90 +++ b/src/physics/grad_shaf_td.F90 @@ -603,7 +603,7 @@ subroutine setup_mfop(self,eq_in) DEBUG_STACK_PUSH self%gs_equil=>eq_in self%gs_device=>eq_in%device -self%ip_target=self%gs_equil%Itor_target +self%ip_target=self%gs_equil%Ip_target self%ip_ratio_target=self%gs_equil%Ip_ratio_target self%f_scale=self%gs_equil%ffp_scale self%p_scale=self%gs_equil%p_scale @@ -632,7 +632,7 @@ subroutine update_mfop(self) self%F=>self%gs_equil%I self%P=>self%gs_equil%P ! Update current target and sync scale factors -self%ip_target=self%gs_equil%Itor_target +self%ip_target=self%gs_equil%Ip_target self%ip_ratio_target=self%gs_equil%Ip_ratio_target self%f_scale=self%gs_equil%ffp_scale self%p_scale=self%gs_equil%p_scale diff --git a/src/physics/grad_shaf_util.F90 b/src/physics/grad_shaf_util.F90 index 2c3e43574..2cfc141ea 100644 --- a/src/physics/grad_shaf_util.F90 +++ b/src/physics/grad_shaf_util.F90 @@ -391,7 +391,7 @@ subroutine gs_save_tokamaker(self,filename) CALL hdf5_write(self%lim_point,filename,'tokamaker/LIM_POINT') CALL hdf5_write(self%diverted,filename,'tokamaker/DIVERTED') !--- -IF(self%Itor_target>0.d0)CALL hdf5_write(self%Itor_target,filename,'tokamaker/IP_TARGET') +IF(self%Ip_target>0.d0)CALL hdf5_write(self%Ip_target,filename,'tokamaker/IP_TARGET') IF(self%Ip_ratio_target>-1.d98)CALL hdf5_write(self%Ip_ratio_target,filename,'tokamaker/IP_RATIO_TARGET') IF(self%R0_target>0.d0)CALL hdf5_write(self%R0_target,filename,'tokamaker/R0_TARGET') IF(self%Z0_target>-1.d98)CALL hdf5_write(self%Z0_target,filename,'tokamaker/Z0_TARGET') @@ -612,7 +612,7 @@ subroutine gs_load_tokamaker(self,filename,error_string) END IF !---Load targets and coil currents IF(hdf5_field_exist(filename,'tokamaker/IP_TARGET'))THEN - CALL hdf5_read(self%Itor_target,filename,'tokamaker/IP_TARGET',success=success) + CALL hdf5_read(self%Ip_target,filename,'tokamaker/IP_TARGET',success=success) IF(.NOT.success)THEN error_string='Failed to read Ip target.' RETURN diff --git a/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py b/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py index 4727c7e75..f3eef9060 100644 --- a/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py +++ b/src/python/OpenFUSIONToolkit/TokaMaker/bootstrap.py @@ -207,15 +207,31 @@ def analyze_bootstrap_edge_spike(psi_N, j_bootstrap, diagnostic_plots=False): psi_edge = psi_N[edge_mask] j_edge = j_bootstrap[edge_mask] - # Find peak in the edge region - peaks, properties = find_peaks(j_edge, height=0.)#0.5*numpy.max(j_edge)) + # Guard against an empty edge slice (psi_N may not reach the edge region): + # numpy.max on a zero-size array raises, where the original code fell through + # to find_peaks returning no peaks. + if j_edge.size == 0: + print("No edge region (psi_N >= 0.7) in profile") + return None + + # Find peak in the edge region (prominence filter suppresses noise + # oscillations that can be misidentified for small bootstrap spikes) + j_edge_max = numpy.max(j_edge) + min_prominence = 0.05 * j_edge_max if j_edge_max > 0 else 0.0 + peaks, properties = find_peaks(j_edge, height=0., prominence=min_prominence) if len(peaks) == 0: print("No clear peak found in the edge region") return None - # Choose peak closest to psi_N = 1 if multiple peaks exist - peak_idx = peaks[numpy.argmax(psi_edge[peaks])] + # Among peaks in the far edge (psi_N > 0.85), pick the tallest; + # fall back to rightmost peak if none are in that zone. + far_edge = psi_edge[peaks] > 0.85 + if numpy.any(far_edge): + far_peaks = peaks[far_edge] + peak_idx = far_peaks[numpy.argmax(j_edge[far_peaks])] + else: + peak_idx = peaks[numpy.argmax(psi_edge[peaks])] peak_psi = psi_edge[peak_idx] peak_height = j_edge[peak_idx] @@ -374,21 +390,19 @@ def solve_jphi(mygs,ffp_prof,pp_prof,Ip_target,pax_target): mygs.solve() def find_optimal_scale(mygs, psi_N, pressure, ffp_prof, pp_prof, j_inductive, - Ip_target, psi_pad, spike_prof=None, find_j0=True, scale_j0=1.0, + Ip_target, psi_pad, spike_prof=None, tolerance=0.01, max_iter=5, diagnostic_plots=False, verbose=True): - r'''! Optimize scaling to match input/output \f$j_0\f$ or \f$I_p\f$ + r'''! Optimize core-\f$j_0\f$ scaling so input/output \f$j_0\f$ match @param mygs Grad-Shafranov solver object @param psi_N Normalized poloidal flux profile \f$\hat{\psi}\f$ @param pressure Plasma pressure profile \f$p(\hat{\psi})\f$ [Pa] @param ffp_prof Toroidal current profile (\f$j_\phi(\hat{\psi})\f$) - @param pp_prof Pressure gradient profile (\f{P'}\f$) + @param pp_prof Pressure gradient profile (\f{P\'}\f$) @param j_inductive Inductive/total toroidal current profile \f$j_{ind}(\hat{\psi})\f$ [A/m^2] @param Ip_target Target plasma current \f$I_p\f$ [A] @param psi_pad Padding for flux surface calculations @param spike_prof Optional bootstrap current profile \f$j_{BS}(\hat{\psi})\f$ [A/m^2] - @param find_j0 If True, rescale core \f$j_\phi\f$; else, rescale \f$I_p\f$ - @param scale_j0 Used if find_j0=False; rescaling factor for \f$j_0\f$ @param tolerance Relative error tolerance for secant method @param max_iter Maximum number of secant iterations @param diagnostic_plots If True, plot input/output \f$j_\phi\f$ @@ -397,15 +411,15 @@ def find_optimal_scale(mygs, psi_N, pressure, ffp_prof, pp_prof, j_inductive, import matplotlib.pyplot as plt n_psi = len(psi_N) - + if spike_prof is None: spike_prof = numpy.zeros_like(j_inductive) - + # We want: Input_J0 - Output_J0 ~ 0 def get_j0_error(scale_val, n): if verbose: print(f"\n--- Checking scale_j0 = {scale_val:.4f} ---") - + matched_input_jphi = scale_val*j_inductive + spike_prof ffp_prof['type'] = 'jphi-linterp' ffp_prof['y'] = matched_input_jphi @@ -413,7 +427,7 @@ def get_j0_error(scale_val, n): pax_target = pressure[0] solve_jphi(mygs,ffp_prof,pp_prof,Ip_target,pax_target) - + # Check Convergence _, f, fp, _, pp = mygs.get_profiles(npsi=n_psi, psi_pad=psi_pad) _, _, ravgs, _, _, _ = mygs.get_q(npsi=n_psi, psi_pad=psi_pad) @@ -430,61 +444,27 @@ def get_j0_error(scale_val, n): plt.ylabel(r'$j_\phi$ [MA/m$^2$]') plt.grid(ls=':') plt.show() - + # Input j_0 at this scale input_j0 = scale_val * j_inductive[0] + spike_prof[0] - + # Output j_0 from solver output_j0 = tmp_jphi[0] - + # Calculate residual (difference) and relative error diff = input_j0 - output_j0 rel_err = abs(diff) / output_j0 - + if verbose: print(f" Input j_0: {input_j0:.4e}") print(f" Output j_0: {output_j0:.4e}") print(f" Mismatch: {rel_err*100:.3f}%") - - return diff, rel_err, tmp_jphi - - # We want: Input_Ip - Output_Ip ~ 0 - def get_Ip_error(scale_val, scale_j0, n): - if verbose: - print(f"\n--- Checking scale_Ip = {scale_val:.4f} ---") - - ffp_prof['type'] = 'jphi-linterp' - ffp_prof['y'] = scale_j0*j_inductive + spike_prof - - scaled_Ip_target = Ip_target*scale_val - pax_target = pressure[0] - solve_jphi(mygs,ffp_prof,pp_prof,scaled_Ip_target,pax_target) - - eq_stats = mygs.get_stats(lcfs_pad=psi_pad) - output_Ip = eq_stats['Ip'] - - # Calculate residual (difference) and relative error - diff = output_Ip - Ip_target - # Use magnitude of Ip_target and guard against zero to avoid sign/zero issues - eps = numpy.finfo(float).eps - denom = max(abs(Ip_target), eps) - rel_err = abs(diff) / denom - - if verbose: - print(f" Input Ip target: {Ip_target/1e6:.4e}") - print(f" Trial Ip target: {Ip_target*scale_val/1e6:.4e}") - print(f" Output Ip: {output_Ip/1e6:.4e}") - print(f" Mismatch: {rel_err*100:.4f}%") - - return diff, rel_err, None + return diff, rel_err, tmp_jphi # --- Step 1: Initial Guess (1.0) --- p0 = 1.0 - if find_j0: - err0, rel_err0, res_jphi = get_j0_error(p0, 0) - else: - err0, rel_err0, res_jphi = get_Ip_error(p0, scale_j0, 0) + err0, rel_err0, res_jphi = get_j0_error(p0, 0) # Check if we got lucky immediately if rel_err0 < tolerance: @@ -493,22 +473,12 @@ def get_Ip_error(scale_val, scale_j0, n): # --- Step 2: Second Guess (Directional) --- # If Input < Output (err0 < 0), we need more current -> Try 1.2 # If Input > Output (err0 > 0), we have too much -> Try 0.8 - - if find_j0: - if err0 < 0: - p1 = 1.2 - else: - p1 = 0.8 + if err0 < 0: + p1 = 1.2 else: - if err0 < 0: - p1 = 1.01 - else: - p1 = 0.99 - - if find_j0: - err1, rel_err1, res_jphi = get_j0_error(p1, 1) - else: - err1, rel_err1, res_jphi = get_Ip_error(p1, scale_j0, 1) + p1 = 0.8 + + err1, rel_err1, res_jphi = get_j0_error(p1, 1) if rel_err1 < tolerance: return p1, res_jphi @@ -518,7 +488,7 @@ def get_Ip_error(scale_val, scale_j0, n): for i in range(max_iter): if verbose: print(f"--- Optimization Iteration {i+1} ---") - + # Avoid division by zero if abs(err1 - err0) < 1e-9: if verbose: @@ -528,25 +498,22 @@ def get_Ip_error(scale_val, scale_j0, n): # Secant formula: x_new = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0)) # This projects the line between (p0, err0) and (p1, err1) to zero. p_new = p1 - err1 * (p1 - p0) / (err1 - err0) - + # Safety clamp: If the projection is wild (negative or huge), constrain it # This acts as a loose bracket to prevent the solver from crashing - if p_new < 0.1: + if p_new < 0.1: p_new = 0.1 - if p_new > 5.0: + if p_new > 5.0: p_new = 5.0 - + # Execute Solver at new point - if find_j0: - err_new, rel_err_new, res_jphi = get_j0_error(p_new, i+2) - else: - err_new, rel_err_new, res_jphi = get_Ip_error(p_new, scale_j0, i+2) + err_new, rel_err_new, res_jphi = get_j0_error(p_new, i+2) if rel_err_new < tolerance: if verbose: print(f"Converged! Optimal scale factor: {p_new:.4f}") return p_new, res_jphi - + # Update points for next step (move window forward) p0, err0 = p1, err1 p1, err1 = p_new, err_new @@ -555,6 +522,7 @@ def get_Ip_error(scale_val, scale_j0, n): print("Max iterations reached. Returning best last effort.") return p1, res_jphi + def calculate_ln_lambda(Te, Ti, ne, ni, Zeff=1.0, electron_lnLambda_model='sauter', ion_lnLambda_model='unity'): @@ -1048,18 +1016,14 @@ def calculate_profiles_and_bootstrap(psi_N, include_jBS): # Find optimal jphi scale final_scale_j0, final_jphi = find_optimal_scale(mygs, psi_N, pressure, ffp_prof, pp_prof, matched_j_inductive, - Ip_target, psi_pad, spike_prof=spike_prof, find_j0=True, + Ip_target, psi_pad, spike_prof=spike_prof, diagnostic_plots=diagnostic_plots, verbose=verbose ) - if verbose: - print('\n >>> Finding optimal Ip scale factor') - # Find optimal Ip_target scale - final_scale_Ip, _ = find_optimal_scale(mygs, - psi_N, pressure, ffp_prof, pp_prof, matched_j_inductive, - Ip_target, psi_pad, spike_prof=spike_prof, find_j0=False, - scale_j0=final_scale_j0, - tolerance=0.001, diagnostic_plots=diagnostic_plots, verbose=verbose - ) + # Ip-scale secant removed: the solver now holds Ip to target + # natively (jphi-linterp cut-cell fix + Ip outer loop in the gs + # solve), so rescaling Ip_target here is redundant. The core-j0 + # secant above is retained. + final_scale_Ip = 1.0 if verbose: print('\n >>> Iterating on H-mode equilibrium solution') diff --git a/src/python/wrappers/tokamaker_f.F90 b/src/python/wrappers/tokamaker_f.F90 index 4edaeca52..f11239649 100644 --- a/src/python/wrappers/tokamaker_f.F90 +++ b/src/python/wrappers/tokamaker_f.F90 @@ -1623,7 +1623,7 @@ SUBROUTINE tokamaker_set_targets(tMaker_ptr,ip_target,ip_ratio_target,pax_target tMaker_obj%gs_equil%pax_target=pax_target*mu0 tMaker_obj%gs_equil%estore_target=estore_target*mu0 tMaker_obj%gs_equil%dflux_target=dflux_target -tMaker_obj%gs_equil%itor_target=ip_target*mu0 +tMaker_obj%gs_equil%ip_target=ip_target*mu0 tMaker_obj%gs_equil%ip_ratio_target=ip_ratio_target END SUBROUTINE tokamaker_set_targets !---------------------------------------------------------------------------------