updates

sandialabs · Jul 16, 2024 · 0fe207e · 0fe207e
1 parent 07ee407
commit 0fe207e
Show file tree

Hide file tree

Showing 29 changed files with 963 additions and 776 deletions.
diff --git a/.github/workflows/build.yml b/.github/workflows/build.yml
@@ -47,6 +47,7 @@ jobs:
       OMPI_CXX: ${{ matrix.cxx-compiler }}
       BUILD_PRETTY_IDENTIFIER: "Linux ${{ matrix.c-compiler }} Python ${{ matrix.py-version }}"
       BUILD_IDENTIFIER: "Linux-${{ matrix.c-compiler }}-${{ matrix.py-version }}"
+      OMP_NUM_THREADS: 1
 
     steps:
     - name: Check out repo
@@ -70,7 +71,7 @@ jobs:
       if: always()
       run: |
         sudo apt-get update
-        sudo apt-get install mpi-default-bin mpi-default-dev libmetis-dev libparmetis-dev libsuitesparse-dev ccache
+        sudo apt-get install -yq mpi-default-bin mpi-default-dev libmetis-dev libparmetis-dev libsuitesparse-dev ccache
 
     - name: Get ccache config dir
       if: always()
@@ -149,7 +150,7 @@ jobs:
       fail-fast: false
     name: Mac clang Python ${{ matrix.py-version }}
 
-    runs-on: macos-latest
+    runs-on: macos-12
     timeout-minutes: 180
     env:
       MPIEXEC_FLAGS: "--allow-run-as-root --oversubscribe"
@@ -158,6 +159,7 @@ jobs:
       SUITESPARSE_LIBRARY_DIR: /usr/local/Cellar/brewsci-suite-sparse/4.5.5_2/lib
       LDFLAGS: -L/usr/local/opt/brewsci-metis/lib -L/usr/local/opt/brewsci-parmetis/lib
       CPPFLAGS: -I/usr/local/opt/brewsci-metis/include -I/usr/local/opt/brewsci-parmetis/include
+      OMP_NUM_THREADS: 1
 
     steps:
     - name: Check out repo

diff --git a/base/PyNucleus_base/utilsFem.py b/base/PyNucleus_base/utilsFem.py
@@ -1356,7 +1356,7 @@ def diffDict(d1, d2, aTol, relTol):
     return diff
 
 
