Resistive singular coupling

.github/workflows/make_test.yml

@@ -69,7 +69,7 @@ jobs:
       PREFIX: ${{ github.workspace }}/deps
       ONEAPI_ROOT: /opt/intel/oneapi
       FC: ifort
-      
+
     steps:
       - name: Checkout
         uses: actions/checkout@v4
@@ -102,7 +102,7 @@ jobs:
         if: steps.netcdf-cache.outputs.cache-hit != 'true'
         run: |
           source /opt/intel/oneapi/setvars.sh
-          
+
           mkdir -p $PREFIX/src
           cd $PREFIX/src
           curl -L -o netcdf-fortran.tar.gz https://github.com/Unidata/netcdf-fortran/archive/refs/tags/v4.6.1.tar.gz
@@ -128,12 +128,12 @@ jobs:
           else
             echo "Cache not found. Built from source."
           fi
-          
+
       - name: Confirm all shared libs
         run: |
           ls -lh $PREFIX/lib
           ls -lh $PREFIX/include
-          
+
   make_gpec_gfortran:
     runs-on: ubuntu-latest
     needs: build-netcdf-gfortran
@@ -159,7 +159,6 @@ jobs:
           git config --file .gitmodules submodule."deps/src/OpenBLAS".url https://github.com/OpenMathLib/OpenBLAS.git
           git submodule sync
           git submodule update --init --recursive "deps/src/harvest"
-
       - name: Install Dependencies
         run: |
           sudo apt update

README.md

@@ -20,7 +20,7 @@ On the PPPL portal computers, navigate to a directory in which you want to do yo
    git checkout -b develop
    git pull origin develop
 
-This will create a directory 'gpec' and with the developmental branch of the repository in it. For more on branches, see the workflow section below. For compile instructions, see the install section. Now you are free to add/edit files in this branch. Be sure to use:: 
+This will create a directory 'gpec' and with the developmental branch of the repository in it. For more on branches, see the workflow section below. For compile instructions, see the install section. Now you are free to add/edit files in this branch. Be sure to use::
 
    git add <new-file>
 
@@ -48,7 +48,7 @@ The highlights are,
   - There are two permanent branches: master and develop
   - The master branch is only updated for at release ready stages
   - New features should be developed in short-lived (days) branches coming off of and merging back to the development branch.
-  
+
 Specific instructions are given in the link above as to exactly how to branch and merge these various branches. For example, the --no-ff option should be used when merging in order to keep branch histories. Just follow the examples and you wont go wrong!
 
 
@@ -76,7 +76,7 @@ The very first time, starting from the root of the repo::
 
     cd docs
     make init
-    
+
 This creates a new gpec-docs directory in the same directory that contains the local repository (not in the repository itself) and checks out a devoted documentation branch. With that set up, the following will build the documentation and push it to the proper place::
 
     make ghpages

coil/coil.F

@@ -53,7 +53,7 @@ MODULE coil_mod
       REAL(r8), DIMENSION(:), POINTER :: shiftx, shifty, shiftz
       REAL(r8), DIMENSION(:), POINTER :: cur
       REAL(r8), DIMENSION(:,:,:), POINTER :: x,y,z
-      END TYPE coil_type      
+      END TYPE coil_type
 
       TYPE(coil_type), DIMENSION(:), POINTER :: coil
       TYPE(spline_type) :: csq
@@ -68,7 +68,7 @@ MODULE coil_mod
 c     declarations.
 c-----------------------------------------------------------------------
       SUBROUTINE coil_read(cdconfile,icoil_num,icoil_name,icoil_cur)
-      
+
       CHARACTER(128), INTENT(IN) :: cdconfile
       INTEGER, OPTIONAL :: icoil_num
       REAL(r8), DIMENSION(500,48), OPTIONAL :: icoil_cur
@@ -301,7 +301,7 @@ SUBROUTINE coil_read(cdconfile,icoil_num,icoil_name,icoil_cur)
                         r = sqrt( z ** 2 + y ** 2)
                         IF(cnpert == 1)THEN  ! want constant tilt angle like tilty case
                            angle = ATAN2(z, y) + tiltx
-                        ELSE 
+                        ELSE
                            ! makes peak tilt angle at y axis where r is maximum, for consistent max vertical displacement
                            ! however, this becomes redundent with tiltx for n=5 which is no good
     !                        angle = ATAN2(z, y) + tiltx*cos((cnpert-1)*(phi-pi/2))
