gyro diamagnetism correction

AppImageBuilder.yml

@@ -100,7 +100,7 @@ AppDir:
     id: HEATAppimage
     name: HEAT_AppImage
     icon: flame
-    version: v1.4.0
+    version: v1.4.1
 #    # Set the python executable as entry point
 #    exec: usr/bin/python3
 #    # Set the application main script path as argument. Use '$@' to forward CLI parameters

source/GUIclass.py

@@ -931,6 +931,7 @@ def bfieldMagnitude(self, PFC):
         PFC.Bmag[:,1] = PFC.centers[:,1] # Y
         PFC.Bmag[:,2] = PFC.centers[:,2] # Z
         PFC.Bmag[:,3] = np.sqrt(PFC.Bxyz[:,0]**2+PFC.Bxyz[:,1]**2+PFC.Bxyz[:,2]**2)
+        PFC.Bsign = np.sign(PFC.ep.g['Bt0'])
         return
 
     def BtraceMultiple(self,data,t):
@@ -1039,6 +1040,7 @@ def gyroTrace(self,x,y,z,t,gyroPhase,gyroDeg,N_gyroSteps,gyroDir,gyroT_eV):
         BR = ep.BRFunc.ev(R,Z)
         BZ = ep.BZFunc.ev(R,Z)
         B = np.sqrt(BR**2 + Bt**2 + BZ**2)
+        Bmag = np.sign(ep.g['Bt0'])
         #Calculate frequencies and gyro radius
         self.GYRO.setupFreqs(B)
         self.GYRO.setupRadius(vPerp)
@@ -1122,9 +1124,12 @@ def loadAccFilters(self, accFilters):
             phiFilterSwitch = False
 
         #for now all PFCs share a single filter state
-        for PFC in self.PFCs:
-            PFC.psiFilterSwitch = psiFilterSwitch
-            PFC.phiFilterSwitch = phiFilterSwitch
+        try:
+            for PFC in self.PFCs:
+                PFC.psiFilterSwitch = psiFilterSwitch
+                PFC.phiFilterSwitch = phiFilterSwitch
+        except:
+            print("No PFCs to set filter switches for.  Skipping.")
         return
 
 
@@ -2149,12 +2154,14 @@ def runOpenFOAM(self):
         #set up OF parts for each PFC part
         for PFC in self.PFCs:
             #check if PFC is a gyroSource Plane
-            if self.GYRO.gyroSourceTag == 'allROI':
-                pass
-            else:
-                if PFC.name in self.GYRO.gyroSources:
-                    print("Not including "+PFC.name+" in thermal analysis")
-                    continue
+            if hasattr(self.GYRO, 'gyroSourceTag'):
+                #check if PFC is a gyroSource Plane
+                if self.GYRO.gyroSourceTag == 'allROI':
+                    pass
+                else:
+                    if PFC.name in self.GYRO.gyroSources:
+                        print("Not including "+PFC.name+" in thermal analysis")
+                        continue
 
             print("Running openFOAM for PFC: "+PFC.name)
             log.info("Running openFOAM for PFC: "+PFC.name)
@@ -2345,6 +2352,20 @@ def runOpenFOAM(self):
             qDiv = np.zeros((len(PFC.centers)))
             ctrs = copy.copy(PFC.centers)*1000.0
 
+
+            #in the future if you want adaptive meshing for each timestep this whole
+            #section needs to be revised.  for now, you can only have 1 boundary
+            #condition set of points.  We take timestep 0 from the PFC
+            #and get the points from a file.  This enables user to swap out
+            #HF_allSources.csv in between HF and Temp runs, to create non-uniform mesh
+            #get the mesh centers from one of the timesteps
+            t = self.MHD.timesteps[0]
+            if self.MHD.shotPath[-1]=='/':
+                HFcsv = self.MHD.shotPath + '{:06d}/'.format(t) + PFC.name + '/HF_allSources.csv'
+            else:
+                HFcsv = self.MHD.shotPath + '/' + '{:06d}/'.format(t) + PFC.name + '/HF.allSources.csv'
+            OFcenters = pd.read_csv(HFcsv).iloc[:,0:3].values
+
             #cycle through timesteps and get HF data from HEAT tree
             for t in OFtimesteps:
                 print("openFOAM timestep: {:d}".format(t))