-def runDriver(path, py, python=None, timeout=600, ranks=None, cacheDir='',
+def runDriver(path, py, python=None, timeout=900, ranks=None, cacheDir='',
               overwriteCache=False,
               aTol=1e-12, relTol=1e-2, extra=None):
     from subprocess import Popen, PIPE, TimeoutExpired
@@ -1481,6 +1481,7 @@ def declareGeneratedProperties(self, props):
             self.declareGeneratedProperty(prop)
 
     def __call__(self):
+        from PyNucleus_fem.DoFMaps import fe_vector, complex_fe_vector
         cached_args = {}
         args = []
         needToBuild = False
@@ -1499,6 +1500,14 @@ def __call__(self):
                             needToBuild = True
                         else:
                             dependencyLogger.log(self.logLevel, 'Values for {} are identical: \'{}\' == \'{}\''.format(prop, oldValue, newValue))
+                    elif isinstance(newValue, (fe_vector, complex_fe_vector)):
+                        cached_args[prop] = newValue.copy()
+                        if oldValue is None or (newValue.toarray() != oldValue.toarray()).any():
+                            dependencyLogger.log(self.logLevel, 'Values for {} differ: \'{}\' != \'{}\', calling \'{}\''.format(prop, oldValue,
+                                                                                                                                newValue, self.fun.__name__))
+                            needToBuild = True
+                        else:
+                            dependencyLogger.log(self.logLevel, 'Values for {} are identical: \'{}\' == \'{}\''.format(prop, oldValue, newValue))
                     elif newValue != oldValue:
                         cached_args[prop] = newValue
                         dependencyLogger.log(self.logLevel, 'Values for {} differ: \'{}\' != \'{}\', calling \'{}\''.format(prop, oldValue,

diff --git a/fem/PyNucleus_fem/algebraicOverlaps.pyx b/fem/PyNucleus_fem/algebraicOverlaps.pyx
@@ -1312,7 +1312,7 @@ cdef class algebraicOverlapManager:
                 v[i] = m
                 m += 1
         self.unique(v)
-        return np.array(v, copy=False, dtype=INDEX)
+        return np.asarray(v, dtype=INDEX)
 
     def __repr__(self):
         s = ''

diff --git a/fem/PyNucleus_fem/distributed_operators_{SCALAR}.pxi b/fem/PyNucleus_fem/distributed_operators_{SCALAR}.pxi
@@ -127,7 +127,7 @@ cdef class {SCALAR_label}CSR_DistributedLinearOperator({SCALAR_label}Distributed
 
         k = 0
         for i in range(self.num_rows):
-            if self.overlap_indices[k] == i:
+            if k < self.overlap_indices.shape[0] and self.overlap_indices[k] == i:
                 k += 1
                 continue
             temp = 0.0
@@ -173,7 +173,7 @@ cdef class {SCALAR_label}CSR_DistributedLinearOperator({SCALAR_label}Distributed
 
         k = 0
         for i in range(self.num_rows):
-            if self.overlap_indices[k] == i:
+            if k < self.overlap_indices.shape[0] and self.overlap_indices[k] == i:
                 k += 1
                 continue
             temp = rhs[i]

diff --git a/fem/PyNucleus_fem/pdeProblems.py b/fem/PyNucleus_fem/pdeProblems.py
@@ -231,7 +231,7 @@ def n(x):
             else:
                 raise NotImplementedError(problem)
         elif domain == 'cube':
-            self.noRef = 5
+            self.noRef = 6
 
             def n(x):
                 if x[2] == 0:

diff --git a/multilevelSolver/PyNucleus_multilevelSolver/connectors.py b/multilevelSolver/PyNucleus_multilevelSolver/connectors.py
@@ -241,7 +241,7 @@ def build(self):
             with self.Timer('Building distribute from \'{}\' to \'{}\''.format(label1, self.label2)):
                 self.algOM.prepareDistributeRepartition(dm)
             if self.debugOverlaps:
-                from PyNucleus_fem import solSin1D, solSin2D, solSin3D
+                from PyNucleus_fem.factories import solSin1D, solSin2D, solSin3D
                 if subdomain.dim == 1:
                     xOld = dm.interpolate(solSin1D)
                 elif subdomain.dim == 2:
@@ -270,7 +270,7 @@ def build(self):
             with self.Timer('Building distribute from \'{}\' to \'{}\''.format(label1, self.label2)):
                 self.algOMnew.prepareDistributeRepartition(dmNew)
             if self.debugOverlaps:
-                from PyNucleus_fem import solSin1D, solSin2D, solSin3D
+                from PyNucleus_fem.factories import solSin1D, solSin2D, solSin3D
                 if subdomainNew.dim == 1:
                     xNew = dmNew.interpolate(solSin1D)
                 elif subdomainNew.dim == 2:
@@ -312,7 +312,7 @@ def build(self):
                 self.algOM.prepareDistributeRepartitionSend(dm)
                 self.algOMnew.prepareDistributeRepartition(dmNew, doSend=False)
             if self.debugOverlaps:
-                from PyNucleus_fem import solSin1D, solSin2D, solSin3D
+                from PyNucleus_fem.factories import solSin1D, solSin2D, solSin3D
                 if subdomain.dim == 1:
                     xOld = dm.interpolate(solSin1D)
                 elif subdomain.dim == 2:

diff --git a/multilevelSolver/PyNucleus_multilevelSolver/geometricMG.py b/multilevelSolver/PyNucleus_multilevelSolver/geometricMG.py
@@ -35,7 +35,7 @@ def paramsForSerialMG(noRef, global_params):
 
 
 def paramsForMG(noRef, onRanks, global_params, dim, element, repartitionFactor=0.05,
-                max_coarse_grid_size=5000):
+                max_coarse_grid_size=4500):
     from . connectors import repartitionConnector
 
     numProcsAvail = len(onRanks)
@@ -61,7 +61,7 @@ def paramsForMG(noRef, onRanks, global_params, dim, element, repartitionFactor=0
         else:
             raise NotImplementedError()
     elif dim == 3:
-        numInitialCells = 6
+        numInitialCells = 48
         if element in ('P1', 1):
             cells2dofsFactor = 1./6.
         elif element in ('P2', 2):

diff --git a/multilevelSolver/PyNucleus_multilevelSolver/hierarchies.py b/multilevelSolver/PyNucleus_multilevelSolver/hierarchies.py
@@ -154,7 +154,7 @@ def build(self):
                             multLvlAlgOvManager.levels[lvl].prepareDistribute()
 
                     if self.params['debugOverlaps']:
-                        from PyNucleus_fem import solSin1D, solSin2D, solSin3D
+                        from PyNucleus_fem.factories import solSin1D, solSin2D, solSin3D
                         for lvl in range(len(self.algebraicLevels)):
                             if self.algebraicLevels[lvl].DoFMap is not None:
                                 dm = self.algebraicLevels[lvl].DoFMap

diff --git a/nl/PyNucleus_nl/interactionDomains.pyx b/nl/PyNucleus_nl/interactionDomains.pyx
@@ -295,7 +295,7 @@ cdef class interactionDomain(parametrizedTwoPointFunction):
         return interactionDomain, (self.horizonFun, self.complement, self.symmetric)
 
     def getLongDescription(self):
-        return ''
+        return 'MISSING DESCRIPTION'
 
 
 cdef class barycenterDomain(interactionDomain):
@@ -659,7 +659,8 @@ cdef class retriangulationDomain(interactionDomain):
                         self.A_Node[self.iterEnd_Node, outside2, outside2] = c2
                         self.vol_Node[self.iterEnd_Node] = bary[outside1]*c2+bary[outside2]*c1-c1*c2
                         self.iterEnd_Node += 1
-                else:
+                elif numIntersections == 2:
+
                     self.A_Node[self.iterEnd_Node, :, :] = 0.
                     self.A_Node[self.iterEnd_Node, inside, inside] = 1
                     self.A_Node[self.iterEnd_Node, outside1, outside1] = c1
@@ -687,6 +688,26 @@ cdef class retriangulationDomain(interactionDomain):
                     self.vol_Node[self.iterEnd_Node] = c2*(1-intersections[1])
                     self.iterEnd_Node += 1
 
+                else:
+
+                    self.A_Node[self.iterEnd_Node, :, :] = 0.
+                    self.A_Node[self.iterEnd_Node, inside, inside] = 1
+                    self.A_Node[self.iterEnd_Node, outside1, outside1] = c1
+                    self.A_Node[self.iterEnd_Node, inside, outside1] = 1-c1
+                    self.A_Node[self.iterEnd_Node, outside2, outside2] = intersections[0]
+                    self.A_Node[self.iterEnd_Node, outside1, outside2] = 1-intersections[0]
+                    self.vol_Node[self.iterEnd_Node] = c1*intersections[0]
+                    self.iterEnd_Node += 1
+
+                    self.A_Node[self.iterEnd_Node, :, :] = 0.
+                    self.A_Node[self.iterEnd_Node, inside, inside] = 1
+                    self.A_Node[self.iterEnd_Node, outside1, outside1] = 1-intersections[0]
+                    self.A_Node[self.iterEnd_Node, outside2, outside1] = intersections[0]
+                    self.A_Node[self.iterEnd_Node, outside2, outside2] = c2
+                    self.A_Node[self.iterEnd_Node, inside, outside2] = 1-c2
+                    self.vol_Node[self.iterEnd_Node] = c2*(1-intersections[0])
+                    self.iterEnd_Node += 1
+
             elif numInside == 2:
                 outside = 0
                 while ind[outside]:
@@ -1561,6 +1582,9 @@ cdef class ellipse_barycenter(linearTransformInteraction):
     def __reduce__(self):
         return ellipse_barycenter, (self.horizonFun, self.a, self.b, self.theta)
 
+    def getLongDescription(self):
+        return '\\chi_{ellipse(x,y;a,b,\\theta)<\\delta}'
+
 
 cdef class ball1_retriangulation(linearTransformInteraction):
     "l1 ball interaction domain"
@@ -1662,7 +1686,7 @@ cdef class ball1_barycenter(linearTransformInteraction):
 
 cdef class ball2_dilation_barycenter(barycenterDomain):
     def __init__(self, sqrtAffineFunction horizonFun):
-        super(ball2_dilation_barycenter, self).__init__(horizonFun, False, False)
+        super(ball2_dilation_barycenter, self).__init__(horizonFun, False, True)
         self.c = horizonFun.c
         self.d = np.linalg.norm(horizonFun.w)
         self.w = np.array(horizonFun.w)/self.d
@@ -1799,10 +1823,13 @@ cdef class ball2_dilation_barycenter(barycenterDomain):
     def __reduce__(self):
         return ballInf_retriangulation, (self.horizonFun, )
 
+    def getLongDescription(self):
+        return '\\chi_{|x-y|_2<\\delta(x)} \\chi_{|x-y|_2<\\delta(y)}'
+
 
 cdef class ball2_dilation_retriangulation(retriangulationDomain):
     def __init__(self, sqrtAffineFunction horizonFun):
-        super(ball2_dilation_retriangulation, self).__init__(horizonFun, False, False)
+        super(ball2_dilation_retriangulation, self).__init__(horizonFun, False, True)
         self.c = horizonFun.c
         self.d = np.linalg.norm(horizonFun.w)
         self.w = np.array(horizonFun.w)/self.d
@@ -1946,11 +1973,11 @@ cdef class ball2_dilation_retriangulation(retriangulationDomain):
         A = -p*0.5
         B = sqrt(A**2-q)
         c = A-B
-        if (c >= 0) and (c <= 1):
+        if (c >= -1e-12) and (c <= 1+1e-12):
             intersections[numIntersections] = c
             numIntersections += 1
         c = A+B
-        if (c >= 0) and (c <= 1):
+        if (c >= -1e-12) and (c <= 1+1e-12):
             intersections[numIntersections] = c
             numIntersections += 1
         return numIntersections
@@ -2030,3 +2057,6 @@ cdef class ball2_dilation_retriangulation(retriangulationDomain):
                     break
             if doExit:
                 break
+
+    def getLongDescription(self):
+        return '\\chi_{|x-y|_2<\\delta(x)} \\chi_{|x-y|_2<\\delta(y)}'
diff --git a/nl/PyNucleus_nl/kernel_params.pxi b/nl/PyNucleus_nl/kernel_params.pxi
@@ -95,3 +95,4 @@ cpdef enum:
     GREENS_3D = 7
     EXPONENTIAL = 8
     POLYNOMIAL = 9
+    PERIODIC_FRACTIONAL = 10
diff --git a/nl/PyNucleus_nl/kernels.py b/nl/PyNucleus_nl/kernels.py
@@ -157,6 +157,12 @@ def getFractionalKernel(dim,
                     phi = fac
         kernel = FractionalKernel(dim_, sFun, horizonFun, interaction, scaling, phi, piecewise=piecewise, boundary=boundary,
                                   derivative=derivative, tempered=tempered, max_horizon=max_horizon)
+
+    from . twoPointFunctions import parametrizedTwoPointFunction
+    if isinstance(kernel.scaling, parametrizedTwoPointFunction):
+        assert kernel.getParamPtrAddr() == kernel.scaling.getParamPtrAddr()
+    if isinstance(kernel.interaction, parametrizedTwoPointFunction):
+        assert kernel.getParamPtrAddr() == kernel.interaction.getParamPtrAddr()
     return kernel