@@ -315,7 +315,7 @@ SUBROUTINE coil_read(cdconfile,icoil_num,icoil_name,icoil_cur)
                      ! petal-like distortions of the shape
                      ! n=0 manually set to rigid shift
                      ! n=1 creates a beaning (one-petal, NOT rigid shift)
-                     ! n=2 creates a peanut shaping (two petal, NOT elliptic), etc. 
+                     ! n=2 creates a peanut shaping (two petal, NOT elliptic), etc.
                      IF(cnpert==0)THEN ! only shift, keep coil same
                         dx = dx + coil_shiftx(ci,cj)
                         dy = dy + coil_shifty(ci,cj)

coil/field.F

@@ -71,7 +71,7 @@ SUBROUTINE field_bs_psi(psi,bmn,wegt,op_start,op_stop, op_verbose)
      $     dx,dy,dz
 
       REAL(r8), DIMENSION(4) :: crzphi_f, crzphi_fx, crzphi_fy
-      
+
       ! splines are 0:cmtheta, but aspl_f_arr can be 1:cmtheta because
       ! we are repeating the first and last points to
       ! make the spline periodic
@@ -105,7 +105,7 @@ SUBROUTINE field_bs_psi(psi,bmn,wegt,op_start,op_stop, op_verbose)
 
 
       area=1.0
-      
+
 c-----------------------------------------------------------------------
 c     for each coil:
 c-----------------------------------------------------------------------
@@ -140,10 +140,10 @@ SUBROUTINE field_bs_psi(psi,bmn,wegt,op_start,op_stop, op_verbose)
                   xa(j,k,iseg)=(coil(i)%x(j,k,iseg+1)+
      $                 coil(i)%x(j,k,iseg))/2.0
                   ya(j,k,iseg)=(coil(i)%y(j,k,iseg+1)+
-     $                 coil(i)%y(j,k,iseg))/2.0      
+     $                 coil(i)%y(j,k,iseg))/2.0
                   za(j,k,iseg)=(coil(i)%z(j,k,iseg+1)+
      $                 coil(i)%z(j,k,iseg))/2.0
-                  
+
                   dx(j,k,iseg)=coil(i)%x(j,k,iseg+1)-coil(i)%x(j,k,iseg)
                   dy(j,k,iseg)=coil(i)%y(j,k,iseg+1)-coil(i)%y(j,k,iseg)
                   dz(j,k,iseg)=coil(i)%z(j,k,iseg+1)-coil(i)%z(j,k,iseg)
@@ -194,26 +194,26 @@ SUBROUTINE field_bs_psi(psi,bmn,wegt,op_start,op_stop, op_verbose)
                dbx=0
                dby=0
                dbz=0
-               
+
                DO j=1,coil(i)%ncoil
                   DO k=1,coil(i)%s
                      DO iseg=1,nseg
                         rx=xobs(itheta,izeta)-xa(j,k,iseg)
                         ry=yobs(itheta,izeta)-ya(j,k,iseg)
                         rz=zobs(itheta,izeta)-za(j,k,iseg)
-                        
+
                         dl=SQRT(rx**2+ry**2+rz**2)
                         dbx=dbx+coil(i)%cur(j)*coil(i)%nw*1e-7/
      $                       dl**3*(dy(j,k,iseg)*rz-dz(j,k,iseg)*ry)
                         dby=dby+coil(i)%cur(j)*coil(i)%nw*1e-7/
      $                       dl**3*(dz(j,k,iseg)*rx-dx(j,k,iseg)*rz)
                         dbz=dbz+coil(i)%cur(j)*coil(i)%nw*1e-7/
-     $                       dl**3*(dx(j,k,iseg)*ry-dy(j,k,iseg)*rx)               
-                        
+     $                       dl**3*(dx(j,k,iseg)*ry-dy(j,k,iseg)*rx)
+
                      ENDDO
                   ENDDO
                ENDDO
-               
+
                bx(itheta,izeta)=dbx
                by(itheta,izeta)=dby
                bz(itheta,izeta)=dbz
@@ -248,7 +248,7 @@ SUBROUTINE field_bs_psi(psi,bmn,wegt,op_start,op_stop, op_verbose)
          tbn=tbn+bn
       ENDDO
 
-         
+
       CALL cspline_alloc(bnzspl,cmzeta,1)
       bnzspl%xs=czeta
       CALL cspline_alloc(bntspl,cmtheta,1)
@@ -265,7 +265,7 @@ SUBROUTINE field_bs_psi(psi,bmn,wegt,op_start,op_stop, op_verbose)
          area=aspl%fsi(cmtheta,1)
          CALL spline_dealloc(aspl)
       ENDIF