@@ -2380,17 +2401,18 @@ def runOpenFOAM(self):
                     else:
                         HFcsv = self.MHD.shotPath + '/' + '{:06d}/'.format(t) + PFC.name + '/HF.allSources.csv'
                     qDiv = pd.read_csv(HFcsv)['HeatFlux'].values
+                    #OFcenters = pd.read_csv(HFcsv).iloc[:,0:3].values
                     #write boundary condition
                     print("Maximum qDiv for this PFC and time: {:f}".format(qDiv.max()))
-                    self.HF.write_openFOAM_boundary(ctrs,qDiv,partDir,OFt)
+                    self.HF.write_openFOAM_boundary(OFcenters,qDiv,partDir,OFt)
                 elif (t < self.MHD.timesteps.min()) or (t > self.MHD.timesteps.max()):
                     #apply zero HF outside of discharge domain (ie tiles cooling)
                     print("OF.timestep: {:d} outside MHD domain".format(t))
                     log.info("OF.timestep: {:d} outside MHD domain".format(t))
-                    qDiv = np.zeros((len(ctrs)))
+                    qDiv = np.zeros((len(OFcenters)))
                     #write boundary condition
                     print("Maximum qDiv for this PFC and time: {:f}".format(qDiv.max()))
-                    self.HF.write_openFOAM_boundary(ctrs,qDiv,partDir,OFt)
+                    self.HF.write_openFOAM_boundary(OFcenters,qDiv,partDir,OFt)
                 else:
                     #boundary using last timestep that we calculated a HF for
                     #(basically a heaviside function in time)
@@ -2444,12 +2466,13 @@ def getOFMinMaxPlots(self):
         pfcNames = []
         for PFC in self.PFCs:
             #check if PFC is a gyroSource Plane
-            if self.GYRO.gyroSourceTag == 'allROI':
-                pass
-            else:
-                if PFC.name in self.GYRO.gyroSources:
-                    print("Not including "+PFC.name+" in MinMax plots")
-                    continue
+            if hasattr(self.GYRO, 'gyroSourceTag'):
+                if self.GYRO.gyroSourceTag == 'allROI' :
+                    pass
+                else:
+                    if PFC.name in self.GYRO.gyroSources:
+                        print("Not including "+PFC.name+" in MinMax plots")
+                        continue
 
             if self.MHD.shotPath[-1]=='/':
                 partDir = self.MHD.shotPath + 'openFoam/heatFoam/'+PFC.name

source/gyroClass.py

@@ -102,6 +102,7 @@ def setupConstants(self, ionMassAMU=2.014):
         self.kB = 8.617e-5 #ev/K
         self.e = 1.602e-19 # C
         self.c = 299792458 #m/s
+        self.diamag = -1 #diamagnetism = -1 for ions, 1 for electrons
 
         self.mass_eV = ionMassAMU * self.AMU
         self.Z=1 #assuming isotopes of hydrogen here
@@ -134,7 +135,7 @@ def setupFreqs(self, B):
         if np.isscalar(self.omegaGyro):
             self.omegaGyro = np.array([self.omegaGyro])
 
-        self.fGyro = self.omegaGyro/(2*np.pi)
+        self.fGyro = np.abs(self.omegaGyro)/(2*np.pi)
         self.TGyro = 1.0/self.fGyro
         return
 
