power balance in gyro HF loop

source/CADClass.py

@@ -380,7 +380,7 @@ def writeMesh2file(self, mesh, label, path='./', resolution=None):
         """
         Writes a mesh object to STL file named by part number.
         If mesh is a list of mesh objects, then write a separate file for
-        each mesh object in the list.  Does not overwrite / clobber
+        each mesh object in the list.  Clobbers if overWriteMask is True
         """
         if resolution == None:
             resolution = self.ROIGridRes
@@ -406,7 +406,7 @@ def writeMesh2file(self, mesh, label, path='./', resolution=None):
             filename = path + label[i] + "___" + resolution +"mm.stl"
             print("Writing mesh file: " + filename)
             log.info("Writing mesh file: " + filename)
-            if os.path.exists(filename):
+            if os.path.exists(filename) and self.overWriteMask == False:
                 pass
             else:
                 mesh[i].write(filename)

source/gyroClass.py

@@ -178,7 +178,7 @@ def setupVelocities(self, N):
             #set upper bound of v*f(v) (note that this cuts off high energy particles)
             self.vMax = 5 * self.vThermal
             #get 100 points to initialize functional form of f(v) (note this is a 2D matrix cause vMax is 2D)
-            self.vScan = np.linspace(0,self.vMax,100).T
+            self.vScan = np.linspace(0,self.vMax,1000).T
             #get velocity slices for each T0
             self.pullEqualProbabilityVelocities()
 
@@ -198,6 +198,8 @@ def pullEqualProbabilityVelocities(self):
         """
         self.vSlices = np.ones((len(self.T0),self.N_vSlice))*np.nan
         self.energySlices = np.zeros((len(self.T0),self.N_vSlice))
+        self.energyIntegrals = np.zeros((len(self.T0),self.N_vSlice))
+        self.energyFracs = np.zeros((len(self.T0),self.N_vSlice))
         for i in range(len(self.T0)):
             #get velocity range for this T0
             v = self.vScan[i,:]
@@ -209,38 +211,30 @@ def pullEqualProbabilityVelocities(self):
             v_cdf = np.cumsum(v_pdf[1:])*np.diff(v)
             v_cdf = np.insert(v_cdf, 0, 0)
             #create bspline interpolators for the cdf and cdf inverse
-            inverseCDF = interp1d(v_cdf, v, kind='cubic')
-            forwardCDF = interp1d(v, v_cdf, kind='cubic')
+            inverseCDF = interp1d(v_cdf, v, kind='linear')
+            forwardCDF = interp1d(v, v_cdf, kind='linear')
             #calculate N_vSlice velocities for each pdf each with equal area (probability)
             cdfMax = v_cdf[-1]
             cdfMin = v_cdf[0]
             sliceWidth = cdfMax / (self.N_vSlice+1)
             cdfSlices = np.linspace(0,1,self.N_vSlice+2)[1:-1]
+            cdfBounds = np.linspace(0,1,self.N_vSlice+1)[1:-1]
             self.vSlices[i,:] = inverseCDF(cdfSlices)
-
+            vBounds = inverseCDF(cdfBounds)
+            vBounds = np.insert(vBounds,0,0)
+            vBounds = np.append(vBounds,self.vMax[i])
             #Now find energies that correspond to these vSlices
             #we integrate: v**2 * f(v)
             #energy pdf (missing 1/2*mass but that gets divided out later anyways )
-            energy = lambda x: x**2 * pdf(x)
-            #if there is only 1 vSlice integrate entire pdf
-            if len(self.vSlices[i]==1):
-                vLo = 0.0
-                vHi = self.vMax[i]
-                self.energySlices[i] = integrate.quad(energy, vLo, vHi)[0]
-            #if there are multiple vSlices use them as integral bounds
-            else:
-                for j in range(len(self.vSlices[i])-1):
-                    if j==0:
-                        vLo = 0.0
-                        vHi = self.vMax[i,0]
-                    elif j==len(self.vSlices[i])-2:
-                        vLo = self.vSlices[i,-1]
-                        vHi = self.vMax[i]
-                    else:
-                        vLo = self.vSlices[i,j-1]
-                        vHi = self.vSlices[i,j]
-
-                    self.energySlices[i,j] = integrate.quad(energy, vLo, vHi)[0]
+            EofV = lambda x: x**2 * pdf(x)
+            #energy slices that correspond to velocity slices
+            self.energySlices[i,:] = EofV(self.vSlices[i,:])
+            #energy integrals and fractions
+            for j in range(self.N_vSlice):
+                self.energyIntegrals[i,j] = integrate.quad(EofV, vBounds[j], vBounds[j+1])[0]
+            energyTotal = self.energyIntegrals[i,:].sum()
+            for j in range(self.N_vSlice):
+                self.energyFracs[i,j] = self.energyIntegrals[i,j] / energyTotal
 
         print("Found N_vPhase velocities of equal probability")
         log.info("Found N_vPhase velocities of equal probability")

source/heatfluxClass.py

@@ -976,53 +976,50 @@ def gyroHF(self, GYRO, PFC):
         """
         print("Calculating gyro orbit heat loads")
         log.info("Calculating gyro orbit heat loads")