-            
+
       DO ipert=1,cmpert
          DO izeta=1,cmzeta
             DO itheta=1,cmtheta
@@ -321,7 +321,7 @@ SUBROUTINE field_bs_rzphi(nr,nz,np,gdr,gdz,cbr,cbz,cbp,op_verbose)
       TYPE(cspline_type) :: bspl
 
       dlmin=2e-3*cro
-      
+
       ! Optional verbosity tells user how many threads we are using
       IF(PRESENT(op_verbose)) THEN
          verbose = op_verbose
@@ -344,7 +344,7 @@ SUBROUTINE field_bs_rzphi(nr,nz,np,gdr,gdz,cbr,cbz,cbp,op_verbose)
      $        dx(coil(i)%ncoil,coil(i)%s,nseg),
      $        dy(coil(i)%ncoil,coil(i)%s,nseg),
      $        dz(coil(i)%ncoil,coil(i)%s,nseg))
-         
+
          IF(verbose) WRITE(*,'(1x,a,i6,a,i4,a,i4,a,i4,a)') "  >"//
      $      "Mapping ", coil(i)%ncoil*coil(i)%s*nseg,
      $      "filaments from the coil "
@@ -359,10 +359,10 @@ SUBROUTINE field_bs_rzphi(nr,nz,np,gdr,gdz,cbr,cbz,cbp,op_verbose)
                   xa(j,k,iseg)=(coil(i)%x(j,k,iseg+1)+
      $                 coil(i)%x(j,k,iseg))/2.0
                   ya(j,k,iseg)=(coil(i)%y(j,k,iseg+1)+
-     $                 coil(i)%y(j,k,iseg))/2.0      
+     $                 coil(i)%y(j,k,iseg))/2.0
                   za(j,k,iseg)=(coil(i)%z(j,k,iseg+1)+
      $                 coil(i)%z(j,k,iseg))/2.0
-                  
+
                   dx(j,k,iseg)=coil(i)%x(j,k,iseg+1)-coil(i)%x(j,k,iseg)
                   dy(j,k,iseg)=coil(i)%y(j,k,iseg+1)-coil(i)%y(j,k,iseg)
                   dz(j,k,iseg)=coil(i)%z(j,k,iseg+1)-coil(i)%z(j,k,iseg)
@@ -387,7 +387,7 @@ SUBROUTINE field_bs_rzphi(nr,nz,np,gdr,gdz,cbr,cbz,cbp,op_verbose)
                   xobs(ir,iz,ip)=gdr(ir,iz)*COS(twopi*phi)
                   yobs(ir,iz,ip)=gdr(ir,iz)*SIN(twopi*phi)
                   zobs(ir,iz,ip)=gdz(ir,iz)
-                  
+
                   dbx=0
                   dby=0
                   dbz=0
@@ -398,20 +398,20 @@ SUBROUTINE field_bs_rzphi(nr,nz,np,gdr,gdz,cbr,cbz,cbp,op_verbose)
                            rx=xobs(ir,iz,ip)-xa(j,k,iseg)
                            ry=yobs(ir,iz,ip)-ya(j,k,iseg)
                            rz=zobs(ir,iz,ip)-za(j,k,iseg)
-                           
+
                            dl=MAXVAL((/SQRT(rx**2+ry**2+rz**2),dlmin/))
-                           
+
                            dbx=dbx+coil(i)%cur(j)*coil(i)%nw*1e-7/
      $                          dl**3*(dy(j,k,iseg)*rz-dz(j,k,iseg)*ry)
                            dby=dby+coil(i)%cur(j)*coil(i)%nw*1e-7/
      $                          dl**3*(dz(j,k,iseg)*rx-dx(j,k,iseg)*rz)
                            dbz=dbz+coil(i)%cur(j)*coil(i)%nw*1e-7/
-     $                          dl**3*(dx(j,k,iseg)*ry-dy(j,k,iseg)*rx)            
-                           
+     $                          dl**3*(dx(j,k,iseg)*ry-dy(j,k,iseg)*rx)
+
                         ENDDO
                      ENDDO
                   ENDDO
-               
+
                   bx(ir,iz,ip)=dbx
                   by(ir,iz,ip)=dby
                   bz(ir,iz,ip)=dbz
@@ -461,15 +461,15 @@ SUBROUTINE field_bs_rzphi(nr,nz,np,gdr,gdz,cbr,cbz,cbp,op_verbose)
       END SUBROUTINE field_bs_rzphi
 
       END MODULE field_mod
-
 
 
 
 
 
 
 
