Merge pull request #339 from PrincetonUniversity/feature/THRIFT_print…

…_DKES_coeffs

Feature/thrift print dkes coeffs

LIBSTELL/Sources/Modules/thrift_globals.f90

@@ -41,6 +41,7 @@ MODULE thrift_globals
       INTEGER, DIMENSION(:), POINTER :: DKES_rundex
       INTEGER, DIMENSION(DKES_NS_MAX) :: DKES_K
       REAL(rprec), DIMENSION(DKES_NSTAR_MAX) :: DKES_Erstar, DKES_Nustar
+      LOGICAL :: save_DKES_coeffs
 
       ! Moved from thrift_runtime
       INTEGER :: nparallel_runs, mboz, nboz

LIBSTELL/Sources/Modules/thrift_input_mod.f90

@@ -36,7 +36,7 @@ MODULE thrift_input_mod
                               freq_ecrh, power_ecrh, &
                               pecrh_aux_t, pecrh_aux_f, ecrh_rc, ecrh_w, &
                               dkes_k, dkes_Erstar, dkes_Nustar, &
-                              etapar_type
+                              etapar_type, save_DKES_coeffs
 
 !-----------------------------------------------------------------------
 !     Subroutines
@@ -90,6 +90,7 @@ SUBROUTINE init_thrift_input
       dkes_k = -1
       dkes_Erstar = 1E10
       dkes_Nustar = 1E10
+      save_DKES_coeffs = .FALSE.
       RETURN
       END SUBROUTINE init_thrift_input
 

PENTA/Sources/penta.f90

@@ -769,8 +769,8 @@ Program penta3
   Else
     Er_test_vals(min_ind) = Er_test_vals(min_ind + 1)/2._rknd
   EndIf
-  Write(*,'(a,i4,a,f10.3)') 'Cannot use Er=0 with log_interp, using Er(',  &
-     min_ind, ') = ', Er_test_vals(min_ind)
+  ! Write(*,'(a,i4,a,f10.3)') 'Cannot use Er=0 with log_interp, using Er(',  &
+  !    min_ind, ') = ', Er_test_vals(min_ind)
 EndIf
 
 ! Loop over Er to get fluxes as a function of Er

PENTA/Sources/read_input_file_mod.f90

@@ -260,19 +260,19 @@ Subroutine read_vmec_file_2(js,run_ident)
     b0 = dsqrt(bsq)
 
     !Print
-    write(*,*) '##################################'
-    write(*,*) '##### GEOMETRICAL QUANTITIES: ####'
-    write(*,*) '##################################'
-    write(*,'(a,g0)') ' Bsq=',bsq
-    write(*,'(a,g0)') ' NOTE!!! USING B0=sqrt(Bsq)=',b0
-    write(*,'(a,g0)') ' psip/2pi=',psip
-    write(*,'(a,g0)') ' chip/2pi=',chip
-    write(*,'(a,g0)') ' psip/2pi=',psip
-    write(*,'(a,g0)') ' dVdr=',vol_p
-    write(*,'(a,g0)') ' btheta=',btheta
-    write(*,'(a,g0)') ' bzeta=',bzeta
-    write(*,*) '##################################'
-    write(*,*) '##################################'
+    ! write(*,*) '##################################'
+    ! write(*,*) '##### GEOMETRICAL QUANTITIES: ####'
+    ! write(*,*) '##################################'
+    ! write(*,'(a,g0)') ' Bsq=',bsq
+    ! write(*,'(a,g0)') ' NOTE!!! USING B0=sqrt(Bsq)=',b0
+    ! write(*,'(a,g0)') ' psip/2pi=',psip
+    ! write(*,'(a,g0)') ' chip/2pi=',chip
+    ! write(*,'(a,g0)') ' psip/2pi=',psip
+    ! write(*,'(a,g0)') ' dVdr=',vol_p
+    ! write(*,'(a,g0)') ' btheta=',btheta
+    ! write(*,'(a,g0)') ' bzeta=',bzeta
+    ! write(*,*) '##################################'
+    ! write(*,*) '##################################'
 
   Endsubroutine read_vmec_file_2
 

