Feature/dkes alpha

DKES/Sources/General/dkes_realspace.f

@@ -22,6 +22,7 @@ MODULE dkes_realspace
 
       ! ADDED BY SAL to ease optimizer use
       INTEGER :: DKES_rad_dex
+      INTEGER, ALLOCATABLE, DIMENSION(:) :: DKES_rundex
       REAL(rprec), ALLOCATABLE, DIMENSION(:) ::
      1   DKES_L11p, DKES_L33p, DKES_L31p,
      2   DKES_L11m, DKES_L33m, DKES_L31m,

STELLOPTV2/ObjectList

@@ -1,4 +1,5 @@
 ObjectFiles = \
+chisq_dkes_alpha.o \
 chisq_dkes_erdiff.o \
 chisq_line_visbrem.o \
 chisq_gamma_c.o \

STELLOPTV2/STELLOPTV2.dep

@@ -384,6 +384,13 @@ chisq_dkes_erdiff.o : \
       $(DKES_DIR)/$(LOCTYPE)/dkes_realspace.o
 
 
+chisq_dkes_alpha.o : \
+      stellopt_runtime.o \
+      stellopt_targets.o \
+      equil_vals.o \
+      $(DKES_DIR)/$(LOCTYPE)/dkes_realspace.o
+
+
 chisq_ne.o : \
       stellopt_runtime.o \
       stellopt_targets.o \

STELLOPTV2/Sources/Chisq/chisq_dkes.f90

@@ -50,7 +50,7 @@ SUBROUTINE chisq_dkes(target,sigma,niter,iflag)
       IF (iflag == 1) WRITE(iunit_out,'(A,2(2X,I3.3))') 'DKES ',ik,15
       IF (iflag == 1) WRITE(iunit_out,'(A)') 'TARGET  SIGMA  VAL  S  NU  ER  L11p  L11m  L33p  L33m  L31p  L31m  SCAL11  SCAL33  SCAL31'
       IF (niter >= 0) THEN
-         ik = 0
+         ik = 0 ! Always the first one
          DO ii = 1, nsd
             IF (sigma(ii) >= bigno) CYCLE
             DO ij = 1, nprof
@@ -79,7 +79,6 @@ SUBROUTINE chisq_dkes(target,sigma,niter,iflag)
             END DO
          END DO
       ELSE
-         nruns_dkes = 0
          DO ii = 1, nsd
             IF (sigma(ii) >= bigno) CYCLE
                DO ij = 1, nprof

STELLOPTV2/Sources/Chisq/chisq_dkes_alpha.f90