-
+
+
 
 
 

dcon/bal.f

@@ -511,7 +511,7 @@ SUBROUTINE bal_pse(theta,u)
 c-----------------------------------------------------------------------
 c     compute alpha series.
 c-----------------------------------------------------------------------
-      CALL spline_eval(bg,theta,0) 
+      CALL spline_eval(bg,theta,0)
       dtheta=theta-theta0
       v(1,1)=v0(1,1)+bg%f(1)/dtheta
       v(2,1)=v0(2,1)+bg%f(2)/dtheta

dcon/dcon.F

@@ -65,8 +65,7 @@ PROGRAM dcon
      $     bin_fmat,out_gmat,bin_gmat,out_kmat,bin_kmat,out_sol,
      $     out_sol_min,out_sol_max,bin_sol,bin_sol_min,bin_sol_max,
      $     out_fl,bin_fl,out_evals,bin_evals,bin_euler,euler_stride,
-     $     bin_vac,ahb_flag,mthsurf0,msol_ahb,netcdf_out,out_fund,
-     $     out_ahg2msc
+     $     bin_vac,ahb_flag,mthsurf0,msol_ahb,netcdf_out,out_fund
 c-----------------------------------------------------------------------
 c     format statements.
 c-----------------------------------------------------------------------
@@ -87,15 +86,6 @@ PROGRAM dcon
       READ(UNIT=in_unit,NML=dcon_output)
       CALL ascii_close(in_unit)
 
-      IF (out_ahg2msc) THEN
-         WRITE(*,*) "WARNING: ahg2msc.out is deprecated and will be " //
-     $          "removed in a future version. Set out_ahg2msc = .FALSE."
-         WRITE(*,*) "         to disable this warning."
-         vac_memory=.FALSE.
-      ELSE
-         vac_memory=.TRUE.
-      ENDIF
-
 c-----------------------------------------------------------------------
 c     set variables
 c-----------------------------------------------------------------------

dcon/dcon_mod.f

@@ -68,8 +68,6 @@ MODULE dcon_mod
       LOGICAL :: node_flag=.FALSE.
       LOGICAL :: res_flag=.FALSE.
       LOGICAL :: ahb_flag=.FALSE.
-      LOGICAL :: out_ahg2msc=.TRUE.
-      LOGICAL :: vac_memory=.FALSE.
 
       INTEGER, PARAMETER :: sol_base=50
       INTEGER :: mlow,mhigh,mpert,mband,nn,nstep=HUGE(0),bin_sol_min,

dcon/dcon_netcdf.f

@@ -154,9 +154,9 @@ SUBROUTINE dcon_netcdf_out(wp,wv,wt,wt0,ep,ev,et)
      $       "Normalized Poloidal Flux") )
       ENDIF
       IF(ALLOCATED(sing_detf))THEN
-         CALL check( nf90_def_dim(ncid, "psi_n_detf", 
+         CALL check( nf90_def_dim(ncid, "psi_n_detf",
      $       SIZE(sing_detf, 2), pd_dim) )
-         CALL check( nf90_def_var(ncid, "psi_n_detf", nf90_double, 
+         CALL check( nf90_def_var(ncid, "psi_n_detf", nf90_double,
      $       pd_dim, pd_id))
          CALL check( nf90_put_att(ncid,pd_id,"long_name",
      $       "Normalized Poloidal Flux") )

dcon/free.f

@@ -143,7 +143,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero,op_netcdf_out)
       ENDIF
       IF(wv_farwall_flag)THEN
          temp=wv
-      ENDIF         
+      ENDIF
 
       farwal_flag=.FALSE. ! self-inductance with the wall.
       kernelsignin=-1.0
@@ -168,7 +168,7 @@ SUBROUTINE free_run(plasma1,vacuum1,total1,nzero,op_netcdf_out)
       IF(wv_farwall_flag)THEN
          wv=temp
       ENDIF
-       
+
       singfac=mlow-nn*qlim+(/(ipert,ipert=0,mpert-1)/)
       DO ipert=1,mpert
          wv(ipert,:)=wv(ipert,:)*singfac
@@ -833,7 +833,7 @@ SUBROUTINE free_wvmats
             wv(:,ipert)=wv(:,ipert)*singfac
          ENDDO
          wvmats%fs(i,:)=RESHAPE(wv,(/mpert**2/))
-         DEALLOCATE(grri,xzpts)         
+         DEALLOCATE(grri,xzpts)
       ENDDO
       CALL unset_dcon_params
       CALL cspline_fit(wvmats,"extrap")