@@ -326,7 +327,7 @@ def singleGyroTrace(self,vPerp,vParallel,gyroPhase,N_gyroSteps,
             xfm = np.vstack([u,v,w]).T
             #get helix path along (proxy) z axis reference frame
             x_helix = self.rGyro*np.cos(self.omegaGyro*t + gyroPhase)
-            y_helix = self.rGyro*np.sin(self.omegaGyro*t + gyroPhase)
+            y_helix = self.diamag*self.rGyro*np.sin(self.omegaGyro*t + gyroPhase)
             z_helix = np.zeros((len(t)))
             #perform rotation to field line reference frame
             helix = np.vstack([x_helix,y_helix,z_helix]).T
@@ -415,7 +416,7 @@ def gyroTraceParallel(self, i, mode='MT'):
         omega = self.omegaGyro[self.GYRO_HLXmap][i]
         theta = self.lastPhase[self.GYRO_HLXmap][i]
         x_helix = rGyro*np.cos(omega*t + theta)
-        y_helix = rGyro*np.sin(omega*t + theta)
+        y_helix = self.diamag*rGyro*np.sin(omega*t + theta)
         z_helix = np.zeros((len(t)))
         #perform rotation to field line reference frame
         helix = np.vstack([x_helix,y_helix,z_helix]).T

source/openFOAMclass.py

@@ -195,9 +195,9 @@ def writeOFtemplateVarFile(self, file, STLpart):
         """
         #determine number of logical cpus and leave 1 for overhead
         if psutil.cpu_count(logical=True) < 2:
-            NCPU = 1
+            self.NCPU = 1
         else:
-            NCPU = psutil.cpu_count(logical=True) - 1
+            self.NCPU = psutil.cpu_count(logical=True) - 1
         with open(file, 'w') as f:
             f.write('xMin {:f};\n'.format(self.xMin))
             f.write('xMax {:f};\n'.format(self.xMax))
@@ -221,7 +221,7 @@ def writeOFtemplateVarFile(self, file, STLpart):
             f.write('zMid {:f};\n'.format(self.zMid))
             f.write('STLfileName '+STLpart+';\n')
             f.write('STLlayerName '+STLpart+'_firstSolid;\n')
-            f.write('NCPU '+str(NCPU)+';\n')
+            f.write('NCPU '+str(self.NCPU)+';\n')
         return
 
     def writeShellScript(self,logFile, parallel=False):
@@ -251,10 +251,11 @@ def writeShellScript(self,logFile, parallel=False):
             self.partDir += '/'
 
         #determine number of logical cpus and leave 1 for overhead
-        if psutil.cpu_count(logical=True) < 2:
-            NCPU = 1
-        else:
-            NCPU = psutil.cpu_count(logical=True) - 1
+        #now done in template file creation method
+        #if psutil.cpu_count(logical=True) < 2:
+        #    self.NCPU = 1
+        #else:
+        #    self.NCPU = psutil.cpu_count(logical=True) - 1
 
 
 #        #write openfoam sourcing script to generate environment variables
@@ -294,12 +295,10 @@ def writeShellScript(self,logFile, parallel=False):
             else:
                 #Run snappyHexMesh across multiple processors
                 #this feature should be added in one day...sigh...for now its a placeholder
-                #you'll have to set up blocks in blockMeshDict or this will not work.
-                #problem is you need to allocate processors in advance.
-                #see: "2.3.11 Running in parallel"
-                #from https://cfd.direct/openfoam/user-guide/v8-damBreak/#x7-610002.3.11
+                #you need to download libscotch-6.0 from ubuntu repo,
+                #then build scotchDecomp from source from src/parallel/decompose/scotchDecomp
                 f.write('decomposePar > ' + logFile + '\n')
-                f.write('mpirun --use-hwthread-cpus -np '+str(NCPU)+' snappyHexMesh -parallel -overwrite > '+logFile+ '\n')
+                f.write('mpirun --use-hwthread-cpus -np '+str(self.NCPU)+' snappyHexMesh -parallel -overwrite > '+logFile+ '\n')
                 f.write('reconstructParMesh -mergeTol 1e-6 -latestTime -constant > '+logFile+ '\n')
 
         os.chmod(file, 0o744)
@@ -427,6 +426,8 @@ def runThermalAnalysis(self):
         """
         print('Running Thermal Analysis')
         log.info('Running Thermal Analysis')
+        print('See HEAT LogFile Tab for Status')
+        log.info('See HEAT LogFile Tab for Status')
 
         from subprocess import run
         #Copy the current environment

source/openFoamTemplates/templateDicts/decomposeParDict.template

@@ -14,14 +14,16 @@ FoamFile
 }
 // * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
 
-numberOfSubdomains  8;
+numberOfSubdomains  |-NCPU-|;
 
-method          simple;
+method          scotch;
 
-simpleCoeffs
-{
-    n               (2 2 2);
-    delta           0.001;
-}
+//OLD method
+//method          simple;
+//simpleCoeffs
+//{
+//    n               (|-NCPU-| 1 1);
+//    delta           0.001;
+//}
 
 // ************************************************************************* //

source/openFoamTemplates/templateDicts/topoSetDict.template

@@ -24,7 +24,8 @@ actions
         source  patchToFace;
         sourceInfo
         {
-            name |-STLfileName-|;
+            //name |-STLfileName-|; //OLD METHOD
+            patch |-STLfileName-|;
         }
     }
 );