THRIFT/Sources/thrift_dkes.f90

@@ -41,14 +41,15 @@ SUBROUTINE thrift_dkes(lscreen,iflag)
 !        ier         Error flag
 !        iunit       File unit number
 !----------------------------------------------------------------------
-      INTEGER :: i, j, k, l, istat, neqs, ier_phi, mystart, myend
+      INTEGER :: i, j, k, l, istat, neqs, ier_phi, mystart, myend, iu_dkes_coeffs
       INTEGER, DIMENSION(:), ALLOCATABLE :: ik_dkes
       REAL(rprec), DIMENSION(:), ALLOCATABLE :: Earr_dkes, nuarr_dkes
       REAL(rprec), DIMENSION(:), POINTER :: f0p1, f0p2, f0m1, f0m2
       REAL(rprec) :: tcpu0, tcpu1, tcpui, tcput, tcpu, tcpua, &
-                     phi_temp, stime, etime
+                     phi_temp, stime, etime, &
+                     dkes_d11_save, dkes_d31_save, dkes_d33_save
       CHARACTER :: tb*1           
-      CHARACTER*50 :: arg1(6)       
+      CHARACTER*50 :: arg1(6)  
       INTEGER :: numargs, iodata, iout_opt, idata, iout
       INTEGER :: m, n ! nmax, mmax (revmoed due to conflict with dkes_realspace
       CHARACTER :: output_file*64, opt_file*64, dkes_input_file*64, temp_str*64
@@ -350,6 +351,29 @@ SUBROUTINE thrift_dkes(lscreen,iflag)
             DKES_D11 = RESHAPE( (DKES_L11p + DKES_L11m)*0.5, shape=(/DKES_NK, DKES_NC, DKES_NE/), order=(/2,3,1/) )
             DKES_D31 = RESHAPE( (DKES_L31p + DKES_L31m)*0.5, shape=(/DKES_NK, DKES_NC, DKES_NE/), order=(/2,3,1/) )
             DKES_D33 = RESHAPE( (DKES_L33p + DKES_L33m)*0.5, shape=(/DKES_NK, DKES_NC, DKES_NE/), order=(/2,3,1/) )
+
+            IF(save_DKES_coeffs) THEN
+               ! Open file
+               WRITE(temp_str,'(i3.3)') mytimestep
+               CALL safe_open(iu_dkes_coeffs, istat, "DKES_coeffs."//TRIM(temp_str), &
+               Trim(Adjustl('unknown')), 'formatted')
+
+               ! Write header
+               WRITE(iu_dkes_coeffs, '(A)') "dkes_k  Er_v           nu_v           D11            D31            D33"
+
+               ! Write data row by row
+               DO l=1,nruns_dkes
+                  dkes_d11_save = (DKES_L11p(l) + DKES_L11m(l))*0.5
+                  dkes_d31_save = (DKES_L31p(l) + DKES_L31m(l))*0.5
+                  dkes_d33_save = (DKES_L33p(l) + DKES_L33m(l))*0.5
+                  WRITE(iu_dkes_coeffs, '(I6, 5E15.7)') ik_dkes(l), Earr_dkes(l), nuarr_dkes(l), dkes_d11_save, dkes_d31_save, dkes_d33_save
+               END DO
+
+               ! Close file
+               CLOSE(iu_dkes_coeffs, iostat=ier)
+               IF (ier /= 0) STOP 'Error closing dkes_coeffs file'              
+            ENDIF
+
          END IF
       END IF
       IF (lscreen) WRITE(6,'(a)') ' -------------------  DKES NEOCLASSICAL CALCULATION DONE  ---------------------'

pySTEL/libstell/penta.py

@@ -12,7 +12,7 @@
 # PENTA Class
 class PENTA:
 
-    def __init__(self, folder_path, plasma=None, Zions=None):
+    def __init__(self, folder_path, plasma=None, Zions=None, lverb=True):
         #folder_path is a path to the folder containnig the following PENTA3 results files:
         # - flows_vs_Er
         # - flows_vs_roa
@@ -22,7 +22,7 @@ def __init__(self, folder_path, plasma=None, Zions=None):
         # - contra_vs_roa
         # - sigmas_vs_roa
 
-        print('\nPENTA class being created...')
+        if(lverb): print('\nPENTA class being created...')
 
         self.folder_path = folder_path
 
@@ -33,25 +33,25 @@ def __init__(self, folder_path, plasma=None, Zions=None):
         elif(plasma is not None and Zions is not None):
             print('Both plasma class and Zions provided. Considering plasma class and discarding Zions array')
             self.list_of_species = plasma.list_of_species
-            print(f'List of species: {self.list_of_species}')
+            if(lverb): print(f'List of species: {self.list_of_species}')
             self.Zions = [plasma.Zcharge[species] for species in self.list_of_species]
             #remove electrons
             self.Zions = self.Zions[1:]
             print(f'Zions={self.Zions}')
         elif(plasma is not None and Zions is None):
             self.list_of_species = plasma.list_of_species
-            print(f'List of species: {self.list_of_species}')
+            if(lverb): print(f'List of species: {self.list_of_species}')
             self.Zions = [plasma.Zcharge[species] for species in self.list_of_species]
             #remove electrons
             self.Zions = self.Zions[1:]
-            print(f'Zions={self.Zions}')
+            if(lverb): print(f'Zions={self.Zions}')
         elif(plasma is None and Zions is not None):
             #check if size of Zions is compatible with number of ions in results files
             self.check_size_Zions(Zions)
             self.list_of_species = ['e'] + [f'i{k}' for k in range(1,len(Zions)+1)]
-            print(f'List of species: {self.list_of_species}')
+            if(lverb): print(f'List of species: {self.list_of_species}')
             self.Zions = Zions
-            print(f'Zions={self.Zions}')
+            if(lverb): print(f'Zions={self.Zions}')
 
         #sets the arrays self.roa_unique and self.Er_search
         self.set_independent_variables()
@@ -64,9 +64,14 @@ def __init__(self, folder_path, plasma=None, Zions=None):
         # - self.uprl[species,root]
         # - self.Jprl[species,root]
         # - self.Gamma[species,root]
+        # - self.QoT[species,root]
         # species is one of the species in self.list_of_species and root should be 
         # 'ion_root', 'electron_root' or 'unstable_root'
         self.set_fluxes_flows_by_root()
+
+        # Sets the dictionary self.root_Maxwell[roa] = 'ion_roots' or 'electron_root'
+        # using Maxwell construction criterium
+        self.set_Maxwell_root()
 
 
     def check_size_Zions(self,Zions):
@@ -90,19 +95,18 @@ def set_independent_variables(self):
 
         penta = np.loadtxt(filename,skiprows=2) 
 
-        self.roa_unique = np.unique(penta[:,0])
-        self.Er_search = np.unique(penta[1,:])
-
+        self.Er_search = np.unique(penta[1,:])          
+
         filename = self.folder_path + '/flows_vs_roa'
 
         penta = np.loadtxt(filename,skiprows=2)
 
+        self.roa_unique = np.unique(penta[:,0])
+
         number_flows_per_species = len(penta[0,3:]) / len(self.list_of_species)
 
         self.Smax = int(number_flows_per_species - 1)
 
-        print(f'Smax={self.Smax}')
-
     def set_variables_by_root(self):
         # sets the dictionaries self.##[root], self.##[root]], self.##[root]],
         # with ## being: roa, Er, Jprl_total, J_BS
@@ -123,10 +127,12 @@ def set_variables_by_root(self):
         self.Er = defaultdict(list)
         self.Jprl_total = defaultdict(list)
         self.JBS = defaultdict(list)
+        self.num_roots = []
 
         i=0
         for _, group in groupby(roa):
             num_roots = len( list(group) )
+            self.num_roots.append(num_roots)
 
             if(num_roots == 1):
                 self.roa['ion_root'].append(roa[i])
@@ -159,10 +165,10 @@ def set_variables_by_root(self):
                 print(f"How come you have {num_roots} roots ??")
                 exit(0)   
 
-            i += num_roots        
+            i += num_roots    
 
     def set_fluxes_flows_by_root(self):
-        # sets the dictionaries self.uprl[species,root], self.Jprl[species,root] and self.Gamma[species,root]
+        # sets the dictionaries self.uprl[species,root], self.Jprl[species,root], self.Gamma[species,root] and self.QoT[species,root]
         # they dictionaries return arrays
         # species is one of the species in self.list_of_species and root should be 'ion_root', 'electron_root' or 'unstable_root'
 
@@ -180,14 +186,16 @@ def set_fluxes_flows_by_root(self):
         penta = np.loadtxt(filename,skiprows=2)
         Jprl_penta = penta[:,3:(3+len(self.list_of_species))]
 
-        #get particle fluxes for all species
+        #get fluxes for all species
         filename = self.folder_path + '/fluxes_vs_roa'
         penta = np.loadtxt(filename,skiprows=2)
         Gamma_penta = penta[:,np.r_[3,5:(5+len(self.Zions))]]
+        QoT_penta = penta[:,np.r_[4,(5+len(self.Zions)):]]
 
         self.uprl = defaultdict(list)
         self.Jprl = defaultdict(list)
         self.Gamma = defaultdict(list)
+        self.QoT = defaultdict(list)
 
         for k,species in enumerate(self.list_of_species):
 
@@ -199,6 +207,7 @@ def set_fluxes_flows_by_root(self):
                     self.uprl[species,'ion_root'].append(uprl0_penta[i,k])
                     self.Jprl[species,'ion_root'].append(Jprl_penta[i,k])
                     self.Gamma[species,'ion_root'].append(Gamma_penta[i,k])
+                    self.QoT[species,'ion_root'].append(QoT_penta[i,k])
                 elif(num_roots ==3 or num_roots>3):
 
                     if(num_roots>3):
@@ -208,21 +217,86 @@ def set_fluxes_flows_by_root(self):
                     self.uprl[species,'ion_root'].append(uprl0_penta[i,k])
                     self.Jprl[species,'ion_root'].append(Jprl_penta[i,k])
                     self.Gamma[species,'ion_root'].append(Gamma_penta[i,k])
+                    self.QoT[species,'ion_root'].append(QoT_penta[i,k])
                     # unstable root
                     self.uprl[species,'unstable_root'].append(uprl0_penta[i+1,k])
                     self.Jprl[species,'unstable_root'].append(Jprl_penta[i+1,k])
                     self.Gamma[species,'unstable_root'].append(Gamma_penta[i+1,k])
+                    self.QoT[species,'unstable_root'].append(QoT_penta[i+1,k])
                     # electron root
                     self.uprl[species,'electron_root'].append(uprl0_penta[i+2,k])
                     self.Jprl[species,'electron_root'].append(Jprl_penta[i+2,k])
                     self.Gamma[species,'electron_root'].append(Gamma_penta[i+2,k])
+                    self.QoT[species,'electron_root'].append(QoT_penta[i+2,k])
                 elif(num_roots ==2):
                     raise ValueError(f"How come you have 2 roots ??") 
                 else:
                     print(f"How come you have {num_roots} roots ??")
                     exit(0)    
                 i += num_roots
 
+    def set_Maxwell_root(self):
+        # sets the dictionary self.root_Maxwell
+        # Of all the ambipolar roots, defines which one will settle using Maxwell criterium
+
+        from collections import defaultdict
+
+        filename = self.folder_path + '/fluxes_vs_Er'
+
+        penta = np.loadtxt(filename,skiprows=2)
+
+        roa = penta[:,0]
+        Er = penta[:,1]
+        gamma_e = penta[:,2]       
+        gamma_i_tot = np.sum(penta[:,3:]*self.Zions,axis=1)
+
+        Jr = gamma_i_tot - gamma_e
+
+        self.root_Maxwell = {}
+        self.Er_Maxw = []
+        self.JBS_Maxw = []
+        self.Gamma_Maxw = defaultdict(list)
+        self.QoT_Maxw = defaultdict(list)
+
+        for k,rho in enumerate(self.roa_unique):
+            num_roots = self.num_roots[k]
+
+            if(num_roots==1):
+                root_type = 'ion_root'
+                self.root_Maxwell[rho] = 'ion_root'
+                idx_Maxw = np.abs(self.roa['ion_root']-rho).argmin()
+            elif(num_roots>=3):
+                # compute integral(Jr.dEr) between electron and ion roots
+                ie = np.abs(self.roa['electron_root']-rho).argmin()
+                electron_root = self.Er['electron_root'][ie]
+                #
+                ii = np.abs(self.roa['ion_root']-rho).argmin()
+                ion_root = self.Er['ion_root'][ii]
+                # find Er closest to electron and ion roots
+                idx_e = np.abs(Er - electron_root).argmin()
+                idx_i = np.abs(Er - ion_root).argmin()
+                # Slice the arrays for the specified range
+                Er_slice = Er[idx_e:idx_i+1]
+                Jr_slice = Jr[idx_e:idx_i+1]
+                # Perform the integration using np.trapz
+                integral = np.trapz(Jr_slice, Er_slice)
+                # Decide root
+                if(integral>0):
+                    self.root_Maxwell[rho] = 'ion_root'
+                    root_type = 'ion_root'
+                    idx_Maxw = ii
+                else:
+                    self.root_Maxwell[rho] = 'electron_root'
+                    root_type = 'electron_root'
+                    idx_Maxw = ie
+
+            # Now save in Maxwell arrays
+            self.Er_Maxw.append(self.Er[root_type][idx_Maxw])
+            self.JBS_Maxw.append(self.JBS[root_type][idx_Maxw])
+            for species in self.list_of_species:
+                self.Gamma_Maxw[species].append(self.Gamma[species,root_type][idx_Maxw])                  
+                self.QoT_Maxw[species].append(self.QoT[species,root_type][idx_Maxw])
+
     def plot_Er_vs_roa(self,which_root='all',plot=True):
         # plots ambipolar Er vs roa
         # which_root is: 'ion_root', 'unstable_root', 'electron_root' or 'all'
@@ -418,7 +492,9 @@ def plot_fluxes_vs_Er(self,roa_user,plot=True):
 
         roa = penta[:,0]
         Er = penta[:,1]
-        gamma_e = penta[:,2]       
+        gamma_e = penta[:,2]  
+
+        roa_unique = np.unique(roa)     
 
         gamma_i_tot = np.sum(penta[:,3:]*self.Zions,axis=1)
 
@@ -442,7 +518,7 @@ def plot_fluxes_vs_Er(self,roa_user,plot=True):
                 print(f'ERROR!! The provided roa={r_user} is outside the interval of PENTA data: [{np.min(roa)},{np.max(roa)}]')
                 exit(0)        
             else:
-                roa_closest = self.roa_unique[ np.argmin(np.abs(self.roa_unique-r_user)) ]
+                roa_closest = roa_unique[ np.argmin(np.abs(roa_unique-r_user)) ]
 
             plt.rc('font', size=16)
             fig=plt.figure(figsize=(8,6))

pySTEL/libstell/plasma.py

@@ -260,7 +260,7 @@ def get_temperature(self,species,rho):
 
         # check if temperature of species has been set
         if(species not in self.temperature):
-            print('ERROR" temperatureof {species} has not been set yet')
+            print('ERROR" temperature of {species} has not been set yet')
             exit(0)
 
         # make sure rho is an array
@@ -460,7 +460,7 @@ def write_plasma_profiles_to_PENTA1(self,rho,filename=None):
                 # Write the row to the file, formatted as space-separated values
                 file.write(" ".join(map(str, row_data)) + '\n')
 
-    def write_plasma_profiles_to_PENTA3(self,rho,filename=None):
+    def write_plasma_profiles_to_PENTA3(self,rho=np.linspace(0,1,200),filename=None):
         # first line: number of rhos
         # from second line: 
         # 1st column: rho=r/a