@@ -0,0 +1,108 @@
+!-----------------------------------------------------------------------
+!     Subroutine:    chisq_dkes_alpha
+!     Authors:       S. Lazerson ([email protected])
+!     Date:          12/19/2022
+!     Description:   This is the chisquared which handles the slope
+!                    of the DKES D*_11 coefficient.
+!-----------------------------------------------------------------------
+      SUBROUTINE chisq_dkes_alpha(target,sigma,niter,iflag)
+!-----------------------------------------------------------------------
+!     Libraries
+!-----------------------------------------------------------------------
+      USE stellopt_runtime
+      USE stellopt_targets
+      USE equil_vals, ONLY: shat
+!DEC$ IF DEFINED (DKES_OPT)
+      USE dkes_realspace, ONLY: DKES_L11p, DKES_L33p, DKES_L31p, &
+                                DKES_L11m, DKES_L33m, DKES_L31m, &
+                                DKES_scal11, DKES_scal33, DKES_scal31, &
+                                DKES_rundex
+!DEC$ ENDIF
+
+!-----------------------------------------------------------------------
+!     Input/Output Variables
+!
+!-----------------------------------------------------------------------
+      IMPLICIT NONE
+      REAL(rprec), INTENT(in)    ::  target(nsd)
+      REAL(rprec), INTENT(in)    ::  sigma(nsd)
+      INTEGER,     INTENT(in)    ::  niter
+      INTEGER,     INTENT(inout) ::  iflag
+
+!-----------------------------------------------------------------------
+!     Local Variables
+!
+!-----------------------------------------------------------------------
+      INTEGER :: ii, ij, ik
+      REAL(rprec) :: fp,fm
+!----------------------------------------------------------------------
+!     BEGIN SUBROUTINE
+!----------------------------------------------------------------------
+      IF (iflag < 0) RETURN
+      ! Count values
+      ik   = COUNT(sigma < bigno)
+      IF (iflag == 1) WRITE(iunit_out,'(A,2(2X,I3.3))') 'DKES_ALPHA ',ik,26
+      IF (iflag == 1) WRITE(iunit_out,'(A)') 'TARGET  SIGMA  VAL  S  NU+  NU- '&
+                // ' ER+  ER-' &
+                // ' L11p+  L11m+  L33p+  L33m+  L31p+  L31m+  SCAL11+  SCAL33+  SCAL31+'& 
+                // ' L11p-  L11m-  L33p-  L33m-  L31p-  L31m-  SCAL11-  SCAL33-  SCAL31-'
+      IF (niter >= 0) THEN
+         ik = COUNT(DKES_rundex < 3)
+         DO ii = 1, nsd
+            IF (sigma(ii) >= bigno) CYCLE
+            DO ij = 1, nprof
+               IF (Ep_DKES_alpha(ij) <= -bigno .or. nup_dkes_alpha(ij) <= -bigno .or. &
+                   Em_DKES_alpha(ij) <= -bigno .or. num_dkes_alpha(ij) <= -bigno) CYCLE
+               mtargets = mtargets + 1
+               targets(mtargets) = target(ii)
+               sigmas(mtargets)  = sigma(ii)
+!DEC$ IF DEFINED (DKES_OPT)
+               ik = ik + 1
+               fp = 0.5*(DKES_L11p(ik)+DKES_L11m(ik))
+               ik = ik + 1
+               fm = 0.5*(DKES_L11p(ik)+DKES_L11m(ik))
+               vals(mtargets) = (LOG10(fp) - LOG10(fm))/(LOG10(nup_dkes_alpha(ij))-LOG10(num_dkes_alpha(ij)))
+               IF (iflag == 1) WRITE(iunit_out,'(26ES22.12E3)') &
+                  targets(mtargets),sigmas(mtargets),vals(mtargets),&
+                  shat(ii), nup_dkes_alpha(ij), num_dkes_alpha(ij), &
+                  Ep_dkes_alpha(ij), Em_dkes_alpha(ij),&
+                  DKES_L11p(ik-1),DKES_L11m(ik-1),DKES_L33p(ik-1),&
+                  DKES_L33m(ik-1),DKES_L31p(ik-1),DKES_L31m(ik-1),&
+                  DKES_scal11(ik-1),DKES_scal33(ik-1),DKES_scal31(ik-1), &
+                  DKES_L11p(ik),DKES_L11m(ik),DKES_L33p(ik),&
+                  DKES_L33m(ik),DKES_L31p(ik),DKES_L31m(ik),&
+                  DKES_scal11(ik),DKES_scal33(ik),DKES_scal31(ik)
+
+!DEC$ ELSE
+               vals(mtargets) = target(ii)
+               IF (iflag == 1) WRITE(iunit_out,'(26ES22.12E3)') &
+                  targets(mtargets),sigmas(mtargets),vals(mtargets),&
+                  shat(ii), nup_dkes_Erdiff(ij), num_dkes_Erdiff(ij), &
+                  Ep_dkes_Erdiff(ij), Em_dkes_erdiff(ij),&
+                  0.0, 0.0, 0.0,&
+                  0.0, 0.0, 0.0,&
+                  0.0, 0.0, 0.0,&
+                  0.0, 0.0, 0.0,&
+                  0.0, 0.0, 0.0,&
+                  0.0, 0.0, 0.0
+!DEC$ ENDIF
+            END DO
+         END DO
+      ELSE
+         DO ii = 1, nsd
+            IF (sigma(ii) >= bigno) CYCLE
+            DO ij = 1, nprof
+               IF (Ep_DKES_alpha(ij) <= -bigno .or. nup_dkes_alpha(ij) <= -bigno .or. &
+                   Em_DKES_alpha(ij) <= -bigno .or. num_dkes_alpha(ij) <= -bigno) CYCLE
+               lbooz(ii) = .TRUE.
+               mtargets = mtargets + 1
+               nruns_dkes = nruns_dkes + 2 ! + and - for each run
+               IF (niter == -2) target_dex(mtargets)=jtarget_dkes_alpha
+            END DO
+         END DO
+      END IF
+      RETURN
+!----------------------------------------------------------------------
+!     END SUBROUTINE
+!----------------------------------------------------------------------
+      END SUBROUTINE chisq_dkes_alpha