-        #get q|| for this PFC surface
-        q = PFC.q[PFC.PFC_GYROmap]
         #get divertor HF
         qDiv = PFC.qDiv[PFC.PFC_GYROmap] / self.elecFrac
         Pdiv = qDiv * PFC.areas[PFC.PFC_GYROmap]
         #Get fractional multipliers for each helical trace
         gyroFrac = 1.0/GYRO.N_gyroPhase
         vPhaseFrac = 1.0/GYRO.N_vPhase
-        #energy: note these units are bad, but get divided out
-        #old method (only converges for high resolution cases)
-        #energies = 1.0/2.0 * GYRO.mass_eV * GYRO.vSlices**2
-        energyTotal = GYRO.energySlices.sum(axis=1)
-        vSliceFrac = np.zeros((len(q),GYRO.N_vSlice))
-
-        for vSlice in range(GYRO.N_vSlice):
-            vSliceFrac[:,vSlice] = GYRO.energySlices[:,vSlice] / energyTotal
-
+#        #energy: note these units are bad, but get divided out
+#        #old method (only converges for high resolution cases)
+#        #energies = 1.0/2.0 * GYRO.mass_eV * GYRO.vSlices**2
+#        energyTotal = GYRO.energySlices.sum(axis=1)
+#        vSliceFrac = np.zeros((len(q),GYRO.N_vSlice))
+#
+#        for vSlice in range(GYRO.N_vSlice):
+#            vSliceFrac[:,vSlice] = GYRO.energySlices[:,vSlice] / energyTotal
+        vSliceFrac = GYRO.energyFracs
         #qMatrix = np.zeros((GYRO.N_gyroPhase,GYRO.N_vPhase,GYRO.N_vSlice,len(q)))
         Pgyro = np.zeros((GYRO.Nt))
         PNaN = 0.0
+        sum=0
+        sum1=0
         #loop through intersect record and redistribute power using multipliers
         for gyroPhase in range(GYRO.N_gyroPhase):
             for vPhase in range(GYRO.N_vPhase):
                 for vSlice in range(GYRO.N_vSlice):
                     idx = GYRO.intersectRecord[gyroPhase,vPhase,vSlice,PFC.CADHOT_GYROmap]
-                    isNan = np.where(np.isnan(idx)==True)[0] #include NaNs (NaNs = no intersection)
-                    notNan = np.where(np.isnan(idx)==False)[0] #dont include NaNs (NaNs = no intersection)
-                    idx1 = idx[~np.isnan(idx)] #indices we map power to
-                    idx1 = idx1.astype(int) #cast as integer
-                    idx2 = idx[np.isnan(idx)] #indices we map power to
-                    idx2 = idx2.astype(int) #cast as integer
-                    #qGyro[idx1] += q[notNan]*GYRO.ionFrac*gyroFrac*vPhaseFrac*vSliceFrac[:,vSlice][notNan]
+                    isNanFrom = np.where(np.isnan(idx)==True)[0] #include NaNs (NaNs = no intersection) index we map from
+                    notNanFrom = np.where(np.isnan(idx)==False)[0] #dont include NaNs (NaNs = no intersection) index we map from
+                    notNanTo = idx[~np.isnan(idx)] #indices we map power to
+                    notNanTo = idx1.astype(int) #cast as integer
+                    isNanTo = idx[np.isnan(idx)] #indices we map power to
+                    isNanTo = idx2.astype(int) #cast as integer
+
                     if len(notNan)>0:
-                        Pgyro[idx1] += Pdiv[notNan]*GYRO.ionFrac*gyroFrac*vPhaseFrac*vSliceFrac[notNan,vSlice]
-                        #Pgyro[idx1] += Pdiv[notNan]*GYRO.ionFrac*gyroFrac*vPhaseFrac*vSliceFrac[notNan,vSlice]
-                        #qGyro[idx1] += q[notNan]*GYRO.ionFrac*gyroFrac*vPhaseFrac*vSliceFrac[notNan,vSlice]
-                        #multiply by hdotn: incident angle of helix * face normal
-                        #qGyro[idx1] *= np.abs(GYRO.hdotn[gyroPhase,vPhase,vSlice,idx1])
+                        #multiple Froms can light up the same To, so we loop
+                        for i in range(len(notNan)):
+                            Pgyro[notNanTo[i]] += Pdiv[notNanFrom[i]]*GYRO.ionFrac*gyroFrac*vPhaseFrac*vSliceFrac[notNanFrom[i],vSlice]
+
                     if len(isNan)>0:
                         PNaN += np.sum(Pdiv[isNan]*GYRO.ionFrac*gyroFrac*vPhaseFrac*vSliceFrac[isNan,vSlice])
 
         #print("\nTEST2")
         #print(GYRO.intersectRecord[0,0,0,1711])
         #print(Pgyro[1711])
 
-
-
         GYRO.gyroPowMatrix += Pgyro
         GYRO.gyroNanPower += PNaN
         return

source/inputs/NSTXU/NSTXU_input.csv

@@ -11,7 +11,7 @@ gridRes, standard
 shot, 204118
 tmin, 0
 tmax, 5000
-nTrace, 10
+nTrace, 2
 dataPath, /home/tom/data/HEAT/
 plasma3Dmask, 0 #1=true, 0=false
 #==============================================================
@@ -38,7 +38,7 @@ fracLO,0.7
 #                Ion Gyro Orbit HF Variables
 #==============================================================
 N_gyroSteps, 5
-gyroDeg, 3
+gyroDeg, 1
 gyroT_eV, 100
 N_vSlice, 1
 N_vPhase, 1