STELLOPTV2/Sources/Chisq/chisq_dkes_erdiff.f90

@@ -15,7 +15,8 @@ SUBROUTINE chisq_dkes_erdiff(target,sigma,niter,iflag)
 !DEC$ IF DEFINED (DKES_OPT)
       USE dkes_realspace, ONLY: DKES_L11p, DKES_L33p, DKES_L31p, &
                                 DKES_L11m, DKES_L33m, DKES_L31m, &
-                                DKES_scal11, DKES_scal33, DKES_scal31
+                                DKES_scal11, DKES_scal33, DKES_scal31, &
+                                DKES_rundex
 !DEC$ ENDIF
 
 !-----------------------------------------------------------------------
@@ -45,7 +46,7 @@ SUBROUTINE chisq_dkes_erdiff(target,sigma,niter,iflag)
                 // ' L11p+  L11m+  L33p+  L33m+  L31p+  L31m+  SCAL11+  SCAL33+  SCAL31+'& 
                 // ' L11p-  L11m-  L33p-  L33m-  L31p-  L31m-  SCAL11-  SCAL33-  SCAL31-'
       IF (niter >= 0) THEN
-         ik = 0
+         ik = COUNT(DKES_rundex < 2)
          DO ii = 1, nsd
             IF (sigma(ii) >= bigno) CYCLE
             mtargets = mtargets + 1
@@ -83,16 +84,11 @@ SUBROUTINE chisq_dkes_erdiff(target,sigma,niter,iflag)
 !DEC$ ENDIF
          END DO
       ELSE
-         nruns_dkes = 0
          DO ii = 1, nsd
             IF (sigma(ii) >= bigno) CYCLE
             lbooz(ii) = .TRUE.
             mtargets = mtargets + 1
-            nruns_dkes = nruns_dkes + 2
-            E_dkes(1) = Ep_DKES_Erdiff
-            E_dkes(2) = Em_DKES_Erdiff
-            nu_dkes(1) = nu_dkes_Erdiff
-            nu_dkes(2) = nu_dkes_Erdiff
+            nruns_dkes = nruns_dkes + 2 ! Ep and Em run
             IF (niter == -2) target_dex(mtargets)=jtarget_dkes_erdiff
          END DO
       END IF

STELLOPTV2/Sources/General/stellopt_dkes.f90

@@ -13,7 +13,11 @@ SUBROUTINE stellopt_dkes(lscreen,iflag)
       USE equil_utils, ONLY: get_equil_phi, nrad, shat, phi_type
       USE stellopt_targets, ONLY: nu_dkes, sigma_dkes, lbooz, nsd, &
                                   E_dkes, nprof, nruns_dkes, &
-                                  sigma_dkes_erdiff
+                                  sigma_dkes_erdiff, Ep_DKES_Erdiff, &
+                                  Em_DKES_Erdiff, Ep_DKES_alpha, &
+                                  Em_DKES_alpha, sigma_dkes_alpha, &
+                                  nu_dkes_Erdiff, num_dkes_alpha, &
+                                  nup_dkes_alpha
       ! NEO LIBRARIES
 !DEC$ IF DEFINED (DKES_OPT)
       USE Vimatrix
@@ -66,6 +70,7 @@ SUBROUTINE stellopt_dkes(lscreen,iflag)
       CALL MPI_BCAST(nruns_dkes,1,MPI_INTEGER,master,MPI_COMM_MYWORLD,ierr_mpi)
 !DEC$ ENDIF
       ! Enter the main loop
+      IF (ALLOCATED(DKES_rundex)) DEALLOCATE(DKES_rundex)
       IF (ALLOCATED(DKES_L11p)) DEALLOCATE(DKES_L11p)
       IF (ALLOCATED(DKES_L33p)) DEALLOCATE(DKES_L33p)
       IF (ALLOCATED(DKES_L31p)) DEALLOCATE(DKES_L31p)
@@ -77,34 +82,83 @@ SUBROUTINE stellopt_dkes(lscreen,iflag)
       IF (ALLOCATED(DKES_scal31)) DEALLOCATE(DKES_scal31)
       ALLOCATE(DKES_L11p(nruns_dkes),DKES_L33p(nruns_dkes),DKES_L31p(nruns_dkes),&
                DKES_L11m(nruns_dkes),DKES_L33m(nruns_dkes),DKES_L31m(nruns_dkes),&
-               DKES_scal11(nruns_dkes),DKES_scal33(nruns_dkes),DKES_scal31(nruns_dkes))
+               DKES_scal11(nruns_dkes),DKES_scal33(nruns_dkes),DKES_scal31(nruns_dkes),&
+               DKES_rundex(nruns_dkes))
       DKES_L11p=0.0; DKES_L33p=0.0; DKES_L31p=0.0
       DKES_L11m=0.0; DKES_L33m=0.0; DKES_L31m=0.0
       DKES_scal11=0.0; DKES_scal33=0.0; DKES_scal31=0.0
+      DKES_rundex=-1;
       ! Setup the helper arrays
       IF (ALLOCATED(ik_dkes)) DEALLOCATE(ik_dkes)
       IF (ALLOCATED(nuarr_dkes)) DEALLOCATE(nuarr_dkes)
       IF (ALLOCATED(Earr_dkes)) DEALLOCATE(Earr_dkes)
       ALLOCATE(ik_dkes(nruns_dkes),nuarr_dkes(nruns_dkes),Earr_dkes(nruns_dkes))
       IF (myworkid == master) THEN
          ik = 0
+         ! First do traditional DKES
          DO ir = 1, nsd
-            IF (sigma_dkes(ir) >= bigno .and. sigma_dkes_erdiff(ir) >= bigno) CYCLE
+            IF (sigma_dkes(ir) >= bigno) CYCLE
             DO ij = 1, nprof
                IF (E_dkes(ij) <= -bigno .or. nu_dkes(ij) <= -bigno) CYCLE
                ik = ik + 1
-               ik_dkes(ik)    = ir
+               ik_dkes(ik) = ir
                nuarr_dkes(ik) = nu_dkes(ij)
-               Earr_dkes(ik)  = E_dkes(ij)
+               Earr_dkes(ik) = E_dkes(ij)
+               DKES_rundex(ik) = 1
             END DO
-         ENDDO
+         END DO
+         ! Now ErDiff
+         DO ir = 1, nsd
+            IF (sigma_dkes_erdiff(ir) >= bigno) CYCLE
+            ik = ik + 1
+            ik_dkes(ik) = ir
+            nuarr_dkes(ik) = nu_dkes_Erdiff
+            Earr_dkes(ik) = Ep_DKES_Erdiff
+            DKES_rundex(ik) = 2
+            ik = ik + 1
+            ik_dkes(ik) = ir
+            nuarr_dkes(ik) = nu_dkes_Erdiff
+            Earr_dkes(ik) = Em_DKES_Erdiff
+            DKES_rundex(ik) = 2
+         END DO
+         ! Now Alpha
+         DO ir = 1, nsd
+            IF (sigma_dkes_alpha(ir) >= bigno) CYCLE
+            DO ij = 1, nprof
+               IF (Ep_DKES_alpha(ij) <= -bigno .or. nup_dkes_alpha(ij) <= -bigno .or. &
+                   Em_DKES_alpha(ij) <= -bigno .or. num_dkes_alpha(ij) <= -bigno) CYCLE
+               ik = ik + 1
+               ik_dkes(ik) = ir
+               nuarr_dkes(ik) = nup_dkes_alpha(ij)
+               Earr_dkes(ik) = Ep_DKES_alpha(ij)
+               DKES_rundex(ik) = 3
+               ik = ik + 1
+               ik_dkes(ik) = ir
+               nuarr_dkes(ik) = num_dkes_alpha(ij)
+               Earr_dkes(ik) = Em_DKES_alpha(ij)
+               DKES_rundex(ik) = 3
+            END DO
+         END DO
+
+
+         !DO ir = 1, nsd
+         !   IF (sigma_dkes(ir) >= bigno .and. sigma_dkes_erdiff(ir) >= bigno) CYCLE
+         !   DO ij = 1, nprof
+         !      IF (E_dkes(ij) <= -bigno .or. nu_dkes(ij) <= -bigno) CYCLE
+         !      ik = ik + 1
+         !      ik_dkes(ik)    = ir
+         !      nuarr_dkes(ik) = nu_dkes(ij)
+         !      Earr_dkes(ik)  = E_dkes(ij)
+         !   END DO
+         !ENDDO
       END IF
       ! Now read the wout file
       CALL read_wout_file(proc_string, ier)
       ! Now bcast the boozer parameters
       CALL bcast_boozer_vars(master, MPI_COMM_MYWORLD, ierr_mpi)
 !DEC$ IF DEFINED (MPI_OPT)
       ierr_mpi = 0
+      CALL MPI_BCAST(DKES_rundex,nruns_dkes,MPI_INTEGER,master,MPI_COMM_MYWORLD,ierr_mpi)
       CALL MPI_BCAST(ik_dkes,nruns_dkes,MPI_INTEGER,master,MPI_COMM_MYWORLD,ierr_mpi)
       CALL MPI_BCAST(nuarr_dkes,nruns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi)
       CALL MPI_BCAST(Earr_dkes,nruns_dkes,MPI_DOUBLE_PRECISION,master,MPI_COMM_MYWORLD,ierr_mpi)

STELLOPTV2/Sources/General/stellopt_fcn.f90

@@ -435,7 +435,7 @@ SUBROUTINE stellopt_fcn(m, n, x, fvec, iflag, ncnt)
 !DEC$ ENDIF
 !DEC$ IF DEFINED (DKES_OPT)
          ctemp_str = 'dkes'
-         IF ((ANY(sigma_dkes < bigno).or.ANY(sigma_dkes_erdiff < bigno)) .and. (iflag>=0)) CALL stellopt_paraexe(ctemp_str,proc_string,lscreen); iflag = ier_paraexe
+         IF ((ANY(sigma_dkes < bigno).or.ANY(sigma_dkes_erdiff < bigno).or.ANY(sigma_dkes_alpha < bigno)) .and. (iflag>=0)) CALL stellopt_paraexe(ctemp_str,proc_string,lscreen); iflag = ier_paraexe
 !DEC$ ENDIF
 
          ! NOTE ALL parallel secondary codes go here

STELLOPTV2/Sources/General/stellopt_load_targets.f90

@@ -251,12 +251,15 @@ SUBROUTINE stellopt_load_targets(m, fvec, iflag, ncnt)
       ! NEO
       IF (ANY(sigma_neo < bigno)) &
          CALL chisq_neo(target_neo, sigma_neo, ncnt,iflag)
-      ! DKES
+      !!!!!!!! DKES Section only move as chunk !!!!!!!!!!!!!!!!!!!!!!!!
+      nruns_dkes = 0
       IF (ANY(sigma_dkes < bigno)) &
          CALL chisq_dkes(target_dkes, sigma_dkes, ncnt,iflag)
-      ! DKES ERdiff
       IF (ANY(sigma_dkes_erdiff < bigno)) &
          CALL chisq_dkes_erdiff(target_dkes_erdiff, sigma_dkes_erdiff, ncnt,iflag)
+      IF (ANY(sigma_dkes_alpha < bigno)) &
+         CALL chisq_dkes_alpha(target_dkes_alpha, sigma_dkes_alpha, ncnt,iflag)
+      !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       ! TXPORT
       IF (ANY(sigma_txport < bigno)) &
          CALL chisq_txport(target_txport, sigma_txport, ncnt,iflag)

STELLOPTV2/Sources/General/stellopt_main.f90

@@ -206,19 +206,27 @@ PROGRAM STELLOPT
       END IF
 
       IF (myworkid == master) THEN
+
          CALL stellopt_optimize
+
+!DEC$ IF DEFINED (MPI_OPT)
+         ! Only the master threads are part of MPI_COMM_STEL
+         ierr_mpi = 0
+         CALL MPI_COMM_FREE(MPI_COMM_STEL, ierr_mpi)
+         IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_FREE_ERR,'stellopt_main: MPI_COMM_STEL',ierr_mpi)
+!DEC$ ENDIF
+
+         ! Get workers
          ltst  = .false.
          tstr1 = 'exit'
          tstr2 = ''
          CALL stellopt_paraexe(tstr1,tstr2,ltst)
+
       END IF
 
       ! All procs (master and workers) will do this part
       ! Clean up
 !DEC$ IF DEFINED (MPI_OPT)
-      ierr_mpi = 0
-      CALL MPI_COMM_FREE(MPI_COMM_STEL, ierr_mpi)
-      IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_FREE_ERR,'stellopt_main',ierr_mpi)
       ierr_mpi = 0
       CALL MPI_FINALIZE(ierr_mpi)
       IF (ierr_mpi /= MPI_SUCCESS) CALL handle_err(MPI_FINE_ERR,'stellopt_main',ierr_mpi)