From 3fc3b98e96b87f976f8acbdb40ff35bd348ac8ff Mon Sep 17 00:00:00 2001 From: Christian Glusa Date: Tue, 18 Jul 2023 17:04:34 -0600 Subject: [PATCH] updates --- .../CSR_LinearOperator_decl_{SCALAR}.pxi | 2 + .../CSR_LinearOperator_{SCALAR}.pxi | 2 + .../DenseLinearOperator_decl_{SCALAR}.pxi | 11 + .../DenseLinearOperator_{SCALAR}.pxi | 87 +++ .../DiagonalLinearOperator_{SCALAR}.pxi | 9 + .../LinearOperator_decl_{SCALAR}.pxi | 15 + .../LinearOperator_{SCALAR}.pxi | 101 +++ .../SSS_LinearOperator_decl_{SCALAR}.pxi | 2 + .../SSS_LinearOperator_{SCALAR}.pxi | 2 + base/PyNucleus_base/utilsFem.py | 1 + drivers/runNonlocal.py | 3 +- drivers/runParallelGMG.py | 2 +- fem/PyNucleus_fem/DoFMaps.pyx | 33 - fem/PyNucleus_fem/factories.py | 2 +- fem/PyNucleus_fem/functions.pyx | 43 ++ fem/PyNucleus_fem/lookupFunction.pxd | 19 + fem/PyNucleus_fem/lookupFunction.pyx | 35 ++ fem/PyNucleus_fem/meshCy.pxd | 6 + fem/PyNucleus_fem/meshCy.pyx | 80 +++ fem/setup.py | 2 + .../PyNucleus_multilevelSolver/hierarchies.py | 2 + nl/PyNucleus_nl/clusterMethodCy.pyx | 2 +- nl/PyNucleus_nl/discretizedProblems.py | 201 +++--- nl/PyNucleus_nl/fractionalLaplacian1D.pyx | 20 +- nl/PyNucleus_nl/fractionalLaplacian2D.pyx | 10 +- nl/PyNucleus_nl/fractionalOrders.pxd | 69 +-- nl/PyNucleus_nl/fractionalOrders.pyx | 576 +++++++----------- nl/PyNucleus_nl/helpers.py | 11 +- nl/PyNucleus_nl/kernelNormalization.pxd | 86 +++ nl/PyNucleus_nl/kernelNormalization.pyx | 418 +++++++++++++ nl/PyNucleus_nl/kernels.py | 21 +- nl/PyNucleus_nl/kernelsCy.pxd | 6 +- nl/PyNucleus_nl/kernelsCy.pyx | 326 ++++++++-- nl/PyNucleus_nl/nonlocalLaplacian.pyx | 261 +++++++- nl/PyNucleus_nl/nonlocalLaplacianBase.pxd | 4 + nl/PyNucleus_nl/nonlocalLaplacianBase.pyx | 77 ++- nl/PyNucleus_nl/nonlocalProblems.py | 199 ++++-- nl/setup.py | 2 + ...nt--elementP0--solvercg-mg--matrixFormatH2 | 16 +- ...-elementP0--solvercg-mg--matrixFormatdense | 16 +- ...nt--elementP1--solvercg-mg--matrixFormatH2 | 16 +- ...-elementP1--solvercg-mg--matrixFormatdense | 16 +- ...mentP1--solvergmres-jacobi--matrixFormatH2 | 16 +- ...tP1--solvergmres-jacobi--matrixFormatdense | 16 +- ...nt--elementP0--solvercg-mg--matrixFormatH2 | 16 +- ...-elementP0--solvercg-mg--matrixFormatdense | 16 +- ...nt--elementP1--solvercg-mg--matrixFormatH2 | 16 +- ...-elementP1--solvercg-mg--matrixFormatdense | 16 +- ...nt--elementP2--solvercg-mg--matrixFormatH2 | 16 +- ...-elementP2--solvercg-mg--matrixFormatdense | 16 +- ...nt--elementP3--solvercg-mg--matrixFormatH2 | 16 +- ...-elementP3--solvercg-mg--matrixFormatdense | 16 +- ...elementP1--solvercg-jacobi--matrixFormatH2 | 12 +- ...oFlux--elementP1--solverlu--matrixFormatH2 | 12 +- ...mentP1--solvergmres-jacobi--matrixFormatH2 | 16 +- ...tP1--solvergmres-jacobi--matrixFormatdense | 16 +- ...nt--elementP0--solvercg-mg--matrixFormatH2 | 6 +- ...-elementP0--solvercg-mg--matrixFormatdense | 6 +- ...nt--elementP1--solvercg-mg--matrixFormatH2 | 6 +- ...-elementP1--solvercg-mg--matrixFormatdense | 6 +- ...mentP1--solvergmres-jacobi--matrixFormatH2 | 6 +- ...tP1--solvergmres-jacobi--matrixFormatdense | 6 +- ...nt--elementP0--solvercg-mg--matrixFormatH2 | 6 +- ...-elementP0--solvercg-mg--matrixFormatdense | 6 +- ...nt--elementP1--solvercg-mg--matrixFormatH2 | 6 +- ...-elementP1--solvercg-mg--matrixFormatdense | 6 +- ...nt--elementP2--solvercg-mg--matrixFormatH2 | 6 +- ...-elementP2--solvercg-mg--matrixFormatdense | 6 +- ...nt--elementP3--solvercg-mg--matrixFormatH2 | 6 +- ...-elementP3--solvercg-mg--matrixFormatdense | 6 +- ...elementP1--solvercg-jacobi--matrixFormatH2 | 6 +- ...oFlux--elementP1--solverlu--matrixFormatH2 | 6 +- ...mentP1--solvergmres-jacobi--matrixFormatH2 | 6 +- ...tP1--solvergmres-jacobi--matrixFormatdense | 6 +- ...empoly-Dirichlet--solverlu--matrixFormatH2 | 9 + ...blempoly-Neumann--solverlu--matrixFormatH2 | 9 + ...empoly-Dirichlet--solverlu--matrixFormatH2 | 9 + ...blempoly-Neumann--solverlu--matrixFormatH2 | 9 + ...empoly-Dirichlet--solverlu--matrixFormatH2 | 9 + ...blempoly-Neumann--solverlu--matrixFormatH2 | 9 + ...oly-Dirichlet--solvercg-mg--matrixFormatH2 | 9 + ...-Dirichlet--solvercg-mg--matrixFormatdense | 8 +- ...oly-Dirichlet--solvercg-mg--matrixFormatH2 | 9 + ...oly-Dirichlet--solvercg-mg--matrixFormatH2 | 9 + ...-Dirichlet--solvercg-mg--matrixFormatdense | 8 +- tests/test_drivers_intFracLapl.py | 42 +- tests/test_fracLapl.py | 4 +- tests/test_h2finiteHorizon.py | 175 +++--- tests/test_kernels.py | 366 +++++++++-- 89 files changed, 2795 insertions(+), 1025 deletions(-) create mode 100644 fem/PyNucleus_fem/lookupFunction.pxd create mode 100644 fem/PyNucleus_fem/lookupFunction.pyx create mode 100644 nl/PyNucleus_nl/kernelNormalization.pxd create mode 100644 nl/PyNucleus_nl/kernelNormalization.pyx create mode 100644 tests/cache_runNonlocal.py--domaininterval--kernelTypeconstant--problempoly-Dirichlet--solverlu--matrixFormatH2 create mode 100644 tests/cache_runNonlocal.py--domaininterval--kernelTypeconstant--problempoly-Neumann--solverlu--matrixFormatH2 create mode 100644 tests/cache_runNonlocal.py--domaininterval--kernelTypefractional--problempoly-Dirichlet--solverlu--matrixFormatH2 create mode 100644 tests/cache_runNonlocal.py--domaininterval--kernelTypefractional--problempoly-Neumann--solverlu--matrixFormatH2 create mode 100644 tests/cache_runNonlocal.py--domaininterval--kernelTypeinverseDistance--problempoly-Dirichlet--solverlu--matrixFormatH2 create mode 100644 tests/cache_runNonlocal.py--domaininterval--kernelTypeinverseDistance--problempoly-Neumann--solverlu--matrixFormatH2 create mode 100644 tests/cache_runNonlocal.py--domainsquare--kernelTypeconstant--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 create mode 100644 tests/cache_runNonlocal.py--domainsquare--kernelTypefractional--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 create mode 100644 tests/cache_runNonlocal.py--domainsquare--kernelTypeinverseDistance--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 diff --git a/base/PyNucleus_base/CSR_LinearOperator_decl_{SCALAR}.pxi b/base/PyNucleus_base/CSR_LinearOperator_decl_{SCALAR}.pxi index a463b77..49b14aa 100644 --- a/base/PyNucleus_base/CSR_LinearOperator_decl_{SCALAR}.pxi +++ b/base/PyNucleus_base/CSR_LinearOperator_decl_{SCALAR}.pxi @@ -26,3 +26,5 @@ cdef class {SCALAR_label}CSR_LinearOperator({SCALAR_label}LinearOperator): cdef {SCALAR}_t getEntry(self, INDEX_t I, INDEX_t J) cdef void setEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t val) cpdef {SCALAR_label}CSR_LinearOperator getBlockDiagonal(self, sparseGraph blocks) + + diff --git a/base/PyNucleus_base/CSR_LinearOperator_{SCALAR}.pxi b/base/PyNucleus_base/CSR_LinearOperator_{SCALAR}.pxi index 26bea96..e73e281 100644 --- a/base/PyNucleus_base/CSR_LinearOperator_{SCALAR}.pxi +++ b/base/PyNucleus_base/CSR_LinearOperator_{SCALAR}.pxi @@ -472,3 +472,5 @@ cdef BOOL_t sort_indices{SCALAR_label}(INDEX_t[::1] indptr, kk -= 1 indices[kk] = j return wasSorted + + diff --git a/base/PyNucleus_base/DenseLinearOperator_decl_{SCALAR}.pxi b/base/PyNucleus_base/DenseLinearOperator_decl_{SCALAR}.pxi index 0b4316e..634b48d 100644 --- a/base/PyNucleus_base/DenseLinearOperator_decl_{SCALAR}.pxi +++ b/base/PyNucleus_base/DenseLinearOperator_decl_{SCALAR}.pxi @@ -26,3 +26,14 @@ cdef class {SCALAR_label}Dense_SubBlock_LinearOperator({SCALAR_label}LinearOpera cdef: dict lookupI, lookupJ public {SCALAR}_t[:, :] data + + +cdef class {SCALAR_label}Dense_VectorLinearOperator({SCALAR_label}VectorLinearOperator): + cdef: + public {SCALAR}_t[:, :, ::1] data + + cdef INDEX_t matvec(self, + {SCALAR}_t[::1] x, + {SCALAR}_t[:, ::1] y) except -1 + cdef void getEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val) + cdef void setEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val) diff --git a/base/PyNucleus_base/DenseLinearOperator_{SCALAR}.pxi b/base/PyNucleus_base/DenseLinearOperator_{SCALAR}.pxi index c13403b..d5a620f 100644 --- a/base/PyNucleus_base/DenseLinearOperator_{SCALAR}.pxi +++ b/base/PyNucleus_base/DenseLinearOperator_{SCALAR}.pxi @@ -156,3 +156,90 @@ cdef class {SCALAR_label}Dense_SubBlock_LinearOperator({SCALAR_label}LinearOpera j = self.lookupJ.get(J, -1) if i >= 0 and j >= 0: self.data[i, j] += val + + +cdef class {SCALAR_label}Dense_VectorLinearOperator({SCALAR_label}VectorLinearOperator): + def __init__(self, + {SCALAR}_t[:, :, ::1] data): + {SCALAR_label}VectorLinearOperator.__init__(self, + data.shape[0], + data.shape[1], + data.shape[2]) + self.data = data + + cdef INDEX_t matvec(self, + {SCALAR}_t[::1] x, + {SCALAR}_t[:, ::1] y) except -1: + cdef: + INDEX_t i, j, k + for i in range(self.num_rows): + for k in range(self.vectorSize): + y[i, k] = 0. + for j in range(self.num_columns): + for k in range(self.vectorSize): + y[i, k] += self.data[i, j, k]*x[j] + return 0 + + cdef INDEX_t matvec_no_overwrite(self, + {SCALAR}_t[::1] x, + {SCALAR}_t[:, ::1] y) except -1: + cdef: + INDEX_t i, j, k + for i in range(self.num_rows): + for j in range(self.num_columns): + for k in range(self.vectorSize): + y[i, k] += self.data[i, j, k]*x[j] + return 0 + + def isSparse(self): + return False + + def toarray(self): + return np.array(self.data, copy=False, dtype={SCALAR}) + + cdef void getEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val): + cdef: + INDEX_t k + for k in range(self.vectorSize): + val[k] = self.data[I, J, k] + + cdef void setEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val): + cdef: + INDEX_t k + for k in range(self.vectorSize): + self.data[I, J, k] = val[k] + + cdef void addToEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val): + cdef: + INDEX_t k + for k in range(self.vectorSize): + self.data[I, J, k] += val[k] + + @staticmethod + def zeros(INDEX_t num_rows, INDEX_t num_columns, INDEX_t vectorSize): + return Dense_LinearOperator(np.zeros((num_rows, num_columns, vectorSize), dtype={SCALAR})) + + @staticmethod + def ones(INDEX_t num_rows, INDEX_t num_columns, INDEX_t vectorSize): + return Dense_LinearOperator(np.ones((num_rows, num_columns, vectorSize), dtype={SCALAR})) + + @staticmethod + def empty(INDEX_t num_rows, INDEX_t num_columns, INDEX_t vectorSize): + return Dense_LinearOperator(uninitialized((num_rows, num_columns, vectorSize), dtype={SCALAR})) + + def getMemorySize(self): + return self.data.shape[0]*self.data.shape[1]*self.data.shape[2]*sizeof({SCALAR}_t) + + def __repr__(self): + sizeInMB = self.getMemorySize() >> 20 + if sizeInMB > 100: + return '<%dx%dx%d %s, %d MB>' % (self.num_rows, + self.num_columns, + self.vectorSize, + self.__class__.__name__, + sizeInMB) + else: + return '<%dx%dx%d %s>' % (self.num_rows, + self.num_columns, + self.vectorSize, + self.__class__.__name__) diff --git a/base/PyNucleus_base/DiagonalLinearOperator_{SCALAR}.pxi b/base/PyNucleus_base/DiagonalLinearOperator_{SCALAR}.pxi index 9489cb8..4738631 100644 --- a/base/PyNucleus_base/DiagonalLinearOperator_{SCALAR}.pxi +++ b/base/PyNucleus_base/DiagonalLinearOperator_{SCALAR}.pxi @@ -22,6 +22,15 @@ cdef class {SCALAR_label}diagonalOperator({SCALAR_label}LinearOperator): y[i] = self.data[i]*x[i] return 0 + cdef INDEX_t matvec_no_overwrite(self, + {SCALAR}_t[::1] x, + {SCALAR}_t[::1] y) except -1: + cdef: + INDEX_t i + for i in range(self.num_rows): + y[i] += self.data[i]*x[i] + return 0 + cdef {SCALAR}_t getEntry(self, INDEX_t I, INDEX_t J): if I == J: return self.data[I] diff --git a/base/PyNucleus_base/LinearOperator_decl_{SCALAR}.pxi b/base/PyNucleus_base/LinearOperator_decl_{SCALAR}.pxi index cef57d2..0ed3737 100644 --- a/base/PyNucleus_base/LinearOperator_decl_{SCALAR}.pxi +++ b/base/PyNucleus_base/LinearOperator_decl_{SCALAR}.pxi @@ -92,3 +92,18 @@ cdef class {SCALAR_label}Product_Linear_Operator({SCALAR_label}LinearOperator): {SCALAR}_t[::1] rhs, {SCALAR}_t[::1] result, BOOL_t simpleResidual=*) + + +cdef class {SCALAR_label}VectorLinearOperator: + cdef: + public INDEX_t num_rows, num_columns, vectorSize + {SCALAR}_t[::1] _diagonal + cdef INDEX_t matvec(self, + {SCALAR}_t[::1] x, + {SCALAR}_t[:, ::1] y) except -1 + cdef INDEX_t matvec_no_overwrite(self, + {SCALAR}_t[::1] x, + {SCALAR}_t[:, ::1] y) except -1 + cdef void addToEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val) + cdef void getEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val) + cdef void setEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val) diff --git a/base/PyNucleus_base/LinearOperator_{SCALAR}.pxi b/base/PyNucleus_base/LinearOperator_{SCALAR}.pxi index 39aee63..0457328 100644 --- a/base/PyNucleus_base/LinearOperator_{SCALAR}.pxi +++ b/base/PyNucleus_base/LinearOperator_{SCALAR}.pxi @@ -494,3 +494,104 @@ cdef class {SCALAR_label}Product_Linear_Operator({SCALAR_label}LinearOperator): def __repr__(self): return '{}*{}'.format(self.A, self.B) + + +cdef class {SCALAR_label}VectorLinearOperator: + def __init__(self, int num_rows, int num_columns, int vectorSize): + self.num_rows = num_rows + self.num_columns = num_columns + self.vectorSize = vectorSize + + cdef INDEX_t matvec(self, + {SCALAR}_t[::1] x, + {SCALAR}_t[:, ::1] y) except -1: + return -1 + + cdef INDEX_t matvec_no_overwrite(self, + {SCALAR}_t[::1] x, + {SCALAR}_t[:, ::1] y) except -1: + return -1 + + def __call__(self, + {SCALAR}_t[::1] x, + {SCALAR}_t[:, ::1] y, + BOOL_t no_overwrite=False): + if no_overwrite: + self.matvec_no_overwrite(x, y) + else: + self.matvec(x, y) + + property shape: + def __get__(self): + return (self.num_rows, self.num_columns, self.vectorSize) + + def isSparse(self): + raise NotImplementedError() + + def to_csr(self): + raise NotImplementedError() + + def to_dense(self): + return Dense_LinearOperator(self.toarray()) + + def toarray(self): + return self.to_csr().toarray() + + def toLinearOperator(self): + def matvec(x): + if x.ndim == 1: + return self*x + elif x.ndim == 2 and x.shape[1] == 1: + if x.flags.c_contiguous: + return self*x[:, 0] + else: + y = np.zeros((x.shape[0]), dtype=x.dtype) + y[:] = x[:, 0] + return self*y + else: + raise NotImplementedError() + + from scipy.sparse.linalg import LinearOperator as ScipyLinearOperator + return ScipyLinearOperator(shape=self.shape, matvec=matvec) + + # def getDenseOpFromApply(self): + # cdef: + # INDEX_t i + # {SCALAR}_t[::1] x = np.zeros((self.shape[1]), dtype={SCALAR}) + # {SCALAR}_t[::1, :] B = np.zeros(self.shape, dtype={SCALAR}, order='F') + # for i in range(self.shape[1]): + # if i > 0: + # x[i-1] = 0. + # x[i] = 1. + # self.matvec(x, B[:, i]) + # return np.ascontiguousarray(B) + + def __getstate__(self): + return + + def __setstate__(self, state): + pass + + def __repr__(self): + return '<%dx%dx%d %s>' % (self.num_rows, + self.num_columns, + self.vectorSize, + self.__class__.__name__) + + cdef void setEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val): + raise NotImplementedError() + + def setEntry_py(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val): + self.setEntry(I, J, val) + + cdef void addToEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val): + raise NotImplementedError() + + def addToEntry_py(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val): + self.addToEntry(I, J, val) + + cdef void getEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val): + raise NotImplementedError() + + def getEntry_py(self, INDEX_t I, INDEX_t J, {SCALAR}_t[::1] val): + return self.getEntry(I, J, val) diff --git a/base/PyNucleus_base/SSS_LinearOperator_decl_{SCALAR}.pxi b/base/PyNucleus_base/SSS_LinearOperator_decl_{SCALAR}.pxi index f8a306f..dcb799d 100644 --- a/base/PyNucleus_base/SSS_LinearOperator_decl_{SCALAR}.pxi +++ b/base/PyNucleus_base/SSS_LinearOperator_decl_{SCALAR}.pxi @@ -22,3 +22,5 @@ cdef class {SCALAR_label}SSS_LinearOperator({SCALAR_label}LinearOperator): {SCALAR}_t[::1] y) except -1 cdef void setEntry(self, INDEX_t I, INDEX_t J, {SCALAR}_t val) cdef {SCALAR}_t getEntry({SCALAR_label}SSS_LinearOperator self, INDEX_t I, INDEX_t J) + + diff --git a/base/PyNucleus_base/SSS_LinearOperator_{SCALAR}.pxi b/base/PyNucleus_base/SSS_LinearOperator_{SCALAR}.pxi index a049923..7c72a68 100644 --- a/base/PyNucleus_base/SSS_LinearOperator_{SCALAR}.pxi +++ b/base/PyNucleus_base/SSS_LinearOperator_{SCALAR}.pxi @@ -326,3 +326,5 @@ cdef class {SCALAR_label}SSS_LinearOperator({SCALAR_label}LinearOperator): self.diagonal[i] *= scaling for i in range(self.data.shape[0]): self.data[i] *= scaling + + diff --git a/base/PyNucleus_base/utilsFem.py b/base/PyNucleus_base/utilsFem.py index 39c8ba3..6b78a1a 100644 --- a/base/PyNucleus_base/utilsFem.py +++ b/base/PyNucleus_base/utilsFem.py @@ -1479,6 +1479,7 @@ def __call__(self): if (newValue != oldValue).any(): needToBuild = True elif newValue != oldValue: + dependencyLogger.log(self.logLevel, 'Values for {} differ: \'{}\' != \'{}\', calling \'{}\''.format(prop, oldValue, newValue, self.fun.__name__)) needToBuild = True except Exception as e: dependencyLogger.log(logging.WARN, 'Cannot compare values {}, {} for property \'{}\', exception {}, force call \'{}\''.format(oldValue, newValue, prop, e, self.fun.__name__)) diff --git a/drivers/runNonlocal.py b/drivers/runNonlocal.py index 29a620c..d585dbf 100755 --- a/drivers/runNonlocal.py +++ b/drivers/runNonlocal.py @@ -6,6 +6,7 @@ # If you want to use this code, please refer to the README.rst and LICENSE files. # ################################################################################### +from mpi4py import MPI import numpy as np from PyNucleus import driver, DIRICHLET, NEUMANN from PyNucleus.nl import (nonlocalPoissonProblem, @@ -15,7 +16,7 @@ description = """Solves a nonlocal Poisson problem with finite horizon.""" -d = driver(description=description) +d = driver(MPI.COMM_WORLD, description=description) p = nonlocalPoissonProblem(d) discrProblem = discretizedNonlocalProblem(d, p) diff --git a/drivers/runParallelGMG.py b/drivers/runParallelGMG.py index 65b4b54..43998fd 100755 --- a/drivers/runParallelGMG.py +++ b/drivers/runParallelGMG.py @@ -174,7 +174,7 @@ A.residual_py(x, rhs, r) resNorm = r.norm(False) rate.add('Rate of convergence P'+label, (resNorm/r0)**(1/numIter), tested=False if label == 'BICGSTAB' else None) - its.add('Number of iterations P'+label, numIter, aTol=1 if label == 'BICGSTAB' else None) + its.add('Number of iterations P'+label, numIter, aTol=2 if label == 'BICGSTAB' else None) res.add('Residual norm P'+label, resNorm) resHist.add('P'+label, residuals, tested=False if label == 'BICGSTAB' else None) diff --git a/fem/PyNucleus_fem/DoFMaps.pyx b/fem/PyNucleus_fem/DoFMaps.pyx index ce4c0b3..89da5b2 100644 --- a/fem/PyNucleus_fem/DoFMaps.pyx +++ b/fem/PyNucleus_fem/DoFMaps.pyx @@ -2295,39 +2295,6 @@ def generateLocalMassMatrix(DoFMap dm, DoFMap dm2=None): return generic_matrix(entries) -cdef class lookupFunction(function): - cdef: - meshBase mesh - public DoFMap dm - public REAL_t[::1] u - public cellFinder2 cellFinder - - def __init__(self, meshBase mesh, DoFMap dm, REAL_t[::1] u, cellFinder2 cF=None): - self.mesh = mesh - self.dm = dm - self.u = u - if cF is None: - self.cellFinder = cellFinder2(self.mesh) - else: - self.cellFinder = cF - - cdef REAL_t eval(self, REAL_t[::1] x): - cdef: - shapeFunction shapeFun - REAL_t val - INDEX_t cellNo, dof, k - cellNo = self.cellFinder.findCell(x) - if cellNo == -1: - return 0. - val = 0. - for k in range(self.dm.dofs_per_element): - dof = self.dm.cell2dof(cellNo, k) - if dof >= 0: - shapeFun = self.dm.localShapeFunctions[k] - val += shapeFun.eval(self.cellFinder.bary)*self.u[dof] - return val - - cdef class elementSizeFunction(function): cdef: meshBase mesh diff --git a/fem/PyNucleus_fem/factories.py b/fem/PyNucleus_fem/factories.py index 2ecce5e..9fe4114 100644 --- a/fem/PyNucleus_fem/factories.py +++ b/fem/PyNucleus_fem/factories.py @@ -34,7 +34,7 @@ radialIndicator, fractalDiffusivity, expDiffusivity, componentVectorFunction) -from . DoFMaps import lookupFunction +from . lookupFunction import lookupFunction rhsFunSin1D = _rhsFunSin1D() diff --git a/fem/PyNucleus_fem/functions.pyx b/fem/PyNucleus_fem/functions.pyx index 264b0cf..52a93e4 100644 --- a/fem/PyNucleus_fem/functions.pyx +++ b/fem/PyNucleus_fem/functions.pyx @@ -173,6 +173,16 @@ cdef class constant(function): def __repr__(self): return '{}'.format(self.value) + def __eq__(self, other): + cdef: + constant o + if isinstance(other, constant): + o = other + return self.value == o.value + else: + return False + + cdef class monomial(function): cdef: @@ -1284,6 +1294,24 @@ cdef class radialIndicator(function): r += (x[i]-self.center[i])*(x[i]-self.center[i]) return r <= self.radius + def __eq__(self, other): + cdef: + radialIndicator o + INDEX_t i + if isinstance(other, radialIndicator): + o = other + if self.radius != o.radius: + return False + if self.centerIsOrigin != o.centerIsOrigin: + return False + for i in range(self.center.shape[0]): + if not self.center[i] == o.center[i]: + return False + return True + else: + return False + + cdef class squareIndicator(function): cdef REAL_t[::1] a, b @@ -1300,6 +1328,21 @@ cdef class squareIndicator(function): return False return True + def __eq__(self, other): + cdef: + squareIndicator o + INDEX_t i + if isinstance(other, squareIndicator): + o = other + for i in range(self.a.shape[0]): + if self.a[i] != o.a[i]: + return False + if self.b[i] != o.b[i]: + return False + return True + else: + return False + cdef class proj(function): cdef: diff --git a/fem/PyNucleus_fem/lookupFunction.pxd b/fem/PyNucleus_fem/lookupFunction.pxd new file mode 100644 index 0000000..0c3c996 --- /dev/null +++ b/fem/PyNucleus_fem/lookupFunction.pxd @@ -0,0 +1,19 @@ +################################################################################### +# Copyright 2021 National Technology & Engineering Solutions of Sandia, # +# LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the # +# U.S. Government retains certain rights in this software. # +# If you want to use this code, please refer to the README.rst and LICENSE files. # +################################################################################### + +from PyNucleus_base.myTypes cimport REAL_t, INDEX_t +from . functions cimport function +from . meshCy cimport meshBase, cellFinder2 +from . DoFMaps cimport DoFMap + + +cdef class lookupFunction(function): + cdef: + meshBase mesh + public DoFMap dm + public REAL_t[::1] u + public cellFinder2 cellFinder diff --git a/fem/PyNucleus_fem/lookupFunction.pyx b/fem/PyNucleus_fem/lookupFunction.pyx new file mode 100644 index 0000000..2c8b15d --- /dev/null +++ b/fem/PyNucleus_fem/lookupFunction.pyx @@ -0,0 +1,35 @@ +################################################################################### +# Copyright 2021 National Technology & Engineering Solutions of Sandia, # +# LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the # +# U.S. Government retains certain rights in this software. # +# If you want to use this code, please refer to the README.rst and LICENSE files. # +################################################################################### + +from . DoFMaps cimport shapeFunction + + +cdef class lookupFunction(function): + def __init__(self, meshBase mesh, DoFMap dm, REAL_t[::1] u, cellFinder2 cF=None): + self.mesh = mesh + self.dm = dm + self.u = u + if cF is None: + self.cellFinder = cellFinder2(self.mesh) + else: + self.cellFinder = cF + + cdef REAL_t eval(self, REAL_t[::1] x): + cdef: + shapeFunction shapeFun + REAL_t val + INDEX_t cellNo, dof, k + cellNo = self.cellFinder.findCell(x) + if cellNo == -1: + return 0. + val = 0. + for k in range(self.dm.dofs_per_element): + dof = self.dm.cell2dof(cellNo, k) + if dof >= 0: + shapeFun = self.dm.localShapeFunctions[k] + val += shapeFun.eval(self.cellFinder.bary)*self.u[dof] + return val diff --git a/fem/PyNucleus_fem/meshCy.pxd b/fem/PyNucleus_fem/meshCy.pxd index 020647e..7c5fd17 100644 --- a/fem/PyNucleus_fem/meshCy.pxd +++ b/fem/PyNucleus_fem/meshCy.pxd @@ -48,6 +48,11 @@ cdef class meshBase: REAL_t[:, ::1] simplexMem, REAL_t[::1] baryMem, REAL_t tol=*) + cdef BOOL_t vertexInCellPtr(self, REAL_t* vertex, + INDEX_t cellNo, + REAL_t[:, ::1] simplexMem, + REAL_t[::1] baryMem, + REAL_t tol=*) cdef void decode_edge(ENCODE_t encodeVal, INDEX_t[::1] e) @@ -119,6 +124,7 @@ cdef class cellFinder2: INDEX_t[::1] key intTuple myKey cdef INDEX_t findCell(self, REAL_t[::1] vertex) + cdef INDEX_t findCellPtr(self, REAL_t* vertex) cdef void getBarycentricCoords1D(REAL_t[:, ::1] simplex, REAL_t[::1] x, REAL_t[::1] bary) cdef void getBarycentricCoords2D(REAL_t[:, ::1] simplex, REAL_t[::1] x, REAL_t[::1] bary) diff --git a/fem/PyNucleus_fem/meshCy.pyx b/fem/PyNucleus_fem/meshCy.pyx index 03ca7bb..e9cef44 100644 --- a/fem/PyNucleus_fem/meshCy.pyx +++ b/fem/PyNucleus_fem/meshCy.pyx @@ -680,6 +680,21 @@ cdef class meshBase: return False return True + cdef BOOL_t vertexInCellPtr(self, REAL_t* vertex, INDEX_t cellNo, REAL_t[:, ::1] simplexMem, REAL_t[::1] baryMem, REAL_t tol=0.): + cdef: + INDEX_t i + self.getSimplex(cellNo, simplexMem) + if self.dim == 1: + getBarycentricCoords1DPtr(simplexMem, vertex, baryMem) + elif self.dim == 2: + getBarycentricCoords2DPtr(simplexMem, vertex, baryMem) + else: + raise NotImplementedError() + for i in range(self.dim+1): + if baryMem[i] < -tol: + return False + return True + def vertexInCell_py(self, REAL_t[::1] vertex, INDEX_t cellNo, REAL_t tol=0.): simplex = uninitialized((self.dim+1, self.dim), dtype=REAL) bary = uninitialized((self.dim+1), dtype=REAL) @@ -2121,6 +2136,51 @@ cdef class cellFinder2: return cellNo return -1 + cdef INDEX_t findCellPtr(self, REAL_t* vertex): + cdef: + INDEX_t j, cellNo, vertexNo, v + set candidates, toCheck = set() + productIterator pit + INDEX_t[::1] keyCenter + for j in range(self.mesh.dim): + self.key[j] = ((vertex[j]-self.x_min[j]) * self.diamInv[j]) + try: + candidates = self.lookup[self.myKey] + except KeyError: + keyCenter = np.array(self.key, copy=True) + pit = productIterator(3, self.mesh.dim) + candidates = set() + pit.reset() + while pit.step(): + for j in range(self.mesh.dim): + self.key[j] = keyCenter[j] + pit.idx[j]-1 + try: + candidates |= self.lookup[self.myKey] + except KeyError: + pass + + # check if the vertex is in any of the cells + for cellNo in candidates: + if self.mesh.vertexInCellPtr(vertex, cellNo, self.simplex, self.bary): + return cellNo + # add neighboring cells of candidate cells + for cellNo in candidates: + for vertexNo in range(self.mesh.dim+1): + v = self.mesh.cells[cellNo, vertexNo] + toCheck |= self.v2c[v] + toCheck -= candidates + for cellNo in toCheck: + if self.mesh.vertexInCellPtr(vertex, cellNo, self.simplex, self.bary): + return cellNo + # allow for some extra room + for cellNo in candidates: + if self.mesh.vertexInCellPtr(vertex, cellNo, self.simplex, self.bary, 1e-15): + return cellNo + for cellNo in toCheck: + if self.mesh.vertexInCellPtr(vertex, cellNo, self.simplex, self.bary, 1e-15): + return cellNo + return -1 + def findCell_py(self, REAL_t[::1] vertex): return self.findCell(vertex) @@ -2145,6 +2205,26 @@ cdef void getBarycentricCoords2D(REAL_t[:, ::1] simplex, REAL_t[::1] x, REAL_t[: bary[2] = 1. - bary[0] - bary[1] +cdef void getBarycentricCoords1DPtr(REAL_t[:, ::1] simplex, REAL_t* x, REAL_t[::1] bary): + cdef: + REAL_t vol + vol = simplex[0, 0]-simplex[1, 0] + bary[0] = (x[0]-simplex[1, 0])/vol + bary[1] = 1.-bary[0] + + +cdef void getBarycentricCoords2DPtr(REAL_t[:, ::1] simplex, REAL_t* x, REAL_t[::1] bary): + cdef: + REAL_t vol + vol = ((simplex[0, 0]-simplex[1,0])*(simplex[2, 1]-simplex[1,1]) - + (simplex[0, 1]-simplex[1,1])*(simplex[2, 0]-simplex[1,0])) + bary[0] = ((x[0]-simplex[1, 0])*(simplex[2, 1]-simplex[1, 1]) - + (x[1]-simplex[1, 1])*(simplex[2, 0]-simplex[1, 0]))/vol + bary[1] = ((x[0]-simplex[2, 0])*(simplex[0, 1]-simplex[2, 1]) - + (x[1]-simplex[2, 1])*(simplex[0, 0]-simplex[2, 0]))/vol + bary[2] = 1. - bary[0] - bary[1] + + def getSubmesh2(mesh, INDEX_t[::1] newCellIndices, INDEX_t num_cells=-1): cdef: INDEX_t[:, ::1] cells, newCells diff --git a/fem/setup.py b/fem/setup.py index 3706859..de5f6aa 100644 --- a/fem/setup.py +++ b/fem/setup.py @@ -55,6 +55,8 @@ sources=[p.folder+"repartitioner.pyx"]) p.addExtension("DoFMaps", sources=[p.folder+"DoFMaps.pyx"]) +p.addExtension("lookupFunction", + sources=[p.folder+"lookupFunction.pyx"]) p.addExtension("quadrature", sources=[p.folder+"quadrature.pyx"]) p.addExtension("meshOverlaps", diff --git a/multilevelSolver/PyNucleus_multilevelSolver/hierarchies.py b/multilevelSolver/PyNucleus_multilevelSolver/hierarchies.py index 9a393a3..1abb760 100644 --- a/multilevelSolver/PyNucleus_multilevelSolver/hierarchies.py +++ b/multilevelSolver/PyNucleus_multilevelSolver/hierarchies.py @@ -78,6 +78,8 @@ def __init__(self, meshLevel, params, comm=None, self.buildType = [NO_BUILD]*(self.params['noRef']+1) elif self.params['assemble'] == 'restrictionProlongation': self.buildType = [RESTRICTION_PROLONGATION_ONLY]*(self.params['noRef']+1) + elif self.params['assemble'] == 'dofmap only last': + self.buildType = [NO_BUILD]*self.params['noRef']+[DOFMAPS] else: raise NotImplementedError() diff --git a/nl/PyNucleus_nl/clusterMethodCy.pyx b/nl/PyNucleus_nl/clusterMethodCy.pyx index 071c5f4..752c34e 100644 --- a/nl/PyNucleus_nl/clusterMethodCy.pyx +++ b/nl/PyNucleus_nl/clusterMethodCy.pyx @@ -1911,7 +1911,7 @@ cdef class productIterator: return True -def assembleFarFieldInteractions(FractionalKernel kernel, dict Pfar, INDEX_t m, DoFMap dm): +def assembleFarFieldInteractions(Kernel kernel, dict Pfar, INDEX_t m, DoFMap dm): cdef: INDEX_t lvl REAL_t[:, ::1] box1, box2, x, y diff --git a/nl/PyNucleus_nl/discretizedProblems.py b/nl/PyNucleus_nl/discretizedProblems.py index 060e335..47b6467 100644 --- a/nl/PyNucleus_nl/discretizedProblems.py +++ b/nl/PyNucleus_nl/discretizedProblems.py @@ -21,8 +21,7 @@ from PyNucleus_fem.DoFMaps import Product_DoFMap from PyNucleus_multilevelSolver import hierarchyManager from copy import copy -from . helpers import (multilevelDirichletCondition, - paramsForFractionalHierarchy) +from . helpers import paramsForFractionalHierarchy from . fractionalOrders import singleVariableUnsymmetricFractionalOrder from . kernelsCy import FRACTIONAL from . nonlocalProblems import (DIRICHLET, @@ -352,18 +351,49 @@ def setDriverArgs(self): self.setDriverFlag('matrixFormat', acceptedValues=['H2', 'sparse', 'sparsified', 'dense'], help='matrix format') self.setDriverFlag('debugAssemblyTimes', False) - def processCmdline(self, params): - matrixFormat = params['matrixFormat'] - kernelType = params['kernelType'] - if matrixFormat.upper() == 'H2' and kernelType != 'fractional': - matrixFormat = 'sparse' - params['matrixFormat'] = matrixFormat - super().processCmdline(params) + @generates(['meshHierarchy', 'finalMesh', + 'dm', 'dmBC', 'dmInterior', + 'R_interior', 'P_interior', + 'R_bc', 'P_bc']) + def buildMeshHierarchy(self, mesh, solverType, domainIndicator, fluxIndicator, noRef, element): + with self.timer('hierarchy - meshes'): + params = {} + params['domain'] = mesh + params['solver'] = solverType + params['tag'] = domainIndicator+fluxIndicator + params['element'] = element + params['keepMeshes'] = 'all' + params['keepAllDoFMaps'] = True + params['buildMass'] = True + params['assemble'] = 'restrictionProlongation' if solverType.find('mg') >= 0 else 'dofmap only last' + params['logging'] = True + if self.debugAssemblyTimes: + from PyNucleus.base.utilsFem import TimerManager + tm = TimerManager(self.driver.logger, comm=self.driver.comm, memoryProfiling=False, loggingSubTimers=True) + params['PLogger'] = tm.PLogger - @generates(['hierarchy', 'bc', 'finalMesh', 'dm', 'dmBC', 'dmInterior', 'R_interior', 'P_interior', 'R_BC', 'P_BC']) - def buildHierarchy(self, mesh, kernel, sArgs, rangedKernel, solverType, matrixFormat, tag, boundaryCondition, domainIndicator, fluxIndicator, + comm = self.driver.comm + onRanks = [0] + if comm is not None and comm.size > 1: + onRanks = range(comm.size) + hierarchies, connectors = paramsForFractionalHierarchy(noRef, params, onRanks) + hM = hierarchyManager(hierarchies, connectors, params, comm) + hM.setup() + self.meshHierarchy = hM + self.finalMesh = hM['fine'].meshLevels[-1].mesh + + self.dmInterior = hM['fine'].algebraicLevels[-1].DoFMap + self.dmBC = self.dmInterior.getComplementDoFMap() + self.dm, self.R_interior, self.R_bc = self.dmInterior.getFullDoFMap(self.dmBC) + self.P_interior = self.R_interior.transpose() + self.P_bc = self.R_bc.transpose() + + @generates('hierarchy') + def buildHierarchy(self, + meshHierarchy, + dm, dmBC, dmInterior, + kernel, sArgs, rangedKernel, solverType, matrixFormat, tag, boundaryCondition, domainIndicator, fluxIndicator, zeroExterior, noRef, eta, target_order, element, quadType, quadTypeBoundary): - assert matrixFormat != 'H2' or (kernel.kernelType == FRACTIONAL), 'Hierarchical matrices are only implemented for fractional kernels' if rangedKernel is not None: hierarchy = self.directlyGetWithoutChecks('hierarchy') if hierarchy is not None: @@ -381,78 +411,84 @@ def buildHierarchy(self, mesh, kernel, sArgs, rangedKernel, solverType, matrixFo assert hierarchy != newHierarchy self.hierarchy = newHierarchy - for prop in ['bc', 'finalMesh', 'dm', 'dmBC', 'dmInterior', 'R_interior', 'P_interior', 'R_BC', 'P_BC']: - setattr(self, prop, self.directlyGetWithoutChecks(prop)) + # for prop in ['bc', 'finalMesh', 'dm', 'dmBC', 'dmInterior', 'R_interior', 'P_interior', 'R_bc', 'P_bc']: + # setattr(self, prop, self.directlyGetWithoutChecks(prop)) return - with self.timer('hierarchy'): - params = {} - params['domain'] = mesh - if rangedKernel is None: - params['kernel'] = kernel - else: - params['kernel'] = rangedKernel - params['solver'] = solverType - params['tag'] = tag - params['element'] = element - params['boundaryCondition'] = boundaryCondition - params['zeroExterior'] = zeroExterior - params['target_order'] = target_order - params['eta'] = eta - params['keepMeshes'] = 'all' - params['keepAllDoFMaps'] = True - params['buildMass'] = True - params['assemble'] = 'ALL' if solverType.find('mg') >= 0 else 'last' - params['dense'] = matrixFormat == 'dense' - params['matrixFormat'] = matrixFormat - params['logging'] = True - if self.debugAssemblyTimes: - from PyNucleus.base.utilsFem import TimerManager - tm = TimerManager(self.driver.logger, comm=self.driver.comm, memoryProfiling=False, loggingSubTimers=True) - params['PLogger'] = tm.PLogger + hM = meshHierarchy - if quadType == 'auto': - quadType = 'classical-refactored' - params['quadType'] = quadType + assemblyParams = {} + if quadType == 'auto': + quadType = 'classical-refactored' + assemblyParams['quadType'] = quadType - if quadTypeBoundary == 'auto': - quadTypeBoundary = 'classical-refactored' - params['quadTypeBoundary'] = quadTypeBoundary + if quadTypeBoundary == 'auto': + quadTypeBoundary = 'classical-refactored' + assemblyParams['quadTypeBoundary'] = quadTypeBoundary + if rangedKernel is None: + assemblyParams['kernel'] = kernel + else: + assemblyParams['kernel'] = rangedKernel + assemblyParams['boundaryCondition'] = boundaryCondition + assemblyParams['zeroExterior'] = zeroExterior + assemblyParams['target_order'] = target_order + assemblyParams['eta'] = eta + assemblyParams['dense'] = matrixFormat == 'dense' + assemblyParams['matrixFormat'] = matrixFormat + + with self.timer('hierarchy - matrices'): + from PyNucleus_multilevelSolver.levels import ASSEMBLY + if solverType.find('mg') >= 0: + for subHierarchy in hM.builtHierarchies: + for level in subHierarchy.algebraicLevels: + level.params.update(assemblyParams) + level.build(ASSEMBLY) + else: + level = hM.builtHierarchies[-1].algebraicLevels[-1] + level.params.update(assemblyParams) + level.build(ASSEMBLY) - comm = self.driver.comm - onRanks = [0] - if comm is not None and comm.size > 1: - onRanks = range(comm.size) - hierarchies, connectors = paramsForFractionalHierarchy(noRef, params, onRanks) - hM = hierarchyManager(hierarchies, connectors, params, comm) - hM.setup() + hierarchy = hM.getLevelList() - if boundaryCondition == DIRICHLET: - bc = multilevelDirichletCondition(hM.getLevelList(), domainIndicator, fluxIndicator) - hierarchy = bc.naturalLevels - dmInterior = bc.naturalDoFMap - dmBC = bc.dirichletDoFMap - else: - hierarchy = hM.getLevelList() - bc = None - dmInterior = hierarchy[-1]['DoFMap'] - dmBC = dmInterior.getComplementDoFMap() if rangedKernel is not None: hierarchy[0]['sArgs'] = sArgs s = kernel.sValue for lvl in range(len(hierarchy)): if 'A' in hierarchy[lvl]: hierarchy[lvl]['A'].set(s, 0) - self.dmBC = dmBC - self.dm, self.R_interior, self.R_bc = dmInterior.getFullDoFMap(dmBC) - self.P_interior = self.R_interior.transpose() - self.P_bc = self.R_bc.transpose() - self.finalMesh = hM['fine'].meshLevels[-1].mesh - self.bc = bc - self.dmInterior = dmInterior self.hierarchy = hierarchy - assert 2*self.finalMesh.h < kernel.horizon.value, "h = {}, horizon = {}".format(self.finalMesh.h, kernel.horizon.value) + if kernel is not None: + assert 2*self.finalMesh.h < kernel.horizon.value, "h = {}, horizon = {}".format(self.finalMesh.h, kernel.horizon.value) + + @generates('A_BC') + def buildBCoperator(self, dmInterior, dmBC, + kernel, sArgs, rangedKernel, solverType, matrixFormat, tag, boundaryCondition, + zeroExterior, noRef, eta, target_order, element, quadType, quadTypeBoundary): + if boundaryCondition == DIRICHLET: + from . helpers import getFracLapl + assemblyParams = {} + if quadType == 'auto': + quadType = 'classical-refactored' + assemblyParams['quadType'] = quadType + + if quadTypeBoundary == 'auto': + quadTypeBoundary = 'classical-refactored' + assemblyParams['quadTypeBoundary'] = quadTypeBoundary + if rangedKernel is None: + assemblyParams['kernel'] = kernel + else: + assemblyParams['kernel'] = rangedKernel + assemblyParams['boundaryCondition'] = boundaryCondition + assemblyParams['zeroExterior'] = zeroExterior + assemblyParams['target_order'] = target_order + assemblyParams['eta'] = eta + assemblyParams['dense'] = matrixFormat == 'dense' + assemblyParams['matrixFormat'] = matrixFormat + assemblyParams['tag'] = tag + self.A_BC = getFracLapl(dmInterior.mesh, dmInterior, dm2=dmBC, **assemblyParams) + else: + self.A_BC = None @generates('mass') def buildMass(self, dm): @@ -466,16 +502,23 @@ def buildMassInterior(self, dmInterior): def getOperators(self, hierarchy): self.A = hierarchy[-1]['A'] + @generates('A_derivative') + def getDerivativeOperator(self, kernel, dmInterior, matrixFormat, eta, target_order): + self.A_derivative = dmInterior.assembleNonlocal(kernel.getDerivativeKernel(derivative=1), matrixFormat=matrixFormat, params={'eta': eta, + 'target_order': target_order}) + @generates('b') - def buildRHS(self, rhs, dim, bc, dirichletData, boundaryCondition, solverType, dmInterior, hierarchy): + def buildRHS(self, rhs, dim, A_BC, dmBC, dirichletData, boundaryCondition, solverType, dmInterior, hierarchy): self.b = dmInterior.assembleRHS(rhs, qr=simplexXiaoGimbutas(3, dim)) - if bc is not None: - bc.setDirichletData(dirichletData) - bc.applyRHScorrection(self.b) + + if A_BC is not None: + assert dmInterior.num_dofs == A_BC.num_rows + assert dmBC.num_dofs == A_BC.num_columns + if dmBC.num_dofs > 0: + self.b -= A_BC*dmBC.interpolate(dirichletData) # pure Neumann condition -> project out nullspace if boundaryCondition in (NEUMANN, HOMOGENEOUS_NEUMANN): - assert bc is None if solverType.find('mg') >= 0: hierarchy[0]['A'] = hierarchy[0]['A'] + Dense_LinearOperator.ones(*hierarchy[0]['A'].shape) else: @@ -500,7 +543,7 @@ def buildSolver(self, solverType, tol, maxiter, hierarchy, kernel): self.solver = solver @generates('modelSolution') - def solve(self, b, bc, dm, dmInterior, P_interior, solver, boundaryCondition, analyticSolution, dirichletData, tol, maxiter): + def solve(self, b, dm, dmInterior, dmBC, P_interior, P_bc, solver, boundaryCondition, analyticSolution, dirichletData, tol, maxiter): uInterior = dmInterior.zeros() with self.timer('solve {}'.format(self.__class__.__name__)): its = solver(b, uInterior) @@ -508,8 +551,8 @@ def solve(self, b, bc, dm, dmInterior, P_interior, solver, boundaryCondition, an resError = (b-solver.A*uInterior).norm(False) if isinstance(solver, iterative_solver): - assert its < maxiter, "Only reached residual error {} > tol = {} in {} iterations".format(resError, tol, its) - # assert resError < tol, "Only reached residual error {} > tol = {}".format(resError, tol) + if its >= maxiter-1: + self.driver.logger.warn("WARNING: Only reached residual error {} > tol = {} in {} iterations".format(resError, tol, its)) # pure Neumann condition -> add nullspace components to match analytic solution if boundaryCondition in (NEUMANN, HOMOGENEOUS_NEUMANN) and analyticSolution is not None: @@ -519,7 +562,7 @@ def solve(self, b, bc, dm, dmInterior, P_interior, solver, boundaryCondition, an u = dm.empty() if boundaryCondition in (DIRICHLET, ): - u.assign(bc.augmentDirichlet(uInterior)) + u.assign(P_interior*uInterior+P_bc*dmBC.interpolate(dirichletData)) else: u.assign(P_interior*uInterior) diff --git a/nl/PyNucleus_nl/fractionalLaplacian1D.pyx b/nl/PyNucleus_nl/fractionalLaplacian1D.pyx index 098a2cb..61a7ed0 100644 --- a/nl/PyNucleus_nl/fractionalLaplacian1D.pyx +++ b/nl/PyNucleus_nl/fractionalLaplacian1D.pyx @@ -30,7 +30,7 @@ ALL.set() cdef class fractionalLaplacian1DZeroExterior(nonlocalLaplacian1D): - def __init__(self, FractionalKernel kernel, meshBase mesh, DoFMap DoFMap, num_dofs=None, **kwargs): + def __init__(self, Kernel kernel, meshBase mesh, DoFMap DoFMap, num_dofs=None, **kwargs): manifold_dim2 = mesh.dim-1 super(fractionalLaplacian1DZeroExterior, self).__init__(kernel, mesh, DoFMap, num_dofs, manifold_dim2=manifold_dim2, **kwargs) self.symmetricCells = False @@ -394,6 +394,8 @@ cdef class fractionalLaplacian1D(nonlocalLaplacian1D): contrib[k] = val*vol + + cdef class fractionalLaplacian1D_nonsym(fractionalLaplacian1D): """The local stiffness matrix @@ -566,9 +568,10 @@ cdef class fractionalLaplacian1D_nonsym(fractionalLaplacian1D): contrib[k] = val*vol + cdef class fractionalLaplacian1D_boundary(fractionalLaplacian1DZeroExterior): def __init__(self, - FractionalKernel kernel, + Kernel kernel, meshBase mesh, DoFMap DoFMap, quad_order_diagonal=None, @@ -577,10 +580,8 @@ cdef class fractionalLaplacian1D_boundary(fractionalLaplacian1DZeroExterior): **kwargs): super(fractionalLaplacian1D_boundary, self).__init__(kernel, mesh, DoFMap, num_dofs, **kwargs) - if self.kernel.variableOrder: - smin, smax = self.kernel.s.min, self.kernel.s.max - else: - smin, smax = self.kernel.sValue, self.kernel.sValue + smin = max(0.5*(-self.kernel.min_singularity-1.), 0.) + smax = max(0.5*(-self.kernel.max_singularity-1.), 0.) if target_order is None: # this is the desired local quadrature error target_order = self.DoFMap.polynomialOrder+1-smin @@ -591,7 +592,7 @@ cdef class fractionalLaplacian1D_boundary(fractionalLaplacian1DZeroExterior): quad_order_diagonal = max(np.ceil(((target_order+1.)*log(self.num_dofs*self.H0)+(2.*smax-1.)*abs(log(self.hmin/self.H0)))/0.8), 2) self.quad_order_diagonal = quad_order_diagonal - if not self.kernel.variableOrder: + if (self.kernel.kernelType != FRACTIONAL) or (not self.kernel.variableOrder): self.getNearQuadRule(COMMON_VERTEX) cdef panelType getQuadOrder(self, @@ -601,7 +602,7 @@ cdef class fractionalLaplacian1D_boundary(fractionalLaplacian1DZeroExterior): cdef: panelType panel, panel2 REAL_t logdh1 = max(log(d/h1), 0.), logdh2 = max(log(d/h2), 0.) - REAL_t s = self.kernel.sValue + REAL_t s = max(0.5*(-self.kernel.getSingularityValue()-1.), 0.) REAL_t h panel = max(ceil(((self.target_order+1.)*log(self.num_dofs*self.H0) + (2.*s-1.)*abs(log(h2/self.H0)) - 2.*s*log(d/h2)) / (logdh1 + 0.8)), @@ -711,6 +712,8 @@ cdef class fractionalLaplacian1D_boundary(fractionalLaplacian1DZeroExterior): y[j] = simplex2[self.perm2[0], j]*qr.nodes[2, m] self.temp[m] = qr.weights[m] * self.kernel.evalPtr(dim, &x[0], &y[0]) + contrib[:] = 0. + for I in range(dofs_per_element): i = self.perm[I] for J in range(I, dofs_per_element): @@ -728,3 +731,4 @@ cdef class fractionalLaplacian1D_boundary(fractionalLaplacian1DZeroExterior): raise NotImplementedError('Unknown panel type: {}'.format(panel)) + diff --git a/nl/PyNucleus_nl/fractionalLaplacian2D.pyx b/nl/PyNucleus_nl/fractionalLaplacian2D.pyx index cb442d6..f47ca73 100644 --- a/nl/PyNucleus_nl/fractionalLaplacian2D.pyx +++ b/nl/PyNucleus_nl/fractionalLaplacian2D.pyx @@ -30,7 +30,7 @@ ALL.set() cdef class fractionalLaplacian2DZeroExterior(nonlocalLaplacian2D): - def __init__(self, FractionalKernel kernel, meshBase mesh, DoFMap DoFMap, num_dofs=None, **kwargs): + def __init__(self, Kernel kernel, meshBase mesh, DoFMap DoFMap, num_dofs=None, **kwargs): manifold_dim2 = mesh.dim-1 super(fractionalLaplacian2DZeroExterior, self).__init__(kernel, mesh, DoFMap, num_dofs, manifold_dim2=manifold_dim2, **kwargs) self.symmetricCells = False @@ -1137,7 +1137,7 @@ cdef class fractionalLaplacian2D_boundary(fractionalLaplacian2DZeroExterior): """ def __init__(self, - FractionalKernel kernel, + Kernel kernel, meshBase mesh, DoFMap DoFMap, target_order=None, @@ -1146,7 +1146,7 @@ cdef class fractionalLaplacian2D_boundary(fractionalLaplacian2DZeroExterior): **kwargs): super(fractionalLaplacian2D_boundary, self).__init__(kernel, mesh, DoFMap, num_dofs, **kwargs) - smax = self.kernel.s.max + smax = max(0.5*(-self.kernel.max_singularity-2.), 0.) if target_order is None: # this is the desired global order wrt to the number of DoFs # target_order = (2.-s)/self.dim @@ -1158,7 +1158,7 @@ cdef class fractionalLaplacian2D_boundary(fractionalLaplacian2DZeroExterior): quad_order_diagonal = max(np.ceil((target_order+0.5+smax)/(0.35)*abs(np.log(self.hmin/self.H0))), 2) self.quad_order_diagonal = quad_order_diagonal - if not self.kernel.variableOrder: + if (self.kernel.kernelType != FRACTIONAL) or (not self.kernel.variableOrder): self.getNearQuadRule(COMMON_EDGE) self.getNearQuadRule(COMMON_VERTEX) @@ -1171,7 +1171,7 @@ cdef class fractionalLaplacian2D_boundary(fractionalLaplacian2DZeroExterior): REAL_t logdh1 = max(log(d/h1), 0.), logdh2 = max(log(d/h2), 0.) REAL_t logh1H0 = abs(log(h1/self.H0)), logh2H0 = abs(log(h2/self.H0)) REAL_t loghminH0 = max(logh1H0, logh2H0) - REAL_t s = self.kernel.sValue + REAL_t s = max(0.5*(-self.kernel.getSingularityValue()-2.), 0.) REAL_t h panel = max(ceil(((0.5*self.target_order+0.25)*log(self.num_dofs*self.H0**2) + loghminH0 + (s-1.)*logh2H0 - s*logdh2) / (max(logdh1, 0) + 0.35)), diff --git a/nl/PyNucleus_nl/fractionalOrders.pxd b/nl/PyNucleus_nl/fractionalOrders.pxd index 645c43d..5d94075 100644 --- a/nl/PyNucleus_nl/fractionalOrders.pxd +++ b/nl/PyNucleus_nl/fractionalOrders.pxd @@ -11,65 +11,16 @@ from PyNucleus_fem.functions cimport function from . twoPointFunctions cimport (twoPointFunction, constantTwoPoint, parametrizedTwoPointFunction) -from . interactionDomains cimport interactionDomain include "kernel_params_decl.pxi" -cdef class constantFractionalLaplacianScaling(constantTwoPoint): - cdef: - INDEX_t dim - REAL_t s, horizon - REAL_t tempered - - -cdef class constantFractionalLaplacianScalingBoundary(constantTwoPoint): - cdef: - INDEX_t dim - REAL_t s, horizon, tempered - - -cdef class constantFractionalLaplacianScalingDerivative(twoPointFunction): - cdef: - INDEX_t dim - REAL_t s - REAL_t horizon - REAL_t horizon2 - BOOL_t normalized - INDEX_t derivative - REAL_t tempered - REAL_t C - REAL_t fac - - -cdef class variableFractionalLaplacianScaling(parametrizedTwoPointFunction): - cdef: - INDEX_t dim - fractionalOrderBase sFun - function horizonFun - REAL_t facInfinite, facFinite - twoPointFunction phi - BOOL_t normalized - INDEX_t derivative - - -cdef class variableFractionalLaplacianScalingBoundary(parametrizedTwoPointFunction): - cdef: - INDEX_t dim - fractionalOrderBase sFun - function horizonFun - REAL_t facInfinite, facFinite - twoPointFunction phi - BOOL_t normalized - - -cdef class variableFractionalLaplacianScalingWithDifferentHorizon(variableFractionalLaplacianScaling): - pass - - cdef class fractionalOrderBase(twoPointFunction): cdef: public REAL_t min, max + public INDEX_t numParameters + cdef void evalGrad(self, REAL_t[::1] x, REAL_t[::1] y, REAL_t[::1] grad) + cdef REAL_t evalGradPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y, INDEX_t vectorSize, REAL_t* grad) cdef class constFractionalOrder(fractionalOrderBase): @@ -86,21 +37,19 @@ cdef class variableFractionalOrder(fractionalOrderBase): cdef class extendedFunction(function): cdef REAL_t eval(self, REAL_t[::1]) cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x) + cdef void evalGrad(self, REAL_t[::1] x, REAL_t[::1] grad) + cdef void evalGradPtr(self, INDEX_t dim, REAL_t* x, INDEX_t vectorSize, REAL_t* grad) cdef class singleVariableUnsymmetricFractionalOrder(variableFractionalOrder): cdef: public extendedFunction sFun + cdef void evalGrad(self, REAL_t[::1] x, REAL_t[::1] y, REAL_t[::1] grad) + cdef REAL_t evalGradPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y, INDEX_t vectorSize, REAL_t* grad) + + cdef class piecewiseConstantFractionalOrder(variableFractionalOrder): cdef: public function blockIndicator REAL_t[:, ::1] sVals - - -cdef class constantIntegrableScaling(constantTwoPoint): - cdef: - kernelType kType - INDEX_t dim - REAL_t horizon - interactionDomain interaction diff --git a/nl/PyNucleus_nl/fractionalOrders.pyx b/nl/PyNucleus_nl/fractionalOrders.pyx index 4b4222b..54bce36 100644 --- a/nl/PyNucleus_nl/fractionalOrders.pyx +++ b/nl/PyNucleus_nl/fractionalOrders.pyx @@ -13,25 +13,15 @@ from libc.math cimport (sin, cos, sinh, cosh, tanh, sqrt, atan, atan2, log, ceil, fabs as abs, M_PI as pi, pow, exp, erf) -from scipy.special.cython_special cimport psi as digamma -from scipy.special.cython_special cimport gamma as cgamma -from PyNucleus_base.myTypes import INDEX, REAL, ENCODE, BOOL +from PyNucleus_base.myTypes import INDEX, REAL, BOOL from PyNucleus_fem.functions cimport constant -from PyNucleus_fem.meshCy cimport meshBase -from PyNucleus_fem.DoFMaps cimport DoFMap +from PyNucleus_fem.meshCy cimport meshBase, cellFinder2 +from PyNucleus_fem.DoFMaps cimport shapeFunction, fe_vector, DoFMap from libc.stdlib cimport malloc from libc.string cimport memcpy -import warnings -from . interactionDomains cimport ball1, ball2, ballInf include "kernel_params.pxi" -cdef REAL_t inf = np.inf - - -cdef inline REAL_t gamma(REAL_t d): - return cgamma(d) - ###################################################################### @@ -58,19 +48,29 @@ cdef enum fracOrderParams: cdef class fractionalOrderBase(twoPointFunction): - def __init__(self, REAL_t smin, REAL_t smax, BOOL_t symmetric): + def __init__(self, REAL_t smin, REAL_t smax, BOOL_t symmetric, INDEX_t numParameters=0): super(fractionalOrderBase, self).__init__(symmetric) self.min = smin self.max = smax + self.numParameters = numParameters cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): raise NotImplementedError() + cdef void evalGrad(self, REAL_t[::1] x, REAL_t[::1] y, REAL_t[::1] grad): + raise NotImplementedError() + + cdef REAL_t evalGradPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y, INDEX_t vectorSize, REAL_t* grad): + raise NotImplementedError() + + def evalGrad_py(self, REAL_t[::1] x, REAL_t[::1] y, REAL_t[::1] grad): + self.evalGrad(x, y, grad) + def __getstate__(self): - return (self.min, self.max, self.symmetric) + return (self.min, self.max, self.symmetric, self.numParameters) def __setstate__(self, state): - fractionalOrderBase.__init__(self, state[0], state[1], state[2]) + fractionalOrderBase.__init__(self, state[0], state[1], state[2], state[3]) cdef class constFractionalOrder(fractionalOrderBase): @@ -95,8 +95,8 @@ cdef class constFractionalOrder(fractionalOrderBase): cdef class variableFractionalOrder(fractionalOrderBase): - def __init__(self, REAL_t smin, REAL_t smax, BOOL_t symmetric): - super(variableFractionalOrder, self).__init__(smin, smax, symmetric) + def __init__(self, REAL_t smin, REAL_t smax, BOOL_t symmetric, INDEX_t numParameters=0): + super(variableFractionalOrder, self).__init__(smin, smax, symmetric, numParameters) self.c_params = malloc(NUM_FRAC_ORDER_PARAMS*OFFSET) cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): @@ -123,10 +123,40 @@ cdef class extendedFunction(function): cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x): pass + cdef void evalGrad(self, REAL_t[::1] x, REAL_t[::1] grad): + raise NotImplementedError() + + cdef void evalGradPtr(self, INDEX_t dim, REAL_t* x, INDEX_t vectorSize, REAL_t* grad): + raise NotImplementedError() + + +cdef class singleVariableTwoPointFunction(twoPointFunction): + cdef: + extendedFunction fun + + def __init__(self, extendedFunction fun): + super(singleVariableTwoPointFunction, self).__init__(False) + self.fun = fun + + cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): + return self.fun.eval(x) + + cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y): + return self.fun.evalPtr(dim, x) + + def __repr__(self): + return '{}'.format(self.fun) + + def __getstate__(self): + return self.fun + + def __setstate__(self, state): + singleVariableTwoPointFunction.__init__(self, state) + cdef class singleVariableUnsymmetricFractionalOrder(variableFractionalOrder): - def __init__(self, extendedFunction sFun, REAL_t smin, REAL_t smax): - super(singleVariableUnsymmetricFractionalOrder, self).__init__(smin, smax, False) + def __init__(self, extendedFunction sFun, REAL_t smin, REAL_t smax, INDEX_t numParameters=0): + super(singleVariableUnsymmetricFractionalOrder, self).__init__(smin, smax, False, numParameters) self.sFun = sFun cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): @@ -135,14 +165,20 @@ cdef class singleVariableUnsymmetricFractionalOrder(variableFractionalOrder): cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y): return self.sFun.evalPtr(dim, x) + cdef void evalGrad(self, REAL_t[::1] x, REAL_t[::1] y, REAL_t[::1] grad): + self.sFun.evalGrad(x, grad) + + cdef REAL_t evalGradPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y, INDEX_t vectorSize, REAL_t* grad): + self.sFun.evalGradPtr(dim, x, vectorSize, grad) + def __repr__(self): return '{}({})'.format(self.__class__.__name__, self.sFun) def __getstate__(self): - return (self.sFun, self.min, self.max) + return (self.sFun, self.min, self.max, self.numParameters) def __setstate__(self, state): - singleVariableUnsymmetricFractionalOrder.__init__(self, state[0], state[1], state[2]) + singleVariableUnsymmetricFractionalOrder.__init__(self, state[0], state[1], state[2], state[3]) cdef REAL_t lambdaFractionalOrderFun(REAL_t *x, REAL_t *y, void *c_params): @@ -328,6 +364,12 @@ cdef class constantExtended(extendedFunction): cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x): return self.value + cdef void evalGrad(self, REAL_t[::1] x, REAL_t[::1] grad): + grad[0] = 1. + + cdef void evalGradPtr(self, INDEX_t dim, REAL_t* x, INDEX_t vectorSize, REAL_t* grad): + grad[0] = 1. + def __repr__(self): return '{}'.format(self.value) @@ -363,7 +405,11 @@ cdef class smoothLeftRight(extendedFunction): cdef class smoothStep(extendedFunction): cdef: - REAL_t sl, sr, r, slope, interface + public REAL_t sl + public REAL_t sr + public REAL_t r + public REAL_t slope + public REAL_t interface def __init__(self, REAL_t sl, REAL_t sr, REAL_t r, REAL_t interface=0.): self.sl = sl @@ -386,6 +432,30 @@ cdef class smoothStep(extendedFunction): return self.sr return self.sl + (self.sr-self.sl) * (3.0*pow((x[0]-self.interface)*self.slope+0.5, 2.0) - 2.0*pow((x[0]-self.interface)*self.slope+0.5, 3.0)) + cdef void evalGrad(self, REAL_t[::1] x, REAL_t[::1] grad): + cdef: + REAL_t fac + if x[0] < self.interface-self.r: + fac = 0. + elif x[0] > self.interface+self.r: + fac = 1. + else: + fac = (3.0*pow((x[0]-self.interface)*self.slope+0.5, 2.0) - 2.0*pow((x[0]-self.interface)*self.slope+0.5, 3.0)) + grad[0] = 1.-fac + grad[1] = fac + + cdef void evalGradPtr(self, INDEX_t dim, REAL_t* x, INDEX_t vectorSize, REAL_t* grad): + cdef: + REAL_t fac + if x[0] < self.interface-self.r: + fac = 0. + elif x[0] > self.interface+self.r: + fac = 1. + else: + fac = (3.0*pow((x[0]-self.interface)*self.slope+0.5, 2.0) - 2.0*pow((x[0]-self.interface)*self.slope+0.5, 3.0)) + grad[0] = 1.-fac + grad[1] = fac + def __repr__(self): return '{}(sl={},sr={},r={},interface={})'.format(self.__class__.__name__, self.sl, self.sr, self.r, self.interface) @@ -415,6 +485,30 @@ cdef class linearStep(extendedFunction): return self.sr return self.sl + self.slope*(x[0]-self.interface+self.r) + cdef void evalGrad(self, REAL_t[::1] x, REAL_t[::1] grad): + cdef: + REAL_t fac + if x[0] < self.interface-self.r: + fac = 0. + elif x[0] > self.interface+self.r: + fac = 1. + else: + fac = 0.5*(x[0]-self.interface+self.r)/self.r + grad[0] = 1-fac + grad[1] = fac + + cdef void evalGradPtr(self, INDEX_t dim, REAL_t* x, INDEX_t vectorSize, REAL_t* grad): + cdef: + REAL_t fac + if x[0] < self.interface-self.r: + fac = 0. + elif x[0] > self.interface+self.r: + fac = 1. + else: + fac = 0.5*(x[0]-self.interface+self.r)/self.r + grad[0] = 1-fac + grad[1] = fac + def __repr__(self): return '{}(sl={},sr={},r={})'.format(self.__class__.__name__, self.sl, self.sr, self.r) @@ -460,13 +554,99 @@ cdef class smoothStepRadial(extendedFunction): return '{}(sl={},sr={},r={},radius={})'.format(self.__class__.__name__, self.sl, self.sr, self.r, self.radius) + +cdef class lookupExtended(extendedFunction): + cdef: + meshBase mesh + public DoFMap dm + public REAL_t[::1] u + public cellFinder2 cellFinder + + def __init__(self, meshBase mesh, DoFMap dm, REAL_t[::1] u, cellFinder2 cF=None): + self.mesh = mesh + self.dm = dm + self.u = u + if cF is None: + self.cellFinder = cellFinder2(self.mesh) + else: + self.cellFinder = cF + + cdef REAL_t eval(self, REAL_t[::1] x): + cdef: + shapeFunction shapeFun + REAL_t val + INDEX_t cellNo, dof, k + cellNo = self.cellFinder.findCell(x) + if cellNo == -1: + return 0. + val = 0. + for k in range(self.dm.dofs_per_element): + dof = self.dm.cell2dof(cellNo, k) + if dof >= 0: + shapeFun = self.dm.localShapeFunctions[k] + val += shapeFun.eval(self.cellFinder.bary)*self.u[dof] + return val + + cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x): + cdef: + shapeFunction shapeFun + REAL_t val + INDEX_t cellNo, dof, k + cellNo = self.cellFinder.findCellPtr(x) + if cellNo == -1: + return 0. + val = 0. + for k in range(self.dm.dofs_per_element): + dof = self.dm.cell2dof(cellNo, k) + if dof >= 0: + shapeFun = self.dm.localShapeFunctions[k] + val += shapeFun.eval(self.cellFinder.bary)*self.u[dof] + return val + + cdef void evalGrad(self, REAL_t[::1] x, REAL_t[::1] grad): + cdef: + shapeFunction shapeFun + INDEX_t cellNo, dof, k + cellNo = self.cellFinder.findCell(x) + if cellNo == -1: + return + for dof in range(self.dm.num_dofs): + grad[dof] = 0. + for k in range(self.dm.dofs_per_element): + dof = self.dm.cell2dof(cellNo, k) + if dof >= 0: + shapeFun = self.dm.localShapeFunctions[k] + grad[dof] += shapeFun.eval(self.cellFinder.bary) + + cdef void evalGradPtr(self, INDEX_t dim, REAL_t* x, INDEX_t vectorSize, REAL_t* grad): + cdef: + shapeFunction shapeFun + INDEX_t cellNo, dof, k + cellNo = self.cellFinder.findCellPtr(x) + if cellNo == -1: + return + for dof in range(self.dm.num_dofs): + grad[dof] = 0. + for k in range(self.dm.dofs_per_element): + dof = self.dm.cell2dof(cellNo, k) + if dof >= 0: + shapeFun = self.dm.localShapeFunctions[k] + grad[dof] += shapeFun.eval(self.cellFinder.bary) + + def __getstate__(self): + return self.mesh, self.dm, self.u + + def __setstate__(self, state): + lookupExtended.__init__(self, state[0], state[1], state[2]) + + cdef class constantNonSymFractionalOrder(singleVariableUnsymmetricFractionalOrder): cdef: public REAL_t value def __init__(self, REAL_t s): sFun = constantExtended(s) - super(constantNonSymFractionalOrder, self).__init__(sFun, s, s) + super(constantNonSymFractionalOrder, self).__init__(sFun, s, s, 1) self.value = s @@ -474,13 +654,13 @@ cdef class smoothedLeftRightFractionalOrder(singleVariableUnsymmetricFractionalO def __init__(self, REAL_t sl, REAL_t sr, REAL_t r=0.1, REAL_t slope=200., REAL_t interface=0.): sFun = smoothStep(sl, sr, r, interface) # sFun = smoothLeftRight(sl, sr, r, slope) - super(smoothedLeftRightFractionalOrder, self).__init__(sFun, min(sl, sr), max(sl, sr)) + super(smoothedLeftRightFractionalOrder, self).__init__(sFun, min(sl, sr), max(sl, sr), 2) cdef class linearLeftRightFractionalOrder(singleVariableUnsymmetricFractionalOrder): def __init__(self, REAL_t sl, REAL_t sr, REAL_t r=0.1, REAL_t interface=0.): sFun = linearStep(sl, sr, r, interface) - super(linearLeftRightFractionalOrder, self).__init__(sFun, min(sl, sr), max(sl, sr)) + super(linearLeftRightFractionalOrder, self).__init__(sFun, min(sl, sr), max(sl, sr), 2) cdef class smoothedInnerOuterFractionalOrder(singleVariableUnsymmetricFractionalOrder): @@ -489,6 +669,23 @@ cdef class smoothedInnerOuterFractionalOrder(singleVariableUnsymmetricFractional super(smoothedInnerOuterFractionalOrder, self).__init__(sFun, min(sl, sr), max(sl, sr)) +cdef class feFractionalOrder(singleVariableUnsymmetricFractionalOrder): + cdef: + fe_vector vec + + def __init__(self, fe_vector vec, REAL_t smin, REAL_t smax): + self.vec = vec + sFun = lookupExtended(vec.dm.mesh, vec.dm, vec) + super(feFractionalOrder, self).__init__(sFun, smin, smax, numParameters=vec.dm.num_dofs) + + def __getstate__(self): + return (self.vec, self.vec.dm, self.min, self.max) + + def __setstate__(self, state): + vec = fe_vector(state[0], state[1]) + feFractionalOrder.__init__(self, vec, state[2], state[3]) + + cdef REAL_t innerOuterFractionalOrderFun(REAL_t *x, REAL_t *y, void *c_params): cdef: INDEX_t dim = getINDEX(c_params, fDIM) @@ -717,330 +914,3 @@ cdef class layersFractionalOrder(variableFractionalOrder): def __repr__(self): numLayers = getINDEX(self.c_params, fR) return '{}(numLayers={})'.format(self.__class__.__name__, numLayers) - - -###################################################################### - -cdef class constantFractionalLaplacianScaling(constantTwoPoint): - def __init__(self, INDEX_t dim, REAL_t s, REAL_t horizon, REAL_t tempered): - self.dim = dim - if 1. < s and s < 2.: - s = s-1. - self.s = s - self.horizon = horizon - self.tempered = tempered - if (self.horizon <= 0.) or (self.s <= 0.) or (self.s >= 1.): - value = np.nan - else: - if horizon < inf: - value = (2.-2*s) * pow(horizon, 2*s-2.) * gamma(0.5*dim)/pow(pi, 0.5*dim) * 0.5 - if dim > 1: - value *= 2. - else: - if (tempered == 0.) or (s == 0.5): - value = 2.0**(2.0*s) * s * gamma(s+0.5*dim)/pow(pi, 0.5*dim)/gamma(1.0-s) * 0.5 - else: - value = gamma(0.5*dim) / abs(gamma(-2*s))/pow(pi, 0.5*dim) * 0.5 * 0.5 - super(constantFractionalLaplacianScaling, self).__init__(value) - - def __getstate__(self): - return (self.dim, self.s, self.horizon, self.tempered) - - def __setstate__(self, state): - constantFractionalLaplacianScaling.__init__(self, state[0], state[1], state[2], state[3]) - - def __repr__(self): - return '{}({},{} -> {})'.format(self.__class__.__name__, self.s, self.horizon, self.value) - - -cdef class constantFractionalLaplacianScalingDerivative(twoPointFunction): - def __init__(self, INDEX_t dim, REAL_t s, REAL_t horizon, BOOL_t normalized, INDEX_t derivative, REAL_t tempered): - super(constantFractionalLaplacianScalingDerivative, self).__init__(True) - - self.dim = dim - self.s = s - self.horizon = horizon - self.normalized = normalized - self.derivative = derivative - self.tempered = tempered - - horizon2 = horizon**2 - self.horizon2 = horizon2 - - if self.normalized: - if horizon2 < inf: - self.C = (2.-2*s) * pow(horizon2, s-1.) * gamma(0.5*dim)/pow(pi, 0.5*dim) * 0.5 - if dim > 1: - self.C *= 2. - else: - if (tempered == 0.) or (s == 0.5): - self.C = 2.0**(2.0*s) * s * gamma(s+0.5*dim) * pow(pi, -0.5*dim) / gamma(1.0-s) * 0.5 - else: - self.C = gamma(0.5*dim) / abs(gamma(-2*s))/pow(pi, 0.5*dim) * 0.5 * 0.5 - else: - self.C = 0.5 - - if self.derivative == 1: - if horizon2 < inf: - self.fac = 1./(s-1.) - else: - self.fac = digamma(s+0.5*self.dim) + digamma(-s) - else: - raise NotImplementedError(self.derivative) - - cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): - cdef: - REAL_t d2 - INDEX_t i - - if self.derivative == 0: - return self.C - elif self.derivative == 1: - d2 = 0. - for i in range(self.dim): - d2 += (x[i]-y[i])*(x[i]-y[i]) - if self.normalized: - if self.horizon2 < inf: - return self.C*(-log(d2/self.horizon2) + self.fac) - else: - return self.C*(-log(0.25*d2) + self.fac) - else: - return -self.C*log(d2) - else: - raise NotImplementedError() - - cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y): - cdef: - REAL_t d2 - INDEX_t i - - if self.derivative == 0: - return self.C - elif self.derivative == 1: - d2 = 0. - for i in range(self.dim): - d2 += (x[i]-y[i])*(x[i]-y[i]) - if self.normalized: - if self.horizon2 < inf: - return self.C*(-log(d2/self.horizon2) + self.fac) - else: - return self.C*(-log(0.25*d2) + self.fac) - else: - return -self.C*log(d2) - else: - raise NotImplementedError() - - def __getstate__(self): - return (self.dim, self.s, self.horizon, self.normalized, self.derivative, self.tempered) - - def __setstate__(self, state): - constantFractionalLaplacianScalingDerivative.__init__(self, state[0], state[1], state[2], state[3], state[4], state[5]) - - -cdef class constantIntegrableScaling(constantTwoPoint): - def __init__(self, kernelType kType, interactionDomain interaction, INDEX_t dim, REAL_t horizon): - self.kType = kType - self.dim = dim - self.interaction = interaction - self.horizon = horizon - if self.horizon <= 0.: - value = np.nan - else: - if kType == INDICATOR: - if dim == 1: - value = 3./horizon**3 / 2. - elif dim == 2: - if isinstance(self.interaction, ball2): - value = 8./pi/horizon**4 / 2. - elif isinstance(self.interaction, ballInf): - value = 3./4./horizon**4 / 2. - else: - raise NotImplementedError() - else: - raise NotImplementedError() - elif kType == PERIDYNAMIC: - if dim == 1: - value = 2./horizon**2 / 2. - elif dim == 2: - if isinstance(self.interaction, ball2): - value = 6./pi/horizon**3 / 2. - else: - raise NotImplementedError() - else: - raise NotImplementedError() - elif kType == GAUSSIAN: - if dim == 1: - # value = 4.0/sqrt(pi)/(horizon/3.)**3 / 2. - value = 4.0/sqrt(pi)/(erf(3.0)-6.0*exp(-9.0)/sqrt(pi))/(horizon/3.0)**3 / 2. - elif dim == 2: - if isinstance(self.interaction, ball2): - # value = 4.0/pi/(horizon/3.0)**4 / 2. - value = 4.0/pi/(1.0-10.0*exp(-9.0))/(horizon/3.0)**4 / 2. - else: - raise NotImplementedError() - else: - raise NotImplementedError() - else: - raise NotImplementedError() - super(constantIntegrableScaling, self).__init__(value) - - def __getstate__(self): - return (self.kType, self.interaction, self.dim, self.horizon) - - def __setstate__(self, state): - constantIntegrableScaling.__init__(self, state[0], state[1], state[2], state[3]) - - def __repr__(self): - return '{}({} -> {})'.format(self.__class__.__name__, self.horizon, self.value) - - -cdef class variableFractionalLaplacianScaling(parametrizedTwoPointFunction): - def __init__(self, BOOL_t symmetric, BOOL_t normalized=True, INDEX_t derivative=0): - super(variableFractionalLaplacianScaling, self).__init__(symmetric) - self.normalized = normalized - self.derivative = derivative - - cdef void setParams(self, void *params): - parametrizedTwoPointFunction.setParams(self, params) - self.dim = getINDEX(self.params, fKDIM) - - cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): - cdef: - REAL_t s = getREAL(self.params, fS) - REAL_t horizon2 = getREAL(self.params, fHORIZON2) - REAL_t C, d2 - INDEX_t i - - if self.normalized: - if horizon2 < inf: - if self.dim == 1: - C = (2.-2*s) * pow(horizon2, s-1.) * 0.5 - elif self.dim == 2: - C = (2.-2*s) * pow(horizon2, s-1.) * 2./pi * 0.5 - elif self.dim == 3: - C = (2.-2*s) * pow(horizon2, s-1.) * 1./pi * 0.5 - else: - raise NotImplementedError() - else: - C = 2.0**(2.0*s) * s * gamma(s+0.5*self.dim) * pow(pi, -0.5*self.dim) / gamma(1.0-s) * 0.5 - else: - C = 0.5 - - if self.derivative == 0: - return C - elif self.derivative == 1: - d2 = 0. - for i in range(self.dim): - d2 += (x[i]-y[i])*(x[i]-y[i]) - if self.normalized: - if horizon2 < inf: - return C*(-log(d2/horizon2) + 1./(s-1.)) - else: - return C*(-log(0.25*d2) + digamma(s+0.5*self.dim) + digamma(-s)) - else: - return -C*log(d2) - else: - raise NotImplementedError() - - cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y): - cdef: - REAL_t s = getREAL(self.params, fS) - REAL_t horizon2 = getREAL(self.params, fHORIZON2) - REAL_t C, d2 - INDEX_t i - - if self.normalized: - if horizon2 < inf: - if self.dim == 1: - C = (2.-2*s) * pow(horizon2, s-1.) * 0.5 - elif self.dim == 2: - C = (2.-2*s) * pow(horizon2, s-1.) * 2./pi * 0.5 - elif self.dim == 3: - C = (2.-2*s) * pow(horizon2, s-1.) * 1./pi * 0.5 - else: - raise NotImplementedError() - else: - C = 2.0**(2.0*s) * s * gamma(s+0.5*self.dim) * pow(pi, -0.5*self.dim) / gamma(1.0-s) * 0.5 - else: - C = 0.5 - - if self.derivative == 0: - return C - elif self.derivative == 1: - d2 = 0. - for i in range(self.dim): - d2 += (x[i]-y[i])*(x[i]-y[i]) - if self.normalized: - if horizon2 < inf: - return C*(-log(d2/horizon2) + 1./(s-1.)) - else: - return C*(-log(0.25*d2) + digamma(s+0.5*self.dim) + digamma(-s)) - else: - return -C*log(d2) - else: - raise NotImplementedError() - - def getScalingWithDifferentHorizon(self): - cdef: - variableFractionalLaplacianScalingWithDifferentHorizon scaling - function horizonFun - BOOL_t horizonFunNull = isNull(self.params, fHORIZONFUN) - if not horizonFunNull: - horizonFun = (((self.params+fHORIZONFUN))[0]) - else: - horizonFun = constant(sqrt(getREAL(self.params, fHORIZON2))) - scaling = variableFractionalLaplacianScalingWithDifferentHorizon(self.symmetric, self.normalized, self.derivative, horizonFun) - return scaling - - def __repr__(self): - return 'variableFractionalLaplacianScaling' - - def __getstate__(self): - return (self.symmetric, self.normalized, self.derivative) - - def __setstate__(self, state): - variableFractionalLaplacianScaling.__init__(self, state[0], state[1], state[2]) - - -###################################################################### - - -cdef class variableFractionalLaplacianScalingWithDifferentHorizon(variableFractionalLaplacianScaling): - def __init__(self, BOOL_t symmetric, BOOL_t normalized, INDEX_t derivative, function horizonFun): - super(variableFractionalLaplacianScalingWithDifferentHorizon, self).__init__(symmetric, normalized, derivative) - self.horizonFun = horizonFun - - cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): - cdef: - void* params - void* paramsModified = malloc(NUM_KERNEL_PARAMS*OFFSET) - REAL_t horizon - horizon = self.horizonFun.eval(x) - params = self.getParams() - memcpy(paramsModified, params, NUM_KERNEL_PARAMS*OFFSET) - setREAL(paramsModified, fHORIZON2, horizon**2) - self.setParams(paramsModified) - scalingValue = variableFractionalLaplacianScaling.eval(self, x, y) - self.setParams(params) - return scalingValue - - cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y): - cdef: - void* params - void* paramsModified = malloc(NUM_KERNEL_PARAMS*OFFSET) - REAL_t horizon - REAL_t[::1] xA = x - horizon = self.horizonFun.eval(xA) - params = self.getParams() - memcpy(paramsModified, params, NUM_KERNEL_PARAMS*OFFSET) - setREAL(paramsModified, fHORIZON2, horizon**2) - self.setParams(paramsModified) - scalingValue = variableFractionalLaplacianScaling.evalPtr(self, dim, x, y) - self.setParams(params) - return scalingValue - - def __getstate__(self): - return (self.symmetric, self.normalized, self.derivative, self.horizonFun) - - def __setstate__(self, state): - variableFractionalLaplacianScalingWithDifferentHorizon.__init__(self, state[0], state[1], state[2], state[3]) diff --git a/nl/PyNucleus_nl/helpers.py b/nl/PyNucleus_nl/helpers.py index e4c0264..7338409 100644 --- a/nl/PyNucleus_nl/helpers.py +++ b/nl/PyNucleus_nl/helpers.py @@ -115,6 +115,8 @@ def processBC(tag, boundaryCondition, kernel): def getFracLapl(mesh, DoFMap, kernel=None, rangedOpParams={}, **kwargs): + if kernel is None and len(rangedOpParams) == 0: + return DoFMap.assembleStiffness(dm2=kwargs.get('dm2', None)) assert kernel is not None or 's' in rangedOpParams, (kernel, rangedOpParams) boundaryCondition = kwargs.get('boundaryCondition', 'Dirichlet') @@ -256,7 +258,14 @@ def getFracLapl(mesh, DoFMap, kernel=None, rangedOpParams={}, **kwargs): params['genKernel'] = kwargs['genKernel'] if kernel is None: kernel = getFractionalKernel(mesh.dim, s, constant(horizon.ranges[0, 0]), scaling=scaling, normalized=normalized) - builder = nonlocalBuilder(mesh, DoFMap, kernel, params, zeroExterior=zeroExterior, comm=comm, logging=logging, PLogger=PLogger) + dm2 = kwargs.pop('dm2', None) + if dm2 is not None and matrixFormat.upper() in ('H2', 'SPARSE', 'SPARSIFIED') and DoFMap.num_boundary_dofs > 0: + # currently not implemented + dm, R_interior, R_bc = DoFMap.getFullDoFMap(dm2) + A = getFracLapl(mesh, dm, kernel, rangedOpParams={}, **kwargs) + A = R_interior*A*R_bc.transpose() + return A + builder = nonlocalBuilder(mesh, DoFMap, kernel, params, zeroExterior=zeroExterior, comm=comm, logging=logging, PLogger=PLogger, dm2=dm2) if diagonal: with timer('Assemble diagonal matrix {}, zeroExterior={}'.format(kernel, zeroExterior)): A = builder.getDiagonal() diff --git a/nl/PyNucleus_nl/kernelNormalization.pxd b/nl/PyNucleus_nl/kernelNormalization.pxd new file mode 100644 index 0000000..3a3b124 --- /dev/null +++ b/nl/PyNucleus_nl/kernelNormalization.pxd @@ -0,0 +1,86 @@ +################################################################################### +# Copyright 2021 National Technology & Engineering Solutions of Sandia, # +# LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the # +# U.S. Government retains certain rights in this software. # +# If you want to use this code, please refer to the README.rst and LICENSE files. # +################################################################################### + +cimport numpy as np +from PyNucleus_base.myTypes cimport INDEX_t, REAL_t, BOOL_t +from PyNucleus_fem.functions cimport function +from . twoPointFunctions cimport (twoPointFunction, + constantTwoPoint, + parametrizedTwoPointFunction) +from . fractionalOrders cimport fractionalOrderBase +from . interactionDomains cimport interactionDomain + +include "kernel_params_decl.pxi" + + +cdef class memoizedFun: + cdef: + dict memory + int hit, miss + cdef REAL_t eval(self, REAL_t x) + + +cdef class constantFractionalLaplacianScaling(constantTwoPoint): + cdef: + INDEX_t dim + REAL_t s, horizon + REAL_t tempered + + +cdef class constantFractionalLaplacianScalingBoundary(constantTwoPoint): + cdef: + INDEX_t dim + REAL_t s, horizon, tempered + + +cdef class constantFractionalLaplacianScalingDerivative(twoPointFunction): + cdef: + INDEX_t dim + REAL_t s + REAL_t horizon + REAL_t horizon2 + BOOL_t normalized + BOOL_t boundary + INDEX_t derivative + REAL_t tempered + REAL_t C + REAL_t fac + + +cdef class variableFractionalLaplacianScaling(parametrizedTwoPointFunction): + cdef: + INDEX_t dim + fractionalOrderBase sFun + function horizonFun + REAL_t facInfinite, facFinite + twoPointFunction phi + BOOL_t normalized + BOOL_t boundary + INDEX_t derivative + memoizedFun digamma + + +cdef class variableFractionalLaplacianScalingBoundary(parametrizedTwoPointFunction): + cdef: + INDEX_t dim + fractionalOrderBase sFun + function horizonFun + REAL_t facInfinite, facFinite + twoPointFunction phi + BOOL_t normalized + + +cdef class variableFractionalLaplacianScalingWithDifferentHorizon(variableFractionalLaplacianScaling): + pass + + +cdef class constantIntegrableScaling(constantTwoPoint): + cdef: + kernelType kType + INDEX_t dim + REAL_t horizon + interactionDomain interaction diff --git a/nl/PyNucleus_nl/kernelNormalization.pyx b/nl/PyNucleus_nl/kernelNormalization.pyx new file mode 100644 index 0000000..d54c678 --- /dev/null +++ b/nl/PyNucleus_nl/kernelNormalization.pyx @@ -0,0 +1,418 @@ +################################################################################### +# Copyright 2021 National Technology & Engineering Solutions of Sandia, # +# LLC (NTESS). Under the terms of Contract DE-NA0003525 with NTESS, the # +# U.S. Government retains certain rights in this software. # +# If you want to use this code, please refer to the README.rst and LICENSE files. # +################################################################################### + +"""Defines normalizations for different types of kernels.""" +import numpy as np +cimport numpy as np +from libc.math cimport (sin, cos, sinh, cosh, tanh, sqrt, atan, atan2, + log, ceil, + fabs as abs, M_PI as pi, pow, + exp, erf) +from scipy.special.cython_special cimport psi as digamma +from scipy.special.cython_special cimport gamma as cgamma + +from PyNucleus_fem.functions cimport constant +from libc.stdlib cimport malloc +from libc.string cimport memcpy +from . interactionDomains cimport ball1, ball2, ballInf + +include "kernel_params.pxi" + +cdef REAL_t inf = np.inf + + +cdef inline REAL_t gamma(REAL_t d): + return cgamma(d) + + +cdef class memoizedFun: + def __init__(self): + self.memory = dict() + self.hit = 0 + self.miss = 0 + + cdef REAL_t eval(self, REAL_t x): + raise NotImplementedError() + + def __call__(self, REAL_t x): + return self.eval(x) + + def stats(self): + print(len(self.memory), self.hit, self.miss) + + +cdef class memoizedDigamma(memoizedFun): + cdef REAL_t eval(self, REAL_t x): + cdef REAL_t val + try: + val = self.memory[x] + self.hit += 1 + return val + except KeyError: + self.miss += 1 + val = digamma(x) + self.memory[x] = val + return val + + +cdef class constantFractionalLaplacianScaling(constantTwoPoint): + def __init__(self, INDEX_t dim, REAL_t s, REAL_t horizon, REAL_t tempered): + self.dim = dim + if 1. < s and s < 2.: + s = s-1. + self.s = s + self.horizon = horizon + self.tempered = tempered + if (self.horizon <= 0.) or (self.s <= 0.) or (self.s >= 1.): + value = np.nan + else: + if horizon < inf: + value = (2.-2*s) * pow(horizon, 2*s-2.) * gamma(0.5*dim)/pow(pi, 0.5*dim) * 0.5 + if dim > 1: + value *= 2. + else: + if (tempered == 0.) or (s == 0.5): + value = 2.0**(2.0*s) * s * gamma(s+0.5*dim)/pow(pi, 0.5*dim)/gamma(1.0-s) * 0.5 + else: + value = gamma(0.5*dim) / abs(gamma(-2*s))/pow(pi, 0.5*dim) * 0.5 * 0.5 + super(constantFractionalLaplacianScaling, self).__init__(value) + + def __getstate__(self): + return (self.dim, self.s, self.horizon, self.tempered) + + def __setstate__(self, state): + constantFractionalLaplacianScaling.__init__(self, state[0], state[1], state[2], state[3]) + + def __repr__(self): + return '{}({},{} -> {})'.format(self.__class__.__name__, self.s, self.horizon, self.value) + + +cdef class constantFractionalLaplacianScalingDerivative(twoPointFunction): + def __init__(self, INDEX_t dim, REAL_t s, REAL_t horizon, BOOL_t normalized, BOOL_t boundary, INDEX_t derivative, REAL_t tempered): + super(constantFractionalLaplacianScalingDerivative, self).__init__(True) + + self.dim = dim + self.s = s + self.horizon = horizon + self.normalized = normalized + self.boundary = boundary + self.derivative = derivative + self.tempered = tempered + + horizon2 = horizon**2 + self.horizon2 = horizon2 + + if self.normalized: + if horizon2 < inf: + self.C = (2.-2*s) * pow(horizon2, s-1.) * gamma(0.5*dim)/pow(pi, 0.5*dim) * 0.5 + if dim > 1: + self.C *= 2. + else: + if (tempered == 0.) or (s == 0.5): + self.C = 2.0**(2.0*s) * s * gamma(s+0.5*dim) * pow(pi, -0.5*dim) / gamma(1.0-s) * 0.5 + else: + self.C = gamma(0.5*dim) / abs(gamma(-2*s))/pow(pi, 0.5*dim) * 0.5 * 0.5 + else: + self.C = 0.5 + + if self.derivative == 1: + if self.normalized: + if horizon2 < inf: + if not self.boundary: + self.fac = -1./(1.-s) + else: + self.fac = -1./(1.-s) - 1./s + else: + if not self.boundary: + self.fac = digamma(s+0.5*self.dim) + digamma(-s) + else: + self.fac = digamma(s+0.5*self.dim) + digamma(1.-s) + else: + if not self.boundary: + self.fac = 0. + else: + self.fac = -1./s + else: + raise NotImplementedError(self.derivative) + + cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): + cdef: + REAL_t d2 + INDEX_t i + + if self.derivative == 0: + return self.C + elif self.derivative == 1: + d2 = 0. + for i in range(self.dim): + d2 += (x[i]-y[i])*(x[i]-y[i]) + if self.normalized: + if self.horizon2 < inf: + return self.C*(-log(d2/self.horizon2) + self.fac) + else: + return self.C*(-log(0.25*d2) + self.fac) + else: + return self.C*(-log(d2) + self.fac) + else: + raise NotImplementedError() + + cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y): + cdef: + REAL_t d2 + INDEX_t i + + if self.derivative == 0: + return self.C + elif self.derivative == 1: + d2 = 0. + for i in range(self.dim): + d2 += (x[i]-y[i])*(x[i]-y[i]) + if self.normalized: + if self.horizon2 < inf: + return self.C*(-log(d2/self.horizon2) + self.fac) + else: + return self.C*(-log(0.25*d2) + self.fac) + else: + return self.C*(-log(d2) + self.fac) + else: + raise NotImplementedError() + + def __getstate__(self): + return (self.dim, self.s, self.horizon, self.normalized, self.boundary, self.derivative, self.tempered) + + def __setstate__(self, state): + constantFractionalLaplacianScalingDerivative.__init__(self, state[0], state[1], state[2], state[3], state[4], state[5], state[6]) + + +cdef class constantIntegrableScaling(constantTwoPoint): + def __init__(self, kernelType kType, interactionDomain interaction, INDEX_t dim, REAL_t horizon): + self.kType = kType + self.dim = dim + self.interaction = interaction + self.horizon = horizon + if self.horizon <= 0.: + value = np.nan + else: + if kType == INDICATOR: + if dim == 1: + value = 3./horizon**3 / 2. + elif dim == 2: + if isinstance(self.interaction, ball2): + value = 8./pi/horizon**4 / 2. + elif isinstance(self.interaction, ballInf): + value = 3./4./horizon**4 / 2. + else: + raise NotImplementedError() + else: + raise NotImplementedError() + elif kType == PERIDYNAMIC: + if dim == 1: + value = 2./horizon**2 / 2. + elif dim == 2: + if isinstance(self.interaction, ball2): + value = 6./pi/horizon**3 / 2. + else: + raise NotImplementedError() + else: + raise NotImplementedError() + elif kType == GAUSSIAN: + if dim == 1: + # value = 4.0/sqrt(pi)/(horizon/3.)**3 / 2. + value = 4.0/sqrt(pi)/(erf(3.0)-6.0*exp(-9.0)/sqrt(pi))/(horizon/3.0)**3 / 2. + elif dim == 2: + if isinstance(self.interaction, ball2): + # value = 4.0/pi/(horizon/3.0)**4 / 2. + value = 4.0/pi/(1.0-10.0*exp(-9.0))/(horizon/3.0)**4 / 2. + else: + raise NotImplementedError() + else: + raise NotImplementedError() + else: + raise NotImplementedError() + super(constantIntegrableScaling, self).__init__(value) + + def __getstate__(self): + return (self.kType, self.interaction, self.dim, self.horizon) + + def __setstate__(self, state): + constantIntegrableScaling.__init__(self, state[0], state[1], state[2], state[3]) + + def __repr__(self): + return '{}({} -> {})'.format(self.__class__.__name__, self.horizon, self.value) + + +cdef class variableFractionalLaplacianScaling(parametrizedTwoPointFunction): + def __init__(self, BOOL_t symmetric, BOOL_t normalized=True, BOOL_t boundary=False, INDEX_t derivative=0): + super(variableFractionalLaplacianScaling, self).__init__(symmetric) + self.normalized = normalized + self.boundary = boundary + self.derivative = derivative + self.digamma = memoizedDigamma() + + cdef void setParams(self, void *params): + parametrizedTwoPointFunction.setParams(self, params) + self.dim = getINDEX(self.params, fKDIM) + + cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): + cdef: + REAL_t s = getREAL(self.params, fS) + REAL_t horizon2 = getREAL(self.params, fHORIZON2) + REAL_t C, d2 + INDEX_t i + + if self.normalized: + if horizon2 < inf: + if self.dim == 1: + C = (2.-2*s) * pow(horizon2, s-1.) * 0.5 + elif self.dim == 2: + C = (2.-2*s) * pow(horizon2, s-1.) * 2./pi * 0.5 + elif self.dim == 3: + C = (2.-2*s) * pow(horizon2, s-1.) * 1./pi * 0.5 + else: + raise NotImplementedError() + else: + C = 2.0**(2.0*s) * s * gamma(s+0.5*self.dim) * pow(pi, -0.5*self.dim) / gamma(1.0-s) * 0.5 + else: + C = 0.5 + + if self.derivative == 0: + return C + elif self.derivative == 1: + d2 = 0. + for i in range(self.dim): + d2 += (x[i]-y[i])*(x[i]-y[i]) + if self.normalized: + if horizon2 < inf: + if not self.boundary: + return C*(-log(d2/horizon2) - 1./(1.-s)) + else: + return C*(-log(d2/horizon2) - 1./(1.-s) - 1./s) + else: + if not self.boundary: + return C*(-log(0.25*d2) + self.digamma.eval(s+0.5*self.dim) + self.digamma.eval(-s)) + else: + return C*(-log(0.25*d2) + self.digamma.eval(s+0.5*self.dim) + self.digamma.eval(1.-s)) + else: + if not self.boundary: + return C*(-log(d2)) + else: + return C*(-log(d2)-1./s) + else: + raise NotImplementedError() + + cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y): + cdef: + REAL_t s = getREAL(self.params, fS) + REAL_t horizon2 = getREAL(self.params, fHORIZON2) + REAL_t C, d2 + INDEX_t i + + if self.normalized: + if horizon2 < inf: + if self.dim == 1: + C = (2.-2*s) * pow(horizon2, s-1.) * 0.5 + elif self.dim == 2: + C = (2.-2*s) * pow(horizon2, s-1.) * 2./pi * 0.5 + elif self.dim == 3: + C = (2.-2*s) * pow(horizon2, s-1.) * 1./pi * 0.5 + else: + raise NotImplementedError() + else: + C = 2.0**(2.0*s) * s * gamma(s+0.5*self.dim) * pow(pi, -0.5*self.dim) / gamma(1.0-s) * 0.5 + else: + C = 0.5 + + if self.derivative == 0: + return C + elif self.derivative == 1: + d2 = 0. + for i in range(self.dim): + d2 += (x[i]-y[i])*(x[i]-y[i]) + if self.normalized: + if horizon2 < inf: + if not self.boundary: + return C*(-log(d2/horizon2) - 1./(1.-s)) + else: + return C*(-log(d2/horizon2) - 1./(1.-s) - 1./s) + else: + if not self.boundary: + return C*(-log(0.25*d2) + self.digamma.eval(s+0.5*self.dim) + self.digamma.eval(-s)) + else: + return C*(-log(0.25*d2) + self.digamma.eval(s+0.5*self.dim) + self.digamma.eval(1.-s)) + else: + if not self.boundary: + return C*(-log(d2)) + else: + return C*(-log(d2)-1./s) + + else: + raise NotImplementedError() + + def getScalingWithDifferentHorizon(self): + cdef: + variableFractionalLaplacianScalingWithDifferentHorizon scaling + function horizonFun + BOOL_t horizonFunNull = isNull(self.params, fHORIZONFUN) + if not horizonFunNull: + horizonFun = (((self.params+fHORIZONFUN))[0]) + else: + horizonFun = constant(sqrt(getREAL(self.params, fHORIZON2))) + scaling = variableFractionalLaplacianScalingWithDifferentHorizon(self.symmetric, self.normalized, self.boundary, self.derivative, horizonFun) + return scaling + + def __repr__(self): + return 'variableFractionalLaplacianScaling' + + def __getstate__(self): + return (self.symmetric, self.normalized, self.boundary, self.derivative) + + def __setstate__(self, state): + variableFractionalLaplacianScaling.__init__(self, state[0], state[1], state[2], state[3]) + + +###################################################################### + + +cdef class variableFractionalLaplacianScalingWithDifferentHorizon(variableFractionalLaplacianScaling): + def __init__(self, BOOL_t symmetric, BOOL_t normalized, BOOL_t boundary, INDEX_t derivative, function horizonFun): + super(variableFractionalLaplacianScalingWithDifferentHorizon, self).__init__(symmetric, normalized, boundary, derivative) + self.horizonFun = horizonFun + + cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): + cdef: + void* params + void* paramsModified = malloc(NUM_KERNEL_PARAMS*OFFSET) + REAL_t horizon + horizon = self.horizonFun.eval(x) + params = self.getParams() + memcpy(paramsModified, params, NUM_KERNEL_PARAMS*OFFSET) + setREAL(paramsModified, fHORIZON2, horizon**2) + self.setParams(paramsModified) + scalingValue = variableFractionalLaplacianScaling.eval(self, x, y) + self.setParams(params) + return scalingValue + + cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y): + cdef: + void* params + void* paramsModified = malloc(NUM_KERNEL_PARAMS*OFFSET) + REAL_t horizon + REAL_t[::1] xA = x + horizon = self.horizonFun.eval(xA) + params = self.getParams() + memcpy(paramsModified, params, NUM_KERNEL_PARAMS*OFFSET) + setREAL(paramsModified, fHORIZON2, horizon**2) + self.setParams(paramsModified) + scalingValue = variableFractionalLaplacianScaling.evalPtr(self, dim, x, y) + self.setParams(params) + return scalingValue + + def __getstate__(self): + return (self.symmetric, self.normalized, self.boundary, self.derivative, self.horizonFun) + + def __setstate__(self, state): + variableFractionalLaplacianScalingWithDifferentHorizon.__init__(self, state[0], state[1], state[2], state[3], state[4]) diff --git a/nl/PyNucleus_nl/kernels.py b/nl/PyNucleus_nl/kernels.py index 5a58fa7..2bf053c 100644 --- a/nl/PyNucleus_nl/kernels.py +++ b/nl/PyNucleus_nl/kernels.py @@ -18,10 +18,11 @@ from . fractionalOrders import (fractionalOrderBase, constFractionalOrder, variableConstFractionalOrder, - constantFractionalLaplacianScaling, - constantFractionalLaplacianScalingDerivative, - variableFractionalLaplacianScaling, - constantIntegrableScaling) + singleVariableTwoPointFunction) +from . kernelNormalization import (constantFractionalLaplacianScaling, + constantFractionalLaplacianScalingDerivative, + variableFractionalLaplacianScaling, + constantIntegrableScaling) from . kernelsCy import (Kernel, FractionalKernel, RangedFractionalKernel, @@ -129,10 +130,13 @@ def getFractionalKernel(dim, if piecewise: warnings.warn('Derivative kernels cannot be piecewise. Switching to piecewise == False.') piecewise = False - scaling = constantFractionalLaplacianScalingDerivative(dim, sFun.value, horizonFun.value, normalized, derivative, tempered) + scaling = constantFractionalLaplacianScalingDerivative(dim, sFun.value, horizonFun.value, normalized, boundary, derivative, tempered) else: symmetric = sFun.symmetric and isinstance(horizonFun, constant) - scaling = variableFractionalLaplacianScaling(symmetric, normalized, derivative) + if piecewise and isinstance(sFun, singleVariableTwoPointFunction): + warnings.warn('Variable s kernels cannot be piecewise. Switching to piecewise == False.') + piecewise = False + scaling = variableFractionalLaplacianScaling(symmetric, normalized, boundary, derivative) if boundary: if isinstance(sFun, (constFractionalOrder, variableConstFractionalOrder)): @@ -154,7 +158,8 @@ def getIntegrableKernel(dim, interaction=None, normalized=True, piecewise=True, - phi=None): + phi=None, + boundary=False): dim_ = _getDim(dim) kType = _getKernelType(kernel) horizonFun = _getHorizon(horizon) @@ -168,7 +173,7 @@ def getIntegrableKernel(dim, raise NotImplementedError() else: scaling = constantTwoPoint(0.5) - return Kernel(dim_, kType=kType, horizon=horizonFun, interaction=interaction, scaling=scaling, phi=phi, piecewise=piecewise) + return Kernel(dim_, kType=kType, horizon=horizonFun, interaction=interaction, scaling=scaling, phi=phi, piecewise=piecewise, boundary=boundary) def getKernel(dim, diff --git a/nl/PyNucleus_nl/kernelsCy.pxd b/nl/PyNucleus_nl/kernelsCy.pxd index 2949727..e3ac4a1 100644 --- a/nl/PyNucleus_nl/kernelsCy.pxd +++ b/nl/PyNucleus_nl/kernelsCy.pxd @@ -34,6 +34,8 @@ cdef class Kernel(twoPointFunction): public BOOL_t variableScaling public BOOL_t variable public BOOL_t piecewise + public BOOL_t boundary + public INDEX_t vectorSize kernel_fun_t kernelFun void *c_kernel_params cdef REAL_t getSingularityValue(self) @@ -46,13 +48,15 @@ cdef class Kernel(twoPointFunction): cdef void evalParams(self, REAL_t[::1] x, REAL_t[::1] y) cdef void evalParamsPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y) cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y) + cdef void evalVector(self, REAL_t[::1] x, REAL_t[::1] y, REAL_t[::1] vec) + cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y) + cdef void evalVectorPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y, INDEX_t vectorSize, REAL_t* vec) cdef class FractionalKernel(Kernel): cdef: public fractionalOrderBase s public BOOL_t variableOrder - public BOOL_t boundary public INDEX_t derivative cdef REAL_t getsValue(self) cdef void setsValue(self, REAL_t s) diff --git a/nl/PyNucleus_nl/kernelsCy.pyx b/nl/PyNucleus_nl/kernelsCy.pyx index 693da83..684ef30 100644 --- a/nl/PyNucleus_nl/kernelsCy.pyx +++ b/nl/PyNucleus_nl/kernelsCy.pyx @@ -16,13 +16,18 @@ cimport numpy as np from PyNucleus_base.myTypes import REAL from PyNucleus_fem.functions cimport constant from . interactionDomains cimport ball1, ball2, ballInf, fullSpace -from . twoPointFunctions cimport productTwoPoint, inverseTwoPoint, productParametrizedTwoPoint +from . twoPointFunctions cimport (constantTwoPoint, + productTwoPoint, + inverseTwoPoint, + productParametrizedTwoPoint) from . fractionalOrders cimport (constFractionalOrder, variableFractionalOrder, - piecewiseConstantFractionalOrder, - constantFractionalLaplacianScaling, - variableFractionalLaplacianScaling, - variableFractionalLaplacianScalingBoundary) + piecewiseConstantFractionalOrder) +from . kernelNormalization cimport (constantFractionalLaplacianScaling, + constantFractionalLaplacianScalingDerivative, + variableFractionalLaplacianScaling, + variableFractionalLaplacianScalingBoundary, + variableFractionalLaplacianScalingWithDifferentHorizon) cdef inline REAL_t gammainc(REAL_t a, REAL_t x): @@ -259,6 +264,28 @@ cdef REAL_t indicatorKernel2D(REAL_t *x, REAL_t *y, void *c_params): return 0. +cdef REAL_t indicatorKernel1Dboundary(REAL_t *x, REAL_t *y, void *c_params): + cdef: + twoPointFunction interaction = (((c_params+fINTERACTION))[0]) + REAL_t C + if interaction.evalPtr(1, x, y) != 0.: + C = getREAL(c_params, fSCALING) + return -C*2.0*sqrt((x[0]-y[0])*(x[0]-y[0])) + else: + return 0. + + +cdef REAL_t indicatorKernel2Dboundary(REAL_t *x, REAL_t *y, void *c_params): + cdef: + twoPointFunction interaction = (((c_params+fINTERACTION))[0]) + REAL_t C + if interaction.evalPtr(2, x, y) != 0.: + C = getREAL(c_params, fSCALING) + return -C*sqrt((x[0]-y[0])*(x[0]-y[0])+(x[1]-y[1])*(x[1]-y[1])) + else: + return 0. + + cdef REAL_t peridynamicKernel1D(REAL_t *x, REAL_t *y, void *c_params): cdef: interactionDomain interaction = (((c_params+fINTERACTION))[0]) @@ -285,6 +312,29 @@ cdef REAL_t peridynamicKernel2D(REAL_t *x, REAL_t *y, void *c_params): return 0. +cdef REAL_t peridynamicKernel1Dboundary(REAL_t *x, REAL_t *y, void *c_params): + cdef: + interactionDomain interaction = (((c_params+fINTERACTION))[0]) + REAL_t C + REAL_t d + if interaction.evalPtr(1, x, y) != 0.: + d = abs(x[0]-y[0]) + C = getREAL(c_params, fSCALING) + return -2.0*C*log(d) + else: + return 0. + + +cdef REAL_t peridynamicKernel2Dboundary(REAL_t *x, REAL_t *y, void *c_params): + cdef: + interactionDomain interaction = (((c_params+fINTERACTION))[0]) + REAL_t C + if interaction.evalPtr(2, x, y) != 0.: + C = getREAL(c_params, fSCALING) + return -2.0*C + else: + return 0. + cdef REAL_t gaussianKernel1D(REAL_t *x, REAL_t *y, void *c_params): cdef: interactionDomain interaction = (((c_params+fINTERACTION))[0]) @@ -313,6 +363,34 @@ cdef REAL_t gaussianKernel2D(REAL_t *x, REAL_t *y, void *c_params): return 0. +cdef REAL_t gaussianKernel1Dboundary(REAL_t *x, REAL_t *y, void *c_params): + cdef: + interactionDomain interaction = (((c_params+fINTERACTION))[0]) + REAL_t C, invD + REAL_t d2 + if interaction.evalPtr(1, x, y) != 0.: + d2 = (x[0]-y[0])*(x[0]-y[0]) + C = getREAL(c_params, fSCALING) + invD = getREAL(c_params, fEXPONENTINVERSE) + return C*sqrt(1./(d2*invD))*gammainc(0.5, d2*invD)*sqrt(d2) + else: + return 0. + + +cdef REAL_t gaussianKernel2Dboundary(REAL_t *x, REAL_t *y, void *c_params): + cdef: + interactionDomain interaction = (((c_params+fINTERACTION))[0]) + REAL_t C, invD + REAL_t d2 + if interaction.evalPtr(2, x, y) != 0.: + d2 = (x[0]-y[0])*(x[0]-y[0]) + (x[1]-y[1])*(x[1]-y[1]) + C = getREAL(c_params, fSCALING) + invD = getREAL(c_params, fEXPONENTINVERSE) + return C*(1./(d2*invD))*gammainc(1.0, d2*invD)*sqrt(d2) + else: + return 0. + + cdef REAL_t updateAndEvalIntegrable(REAL_t *x, REAL_t *y, void *c_params): cdef: INDEX_t dim = getINDEX(c_params, fKDIM) @@ -365,14 +443,16 @@ cdef REAL_t updateAndEvalFractional(REAL_t *x, REAL_t *y, void *c_params): cdef class Kernel(twoPointFunction): """A kernel functions that can be used to define a nonlocal operator.""" - def __init__(self, INDEX_t dim, kernelType kType, function horizon, interactionDomain interaction, twoPointFunction scaling, twoPointFunction phi, BOOL_t piecewise=True): + def __init__(self, INDEX_t dim, kernelType kType, function horizon, interactionDomain interaction, twoPointFunction scaling, twoPointFunction phi, BOOL_t piecewise=True, BOOL_t boundary=False, INDEX_t vectorSize=1): cdef: parametrizedTwoPointFunction parametrizedScaling int i self.dim = dim + self.vectorSize = vectorSize self.kernelType = kType self.piecewise = piecewise + self.boundary = boundary self.c_kernel_params = malloc(NUM_KERNEL_PARAMS*OFFSET) for i in range(NUM_KERNEL_PARAMS): @@ -429,45 +509,85 @@ cdef class Kernel(twoPointFunction): self.variable = self.variableHorizon or self.variableScaling if self.piecewise: - if dim == 1: - if self.kernelType == INDICATOR: - self.kernelFun = indicatorKernel1D - elif self.kernelType == PERIDYNAMIC: - self.kernelFun = peridynamicKernel1D - elif self.kernelType == GAUSSIAN: - self.kernelFun = gaussianKernel1D - elif dim == 2: - if self.kernelType == INDICATOR: - self.kernelFun = indicatorKernel2D - elif self.kernelType == PERIDYNAMIC: - self.kernelFun = peridynamicKernel2D - elif self.kernelType == GAUSSIAN: - self.kernelFun = gaussianKernel2D - elif dim == 3: - pass + if not self.boundary: + if dim == 1: + if self.kernelType == INDICATOR: + self.kernelFun = indicatorKernel1D + elif self.kernelType == PERIDYNAMIC: + self.kernelFun = peridynamicKernel1D + elif self.kernelType == GAUSSIAN: + self.kernelFun = gaussianKernel1D + elif dim == 2: + if self.kernelType == INDICATOR: + self.kernelFun = indicatorKernel2D + elif self.kernelType == PERIDYNAMIC: + self.kernelFun = peridynamicKernel2D + elif self.kernelType == GAUSSIAN: + self.kernelFun = gaussianKernel2D + elif dim == 3: + pass + else: + raise NotImplementedError() else: - raise NotImplementedError() + if dim == 1: + if self.kernelType == INDICATOR: + self.kernelFun = indicatorKernel1Dboundary + elif self.kernelType == PERIDYNAMIC: + self.kernelFun = peridynamicKernel1Dboundary + elif self.kernelType == GAUSSIAN: + self.kernelFun = gaussianKernel1Dboundary + elif dim == 2: + if self.kernelType == INDICATOR: + self.kernelFun = indicatorKernel2Dboundary + elif self.kernelType == PERIDYNAMIC: + self.kernelFun = peridynamicKernel2Dboundary + elif self.kernelType == GAUSSIAN: + self.kernelFun = gaussianKernel2Dboundary + elif dim == 3: + pass + else: + raise NotImplementedError() else: self.kernelFun = updateAndEvalIntegrable - if dim == 1: - if self.kernelType == INDICATOR: - setFun(self.c_kernel_params, fEVAL, indicatorKernel1D) - elif self.kernelType == PERIDYNAMIC: - setFun(self.c_kernel_params, fEVAL, peridynamicKernel1D) - elif self.kernelType == GAUSSIAN: - setFun(self.c_kernel_params, fEVAL, gaussianKernel1D) - elif dim == 2: - if self.kernelType == INDICATOR: - setFun(self.c_kernel_params, fEVAL, indicatorKernel2D) - elif self.kernelType == PERIDYNAMIC: - setFun(self.c_kernel_params, fEVAL, peridynamicKernel2D) - elif self.kernelType == GAUSSIAN: - setFun(self.c_kernel_params, fEVAL, gaussianKernel2D) - elif dim == 3: - pass + if not self.boundary: + if dim == 1: + if self.kernelType == INDICATOR: + setFun(self.c_kernel_params, fEVAL, indicatorKernel1D) + elif self.kernelType == PERIDYNAMIC: + setFun(self.c_kernel_params, fEVAL, peridynamicKernel1D) + elif self.kernelType == GAUSSIAN: + setFun(self.c_kernel_params, fEVAL, gaussianKernel1D) + elif dim == 2: + if self.kernelType == INDICATOR: + setFun(self.c_kernel_params, fEVAL, indicatorKernel2D) + elif self.kernelType == PERIDYNAMIC: + setFun(self.c_kernel_params, fEVAL, peridynamicKernel2D) + elif self.kernelType == GAUSSIAN: + setFun(self.c_kernel_params, fEVAL, gaussianKernel2D) + elif dim == 3: + pass + else: + raise NotImplementedError() else: - raise NotImplementedError() + if dim == 1: + if self.kernelType == INDICATOR: + setFun(self.c_kernel_params, fEVAL, indicatorKernel1Dboundary) + elif self.kernelType == PERIDYNAMIC: + setFun(self.c_kernel_params, fEVAL, peridynamicKernel1Dboundary) + elif self.kernelType == GAUSSIAN: + setFun(self.c_kernel_params, fEVAL, gaussianKernel1Dboundary) + elif dim == 2: + if self.kernelType == INDICATOR: + setFun(self.c_kernel_params, fEVAL, indicatorKernel2Dboundary) + elif self.kernelType == PERIDYNAMIC: + setFun(self.c_kernel_params, fEVAL, peridynamicKernel2Dboundary) + elif self.kernelType == GAUSSIAN: + setFun(self.c_kernel_params, fEVAL, gaussianKernel2Dboundary) + elif dim == 3: + pass + else: + raise NotImplementedError() @property def singularityValue(self): @@ -561,18 +681,32 @@ cdef class Kernel(twoPointFunction): cdef REAL_t eval(self, REAL_t[::1] x, REAL_t[::1] y): return self.kernelFun(&x[0], &y[0], self.c_kernel_params) + cdef void evalVector(self, REAL_t[::1] x, REAL_t[::1] y, REAL_t[::1] vec): + vec[0] = self.kernelFun(&x[0], &y[0], self.c_kernel_params) + cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y): return self.kernelFun(x, y, self.c_kernel_params) + cdef void evalVectorPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y, INDEX_t vectorSize, REAL_t* vec): + vec[0] = self.kernelFun(x, y, self.c_kernel_params) + def __call__(self, REAL_t[::1] x, REAL_t[::1] y, BOOL_t callEvalParams=True): "Evaluate the kernel." if self.piecewise and callEvalParams: self.evalParams(x, y) return self.kernelFun(&x[0], &y[0], self.c_kernel_params) + def evalVector_py(self, REAL_t[::1] x, REAL_t[::1] y, REAL_t[::1] vec, BOOL_t callEvalParams=True): + "Evaluate the kernel." + if self.piecewise and callEvalParams: + self.evalParams(x, y) + self.evalVector(x, y, vec) + def getModifiedKernel(self, function horizon=None, twoPointFunction scaling=None): + cdef: + Kernel newKernel if horizon is None: horizon = self.horizon interaction = self.interaction @@ -584,6 +718,7 @@ cdef class Kernel(twoPointFunction): scaling = self.scaling from . kernels import getKernel newKernel = getKernel(dim=self.dim, kernel=self.kernelType, horizon=horizon, interaction=interaction, scaling=scaling, piecewise=self.piecewise) + setREAL(newKernel.c_kernel_params, fEXPONENTINVERSE, getREAL(self.c_kernel_params, fEXPONENTINVERSE)) return newKernel def getComplementKernel(self): @@ -602,13 +737,13 @@ cdef class Kernel(twoPointFunction): kernelName = 'Gaussian' else: raise NotImplementedError() - return "{}({}, {}, {})".format(self.__class__.__name__, kernelName, repr(self.interaction), self.scaling) + return "{}({}{}, {}, {})".format(self.__class__.__name__, kernelName, '' if not self.boundary else '-boundary', repr(self.interaction), self.scaling) def __getstate__(self): - return (self.dim, self.kernelType, self.horizon, self.interaction, self.scaling, self.phi, self.piecewise) + return (self.dim, self.kernelType, self.horizon, self.interaction, self.scaling, self.phi, self.piecewise, self.boundary) def __setstate__(self, state): - Kernel.__init__(self, state[0], state[1], state[2], state[3], state[4], state[5], state[6]) + Kernel.__init__(self, state[0], state[1], state[2], state[3], state[4], state[5], state[6], state[7]) def plot(self, x0=None): "Plot the kernel function." @@ -659,6 +794,30 @@ cdef class Kernel(twoPointFunction): plt.xlabel('$x_1-y_1$') plt.ylabel('$x_2-y_2$') + def getBoundaryKernel(self): + "Get the boundary kernel. This is the kernel that corresponds to the elimination of a subdomain via Gauss theorem." + cdef: + Kernel newKernel + from copy import deepcopy + + scaling = deepcopy(self.scaling) + if self.phi is not None: + phi = deepcopy(self.phi) + else: + phi = None + + from . kernels import getIntegrableKernel + newKernel = getIntegrableKernel(kernel=self.kernelType, + dim=self.dim, + horizon=deepcopy(self.horizon), + interaction=None, + scaling=scaling, + phi=phi, + piecewise=self.piecewise, + boundary=True) + setREAL(newKernel.c_kernel_params, fEXPONENTINVERSE, getREAL(self.c_kernel_params, fEXPONENTINVERSE)) + return newKernel + cdef class FractionalKernel(Kernel): @@ -677,10 +836,16 @@ cdef class FractionalKernel(Kernel): REAL_t tempered=0.): cdef: parametrizedTwoPointFunction parametrizedScaling - super(FractionalKernel, self).__init__(dim, FRACTIONAL, horizon, interaction, scaling, phi, piecewise) + if derivative == 0: + vectorSize = 1 + elif derivative == 1: + vectorSize = s.numParameters + else: + vectorSize = 1 + + super(FractionalKernel, self).__init__(dim, FRACTIONAL, horizon, interaction, scaling, phi, piecewise, boundary, vectorSize) self.symmetric = s.symmetric and isinstance(horizon, constant) and scaling.symmetric - self.boundary = boundary self.derivative = derivative self.temperedValue = tempered @@ -922,6 +1087,32 @@ cdef class FractionalKernel(Kernel): scalingValue = self.scaling.evalPtr(dim, x, y) self.setScalingValue(scalingValue) + cdef void evalVector(self, REAL_t[::1] x, REAL_t[::1] y, REAL_t[::1] vec): + cdef: + INDEX_t i + REAL_t fac + if self.derivative == 0: + vec[0] = self.kernelFun(&x[0], &y[0], self.c_kernel_params) + elif self.derivative == 1: + fac = self.kernelFun(&x[0], &y[0], self.c_kernel_params) + self.s.evalGrad(x, y, vec) + for i in range(self.vectorSize): + vec[i] *= fac + + cdef void evalVectorPtr(self, INDEX_t dim, REAL_t* x, REAL_t* y, INDEX_t vectorSize, REAL_t* vec): + cdef: + INDEX_t i + REAL_t fac + if self.derivative == 0: + vec[0] = self.kernelFun(x, y, self.c_kernel_params) + elif self.derivative == 1: + fac = self.kernelFun(x, y, self.c_kernel_params) + # print(fac) + self.s.evalGradPtr(dim, x, y, vectorSize, vec) + # print(vec[0]) + for i in range(vectorSize): + vec[i] *= fac + def getModifiedKernel(self, fractionalOrderBase s=None, function horizon=None, @@ -948,28 +1139,43 @@ cdef class FractionalKernel(Kernel): if scaling is None: scaling = self.scaling from . kernels import getFractionalKernel - newKernel = getFractionalKernel(dim=self.dim, s=s, horizon=horizon, interaction=interaction, scaling=scaling, piecewise=self.piecewise) + newKernel = getFractionalKernel(dim=self.dim, s=s, horizon=horizon, interaction=interaction, scaling=scaling, piecewise=self.piecewise, boundary=self.boundary, derivative=self.derivative) return newKernel def getBoundaryKernel(self): "Get the boundary kernel. This is the kernel that corresponds to the elimination of a subdomain via Gauss theorem." + cdef: + constantFractionalLaplacianScaling scalConst + constantFractionalLaplacianScalingDerivative scal + variableFractionalLaplacianScalingWithDifferentHorizon scalVarDiffHorizon + variableFractionalLaplacianScaling scalVar + from copy import deepcopy s = deepcopy(self.s) if not self.variableOrder: fac = constantTwoPoint(1/s.value) else: fac = inverseTwoPoint(s) - if self.phi is not None: - phi = fac*deepcopy(self.phi) + phi = fac + + if isinstance(self.scaling, constantFractionalLaplacianScalingDerivative): + scal = self.scaling + scaling = constantFractionalLaplacianScalingDerivative(scal.dim, scal.s, scal.horizon, scal.normalized, True, scal.derivative, scal.tempered) + elif isinstance(self.scaling, variableFractionalLaplacianScalingWithDifferentHorizon): + scalVarDiffHorizon = self.scaling + scaling = variableFractionalLaplacianScalingWithDifferentHorizon(scalVarDiffHorizon.symmetric, scalVarDiffHorizon.normalized, True, scalVarDiffHorizon.derivative, scalVarDiffHorizon.horizonFun) + elif isinstance(self.scaling, variableFractionalLaplacianScaling): + scalVar = self.scaling + scaling = variableFractionalLaplacianScaling(scalVar.symmetric, scalVar.normalized, True, scalVar.derivative) else: - phi = fac + scaling = deepcopy(self.scaling) from . kernels import getFractionalKernel newKernel = getFractionalKernel(dim=self.dim, s=s, horizon=deepcopy(self.horizon), interaction=None, - scaling=deepcopy(self.scaling), + scaling=scaling, phi=phi, piecewise=self.piecewise, boundary=True, @@ -981,6 +1187,22 @@ cdef class FractionalKernel(Kernel): newKernel = getFractionalKernel(dim=self.dim, s=self.s, horizon=self.horizon, interaction=self.interaction.getComplement(), scaling=self.scaling, piecewise=self.piecewise) return newKernel + def getDerivativeKernel(self): + cdef: + constantFractionalLaplacianScaling scal + variableFractionalLaplacianScaling scalVar + if isinstance(self.scaling, constantFractionalLaplacianScaling): + scal = self.scaling + scaling = constantFractionalLaplacianScalingDerivative(scal.dim, scal.s, scal.horizon, True, False, 1, scal.tempered) + elif isinstance(self.scaling, variableFractionalLaplacianScaling): + scalVar = self.scaling + scaling = variableFractionalLaplacianScaling(scalVar.symmetric, scalVar.normalized, scalVar.boundary, 1) + elif isinstance(self.scaling, constantTwoPoint): + scaling = constantFractionalLaplacianScalingDerivative(self.dim, self.sValue, np.nan, False, self.boundary, 1, 0.) + else: + raise NotImplementedError() + return FractionalKernel(self.dim, self.s, self.horizon, self.interaction, scaling, self.phi, False, self.boundary, 1, self.temperedValue) + def __repr__(self): if self.temperedValue != 0.: return "kernel(tempered-fractional{}, {}, {}, {}, {})".format('' if not self.boundary else '-boundary', self.s, repr(self.interaction), self.scaling, self.temperedValue) @@ -988,10 +1210,10 @@ cdef class FractionalKernel(Kernel): return "kernel(fractional{}, {}, {}, {})".format('' if not self.boundary else '-boundary', self.s, repr(self.interaction), self.scaling) def __getstate__(self): - return (self.dim, self.s, self.horizon, self.interaction, self.scaling, self.phi, self.piecewise) + return (self.dim, self.s, self.horizon, self.interaction, self.scaling, self.phi, self.piecewise, self.boundary, self.derivative, self.temperedValue) def __setstate__(self, state): - FractionalKernel.__init__(self, state[0], state[1], state[2], state[3], state[4], state[5], state[6]) + FractionalKernel.__init__(self, state[0], state[1], state[2], state[3], state[4], state[5], state[6], state[7], state[8], state[9]) cdef class RangedFractionalKernel(FractionalKernel): diff --git a/nl/PyNucleus_nl/nonlocalLaplacian.pyx b/nl/PyNucleus_nl/nonlocalLaplacian.pyx index 520a2b6..edbfc57 100644 --- a/nl/PyNucleus_nl/nonlocalLaplacian.pyx +++ b/nl/PyNucleus_nl/nonlocalLaplacian.pyx @@ -29,6 +29,8 @@ from PyNucleus_base.sparsityPattern cimport sparsityPattern from PyNucleus_base.linear_operators cimport (CSR_LinearOperator, SSS_LinearOperator, Dense_LinearOperator, + VectorLinearOperator, + Dense_VectorLinearOperator, Dense_SubBlock_LinearOperator, diagonalOperator, TimeStepperLinearOperator, @@ -42,8 +44,8 @@ from . fractionalOrders cimport (fractionalOrderBase, constFractionalOrder, piecewiseConstantFractionalOrder, variableFractionalOrder, - singleVariableUnsymmetricFractionalOrder, - variableFractionalLaplacianScaling) + singleVariableUnsymmetricFractionalOrder) +from . kernelNormalization cimport variableFractionalLaplacianScaling from . kernels import getFractionalKernel from . clusterMethodCy import (assembleFarFieldInteractions, getDoFBoxesAndCells, @@ -618,6 +620,98 @@ cdef class IndexManager: return s +cdef class IndexManagerVector(IndexManager): + cdef: + VectorLinearOperator vecA + INDEX_t vectorSize + + def __init__(self, DoFMap dm, VectorLinearOperator A=None, cellPairIdentifierSize=1, indexSet myDofs=None, sparsityPattern sP=None): + super(IndexManagerVector, self).__init__(dm, None, cellPairIdentifierSize, myDofs, sP) + self.vecA = A + self.vectorSize = A.vectorSize + + cdef inline void addToMatrixElemSymVector(self, REAL_t[:, ::1] contrib, REAL_t fac): + cdef: + INDEX_t k, p, q, I, J + for p in range(contrib.shape[0]): + for q in range(self.vectorSize): + contrib[p, q] *= fac + k = 0 + for p in range(self.dm.dofs_per_element): + I = self.localDoFs[p] + if I >= 0: + self.vecA.addToEntry(I, I, contrib[k, :]) + k += 1 + for q in range(p+1, self.dm.dofs_per_element): + J = self.localDoFs[q] + if J >= 0: + self.vecA.addToEntry(I, J, contrib[k, :]) + self.vecA.addToEntry(J, I, contrib[k, :]) + k += 1 + else: + k += self.dm.dofs_per_element-p + + cdef inline void addToMatrixElemVector(self, REAL_t[:, ::1] contrib, REAL_t fac): + cdef: + INDEX_t k, p, q, I, J + for p in range(contrib.shape[0]): + for q in range(self.vectorSize): + contrib[p, q] *= fac + k = 0 + for p in range(self.dm.dofs_per_element): + I = self.localDoFs[p] + if I >= 0: + for q in range(self.dm.dofs_per_element): + J = self.localDoFs[q] + if J >= 0: + self.vecA.addToEntry(I, J, contrib[k, :]) + k += 1 + else: + k += self.dm.dofs_per_element + + cdef inline void addToMatrixElemElemSymVector(self, REAL_t[:, ::1] contrib, REAL_t fac): + # Add symmetric 'contrib' to elements i and j in symmetric fashion + cdef: + INDEX_t k, p, q, I, J + for p in range(contrib.shape[0]): + for q in range(self.vectorSize): + contrib[p, q] *= fac + k = 0 + for p in range(2*self.dm.dofs_per_element): + I = self.localDoFs[p] + if I >= 0: + self.vecA.addToEntry(I, I, contrib[k, :]) + k += 1 + for q in range(p+1, 2*self.dm.dofs_per_element): + J = self.localDoFs[q] + if J >= 0: + self.vecA.addToEntry(I, J, contrib[k, :]) + self.vecA.addToEntry(J, I, contrib[k, :]) + k += 1 + else: + k += 2*self.dm.dofs_per_element-p + + cdef inline void addToMatrixElemElemVector(self, REAL_t[:, ::1] contrib, REAL_t fac): + # Add general 'contrib' to elements i and j + cdef: + INDEX_t k, p, q, I, J + for p in range(contrib.shape[0]): + for q in range(self.vectorSize): + contrib[p, q] *= fac + k = 0 + for p in range(2*self.dm.dofs_per_element): + I = self.localDoFs[p] + if I >= 0: + for q in range(2*self.dm.dofs_per_element): + J = self.localDoFs[q] + if J >= 0: + self.vecA.addToEntry(I, J, contrib[k, :]) + k += 1 + else: + k += 2*self.dm.dofs_per_element + + + # These functions are used by getEntry cdef inline MASK_t getElemSymMask(DoFMap DoFMap, INDEX_t cellNo1, INDEX_t I, INDEX_t J): @@ -758,6 +852,7 @@ cdef class nonlocalBuilder: else: self.contribZeroExterior = uninitialized((0), dtype=REAL) + if PLogger is not None: self.PLogger = PLogger else: @@ -851,7 +946,30 @@ cdef class nonlocalBuilder: else: raise NotImplementedError() else: - local_matrix = None + assert isinstance(self.kernel.horizon, constant) + if infHorizon: + kernelInfHorizon = self.kernel.getModifiedKernel(horizon=constant(np.inf)) + else: + kernelInfHorizon = self.kernel + if quadType == 'classical-refactored': + kernelBoundary = kernelInfHorizon.getBoundaryKernel() + if self.mesh.dim == 1: + local_matrix = fractionalLaplacian1D_boundary(kernelBoundary, + mesh=self.mesh, + DoFMap=self.dm, + target_order=target_order) + elif self.mesh.dim == 2: + if not isinstance(self.dm, Product_DoFMap): + local_matrix = fractionalLaplacian2D_boundary(kernelBoundary, + mesh=self.mesh, + DoFMap=self.dm, + target_order=target_order) + else: + raise NotImplementedError() + else: + raise NotImplementedError() + else: + local_matrix = None return local_matrix def getSparse(self, BOOL_t returnNearField=False, str prefix=''): @@ -1228,6 +1346,113 @@ cdef class nonlocalBuilder: LOGGER.warning('Not converting dense to sparse matrix, since only {}% of entries are zero.'.format(100.*ratio)) return A + def getDenseVector(self, BOOL_t trySparsification=False): + cdef: + INDEX_t cellNo1, cellNo2 + VectorLinearOperator A = None + REAL_t[:, ::1] contrib = self.contribVector, contribZeroExterior = self.contribZeroExteriorVector + INDEX_t start, end + meshBase surface + IndexManagerVector iM + INDEX_t i, j, explicitZerosRow + np.int64_t explicitZeros + REAL_t[:, ::1] data + REAL_t sparsificationThreshold = 0.8 + BOOL_t symmetricLocalMatrix = self.local_matrix.symmetricLocalMatrix + BOOL_t symmetricCells = self.local_matrix.symmetricCells + MASK_t mask + + if self.comm: + start = np.ceil(self.mesh.num_cells*self.comm.rank/self.comm.size) + end = np.ceil(self.mesh.num_cells*(self.comm.rank+1)/self.comm.size) + else: + start = 0 + end = self.mesh.num_cells + + + if self.dm2 is None: + A = Dense_VectorLinearOperator(np.zeros((self.dm.num_dofs, self.dm.num_dofs, self.kernel.vectorSize), dtype=REAL)) + else: + A = Dense_VectorLinearOperator(np.zeros((self.dm.num_dofs, self.dm2.num_dofs, self.kernel.vectorSize), dtype=REAL)) + + if self.dm2 is None: + iM = IndexManagerVector(self.dm, A) + else: + LOGGER.warning('Efficiency of assembly with 2 DoFMaps is bad.') + dmCombined = self.dm.combine(self.dm2) + B = SubMatrixAssemblyOperator(A, + np.arange(self.dm.num_dofs, dtype=INDEX), + np.arange(self.dm.num_dofs, self.dm.num_dofs+self.dm2.num_dofs, dtype=INDEX)) + iM = IndexManagerVector(dmCombined, B) + + # Omega x Omega + with self.PLogger.Timer('interior'): + for cellNo1 in range(start, end): + self.local_matrix.setCell1(cellNo1) + for cellNo2 in range(cellNo1, self.mesh.num_cells): + self.local_matrix.setCell2(cellNo2) + if iM.getDoFsElemElem(cellNo1, cellNo2): + continue + panel = self.local_matrix.getPanelType() + if cellNo1 == cellNo2: + if panel != IGNORED: + self.local_matrix.evalVector(contrib, panel) + if symmetricLocalMatrix: + iM.addToMatrixElemElemSymVector(contrib, 1.) + else: + iM.addToMatrixElemElemVector(contrib, 1.) + else: + if symmetricCells: + if panel != IGNORED: + self.local_matrix.evalVector(contrib, panel) + # If the kernel is symmetric, the contributions from (cellNo1, cellNo2) and (cellNo2, cellNo1) + # are the same. We multiply by 2 to account for the contribution from cells (cellNo2, cellNo1). + if symmetricLocalMatrix: + iM.addToMatrixElemElemSymVector(contrib, 2.) + else: + iM.addToMatrixElemElemVector(contrib, 2.) + else: + if panel != IGNORED: + self.local_matrix.evalVector(contrib, panel) + if symmetricLocalMatrix: + iM.addToMatrixElemElemSymVector(contrib, 1.) + else: + iM.addToMatrixElemElemVector(contrib, 1.) + self.local_matrix.swapCells() + panel = self.local_matrix.getPanelType() + if panel != IGNORED: + if iM.getDoFsElemElem(cellNo2, cellNo1): + continue + self.local_matrix.evalVector(contrib, panel) + if symmetricLocalMatrix: + iM.addToMatrixElemElemSymVector(contrib, 1.) + else: + iM.addToMatrixElemElemVector(contrib, 1.) + self.local_matrix.swapCells() + + # Omega x Omega^C + if self.zeroExterior: + with self.PLogger.Timer('zeroExterior'): + surface = self.mesh.get_surface_mesh() + + self.local_matrix_zeroExterior.setMesh2(surface) + + for cellNo1 in range(start, end): + iM.getDoFsElem(cellNo1) + mask = iM.getElemSymMask() + self.local_matrix_zeroExterior.setCell1(cellNo1) + for cellNo2 in range(surface.num_cells): + self.local_matrix_zeroExterior.setCell2(cellNo2) + panel = self.local_matrix_zeroExterior.getPanelType() + self.local_matrix_zeroExterior.evalVector(contribZeroExterior, panel, mask) + # if local_matrix_zeroExterior.symmetricLocalMatrix: + iM.addToMatrixElemSymVector(contribZeroExterior, 1.) + # else: + # raise NotImplementedError() + if self.comm: + self.comm.Allreduce(MPI.IN_PLACE, A.data) + return A + cpdef REAL_t getEntryCluster(self, INDEX_t I, INDEX_t J): cdef: tree_node n1, n2, n3 @@ -1627,7 +1852,10 @@ cdef class nonlocalBuilder: vol = 2*np.pi * self.kernel.horizonValue else: raise NotImplementedError() - coeff = constant(-vol*self.kernel.scalingValue*pow(self.kernel.horizonValue, 1-self.mesh.dim-2*self.kernel.sValue)/self.kernel.sValue) + x = np.zeros((self.mesh.dim), dtype=REAL) + y = np.zeros((self.mesh.dim), dtype=REAL) + y[0] = self.kernel.horizonValue + coeff = constant(-vol*self.local_matrix_zeroExterior.kernel(x, y)) qr = simplexXiaoGimbutas(2, self.mesh.dim) if self.mesh.dim == 1: mass = mass_1d_sym_scalar_anisotropic(coeff, self.dm, qr) @@ -2447,7 +2675,6 @@ cdef class nonlocalBuilder: cdef: meshBase mesh = self.mesh DoFMap DoFMap = self.dm - fractionalOrderBase s = self.kernel.s REAL_t[:, :, ::1] boxes = None, local_boxes sparseGraph cells = None REAL_t[:, ::1] coords = None @@ -2457,17 +2684,9 @@ cdef class nonlocalBuilder: BOOL_t forceUnsymmetric, doDistributedAssembly = False, assembleOnRoot = True, localFarFieldIndexing = False refinementParams refParams CSR_LinearOperator lclR - assert isinstance(self.kernel, FractionalKernel), 'H2 is only implemented for fractional kernels' refParams = self.getH2RefinementParams() - LOGGER.info('interpolation_order: {}, maxLevels: {}, minClusterSize: {}, minMixedClusterSize: {}, minFarFieldBlockSize: {}, eta: {}'.format(refParams.interpolation_order, - refParams.maxLevels, - refParams.minSize, - refParams.minMixedSize, - refParams.farFieldInteractionSize, - refParams.eta)) - forceUnsymmetric = self.params.get('forceUnsymmetric', False) doDistributedAssembly = self.comm is not None and self.comm.size > 1 and DoFMap.num_dofs > self.comm.size assembleOnRoot = self.params.get('assembleOnRoot', True) @@ -2491,6 +2710,13 @@ cdef class nonlocalBuilder: lenPfar = self.comm.bcast(lenPfar) if lenPfar > 0: + LOGGER.info('interpolation_order: {}, maxLevels: {}, minClusterSize: {}, minMixedClusterSize: {}, minFarFieldBlockSize: {}, eta: {}'.format(refParams.interpolation_order, + refParams.maxLevels, + refParams.minSize, + refParams.minMixedSize, + refParams.farFieldInteractionSize, + refParams.eta)) + # get near field matrix with self.PLogger.Timer('near field'): Anear = self.assembleClusters(Pnear, jumps=jumps, forceUnsymmetric=forceUnsymmetric, myRoot=myRoot, doDistributedAssembly=doDistributedAssembly) @@ -2503,12 +2729,12 @@ cdef class nonlocalBuilder: with self.PLogger.Timer('leaf values'): # get leave values - if (s.min > 0) and (s.max < 1): + if self.kernel.max_singularity > -self.kernel.dim-2: if not localFarFieldIndexing: root.enterLeafValues(mesh, DoFMap, refParams.interpolation_order, boxes, self.comm, assembleOnRoot=assembleOnRoot) else: myRoot.enterLeafValues(mesh, local_dm, refParams.interpolation_order, local_boxes, local=True) - elif (s.min > 1) and (s.max < 2): + elif (self.kernel.min_singularity < -self.kernel.dim-2) and (self.kernel.max_singularity > -self.kernel.dim-4): if not localFarFieldIndexing: root.enterLeafValuesGrad(mesh, DoFMap, refParams.interpolation_order, boxes, self.comm) else: @@ -2536,14 +2762,13 @@ cdef class nonlocalBuilder: h2 = DistributedH2Matrix_localData(local_h2, Pnear, self.comm, self.dm, local_dm, lclR, lclP) else: h2 = nullOperator(self.dm.num_dofs, self.dm.num_dofs) - + LOGGER.info('{}'.format(h2)) elif len(Pnear) == 0: h2 = nullOperator(self.dm.num_dofs, self.dm.num_dofs) else: - LOGGER.info('Near field pairs: {}, far field pairs: {}, tree nodes: {}'.format(len(Pnear), lenPfar, root.nodes)) + LOGGER.info('Cannot assemble H2 operator, assembling dense matrix instead') with self.PLogger.Timer('dense operator'): h2 = self.getDense() - LOGGER.info('{}'.format(h2)) if returnNearField: if returnTree: return h2, Pnear, root diff --git a/nl/PyNucleus_nl/nonlocalLaplacianBase.pxd b/nl/PyNucleus_nl/nonlocalLaplacianBase.pxd index 03fda04..3c14fca 100644 --- a/nl/PyNucleus_nl/nonlocalLaplacianBase.pxd +++ b/nl/PyNucleus_nl/nonlocalLaplacianBase.pxd @@ -87,6 +87,10 @@ cdef class double_local_matrix_t: REAL_t[::1] contrib, panelType panel, MASK_t mask=*) + cdef void evalVector(self, + REAL_t[:, ::1] contrib, + panelType panel, + MASK_t mask=*) cdef panelType getQuadOrder(self, const REAL_t h1, const REAL_t h2, diff --git a/nl/PyNucleus_nl/nonlocalLaplacianBase.pyx b/nl/PyNucleus_nl/nonlocalLaplacianBase.pyx index 66cfa04..f78052a 100644 --- a/nl/PyNucleus_nl/nonlocalLaplacianBase.pyx +++ b/nl/PyNucleus_nl/nonlocalLaplacianBase.pyx @@ -309,11 +309,22 @@ cdef class double_local_matrix_t: MASK_t mask=ALL): raise NotImplementedError() + cdef void evalVector(self, + REAL_t[:, ::1] contrib, + panelType panel, + MASK_t mask=ALL): + raise NotImplementedError() + def eval_py(self, REAL_t[::1] contrib, panel): self.eval(contrib, panel, ALL) + def evalVector_py(self, + REAL_t[:, ::1] contrib, + panel): + self.evalVector(contrib, panel, ALL) + cdef panelType getQuadOrder(self, const REAL_t h1, const REAL_t h2, @@ -901,6 +912,14 @@ cdef class nonlocalLaplacian(double_local_matrix_t): INDEX_t dim = simplex1.shape[1] BOOL_t cutElements = False REAL_t w + REAL_t c1, c2, val2, PHI_I_0, PHI_I_1, PSI_J + transformQuadratureRule qr0trans, qr1trans + INDEX_t dofs_per_element, numQuadNodes0, numQuadNodes1 + REAL_t a_b1[3] + REAL_t a_A1[3][3] + REAL_t a_A2[3][3] + REAL_t[::1] b1 + REAL_t[:, ::1] A1, A2 if self.kernel.finiteHorizon: # check if the horizon might cut the elements @@ -948,7 +967,63 @@ cdef class nonlocalLaplacian(double_local_matrix_t): contrib[k] = val*vol k += 1 else: - raise NotImplementedError() + if panel < 0: + sQR = (self.distantQuadRulesPtr[MAX_PANEL+panel]) + else: + sQR = (self.distantQuadRulesPtr[panel]) + qr2 = (sQR.qr) + if sQR.qrTransformed0 is not None: + qr0trans = sQR.qrTransformed0 + else: + qr0 = qr2.rule1 + qr0trans = transformQuadratureRule(qr0) + sQR.qrTransformed0 = qr0trans + if sQR.qrTransformed1 is not None: + qr1trans = sQR.qrTransformed1 + else: + qr1 = qr2.rule2 + qr1trans = transformQuadratureRule(qr1) + sQR.qrTransformed1 = qr1trans + + numQuadNodes0 = qr0trans.num_nodes + numQuadNodes1 = qr1trans.num_nodes + + vol = vol1*vol2 + dofs_per_element = self.DoFMap.dofs_per_element + + A1 = a_A1 + b1 = a_b1 + A2 = a_A2 + + self.kernel.interaction.startLoopSubSimplices_Simplex(simplex1, simplex2) + while self.kernel.interaction.nextSubSimplex_Simplex(A1, b1, &c1): + qr0trans.setAffineBaryTransform(A1, b1) + qr0trans.nodesInGlobalCoords(simplex1, self.x) + for i in range(qr0trans.num_nodes): + self.kernel.interaction.startLoopSubSimplices_Node(self.x[i, :], simplex2) + while self.kernel.interaction.nextSubSimplex_Node(A2, &c2): + qr1trans.setLinearBaryTransform(A2) + qr1trans.nodesInGlobalCoords(simplex2, self.y) + for j in range(qr1trans.num_nodes): + w = qr0trans.weights[i]*qr1trans.weights[j]*c1 * c2 * vol + val = w*self.kernel.evalPtr(dim, &self.x[i, 0], &self.y[j, 0]) + val2 = w*self.kernel.evalPtr(dim, &self.y[j, 0], &self.x[i, 0]) + k = 0 + for I in range(2*dofs_per_element): + if I < dofs_per_element: + PHI_I_0 = self.getLocalShapeFunction(I).evalStrided(&qr0trans.nodes[0, i], numQuadNodes0) + PHI_I_1 = 0. + else: + PHI_I_0 = 0. + PHI_I_1 = self.getLocalShapeFunction(I-dofs_per_element).evalStrided(&qr1trans.nodes[0, j], numQuadNodes1) + for J in range(2*dofs_per_element): + if mask[k]: + if J < dofs_per_element: + PSI_J = self.getLocalShapeFunction(J).evalStrided(&qr0trans.nodes[0, i], numQuadNodes0) + else: + PSI_J = -self.getLocalShapeFunction(J-dofs_per_element).evalStrided(&qr1trans.nodes[0, j], numQuadNodes1) + contrib[k] += (val * PHI_I_0 - val2 * PHI_I_1) * PSI_J + k += 1 cdef void addQuadRule_boundary(self, panelType panel): cdef: diff --git a/nl/PyNucleus_nl/nonlocalProblems.py b/nl/PyNucleus_nl/nonlocalProblems.py index f293d17..55d83c1 100644 --- a/nl/PyNucleus_nl/nonlocalProblems.py +++ b/nl/PyNucleus_nl/nonlocalProblems.py @@ -29,7 +29,7 @@ from PyNucleus_fem import (PHYSICAL, NO_BOUNDARY, DIRICHLET, HOMOGENEOUS_DIRICHLET, NEUMANN, HOMOGENEOUS_NEUMANN, - NORM) + NORM, dofmapFactory) from PyNucleus_fem.factories import functionFactory, rhsFractional2D_nonPeriodic from scipy.special import gamma as Gamma, binom from . twoPointFunctions import (constantTwoPoint, @@ -52,18 +52,35 @@ smoothedInnerOuterFractionalOrder, islandsFractionalOrder, layersFractionalOrder, - singleVariableUnsymmetricFractionalOrder) + singleVariableUnsymmetricFractionalOrder, + feFractionalOrder) from . kernelsCy import (getKernelEnum, FRACTIONAL, INDICATOR, PERIDYNAMIC, GAUSSIAN, ) from . kernels import (getFractionalKernel, getIntegrableKernel, getKernel) +from copy import deepcopy + + +class fractionalOrderFactoryClass(factory): + def build(self, name, *args, **kwargs): + dm = None + if 'dm' in kwargs: + dm = kwargs.pop('dm') + if dm is not None: + s = self.build(name, *args, **kwargs) + assert isinstance(s, (constFractionalOrder, variableConstFractionalOrder, + constantNonSymFractionalOrder, singleVariableUnsymmetricFractionalOrder)) + sVec = dm.interpolate(s.fixedY(np.zeros((dm.mesh.dim), dtype=REAL))) + return super().build('fe', sVec, s.min, s.max) + else: + return super().build(name, *args, **kwargs) -fractionalOrderFactory = factory() +fractionalOrderFactory = fractionalOrderFactoryClass() fractionalOrderFactory.register('constant', constFractionalOrder, aliases=['const']) -fractionalOrderFactory.register('varConst', variableConstFractionalOrder, aliases=['constVar']) +fractionalOrderFactory.register('varConst', variableConstFractionalOrder, aliases=['constVar', 'constantSym']) fractionalOrderFactory.register('leftRight', leftRightFractionalOrder, aliases=['twoDomain']) fractionalOrderFactory.register('linearLeftRightNonSym', linearLeftRightFractionalOrder) fractionalOrderFactory.register('smoothedLeftRight', smoothedLeftRightFractionalOrder, params={'r': 0.1, 'slope': 200.}, aliases=['twoDomainNonSym']) @@ -72,6 +89,7 @@ fractionalOrderFactory.register('innerOuterNonSym', smoothedInnerOuterFractionalOrder) fractionalOrderFactory.register('islands', islandsFractionalOrder, params={'r': 0.1, 'r2': 0.6}) fractionalOrderFactory.register('layers', layersFractionalOrder) +fractionalOrderFactory.register('fe', feFractionalOrder) twoPointFunctionFactory = factory() twoPointFunctionFactory.register('constant', constantTwoPoint, aliases=['const', 'constantTwoPoint']) @@ -110,25 +128,31 @@ def build(self, name, kernel, boundaryCondition, noRef=0, useMulti=False, **kwar if 'skipMesh' in kwargs: skipMesh = kwargs.pop('skipMesh') + if kernel is None: + horizonValue = 0. + else: + horizonValue = kernel.horizon.value + domainIndicator, boundaryIndicator, interactionIndicator = super(nonlocalMeshFactoryClass, self).build(name, **kwargs) if boundaryCondition == HOMOGENEOUS_DIRICHLET: - if kernel.horizon.value == np.inf: - if kernel.s.max < 0.5: - tag = NO_BOUNDARY - else: - tag = PHYSICAL + if horizonValue == np.inf: + # if kernel.s.max < 0.5: + # tag = NO_BOUNDARY + # else: + # tag = PHYSICAL + tag = PHYSICAL zeroExterior = True else: tag = domainIndicator zeroExterior = False - hasInteractionDomain = kernel.horizon.value < np.inf + hasInteractionDomain = 0 < horizonValue < np.inf elif boundaryCondition == HOMOGENEOUS_NEUMANN: tag = NO_BOUNDARY zeroExterior = False hasInteractionDomain = False elif boundaryCondition == DIRICHLET: - if kernel.horizon.value == np.inf: + if horizonValue == np.inf: if kernel.s.max < 0.5: tag = NO_BOUNDARY else: @@ -137,9 +161,9 @@ def build(self, name, kernel, boundaryCondition, noRef=0, useMulti=False, **kwar else: tag = NO_BOUNDARY zeroExterior = False - hasInteractionDomain = True + hasInteractionDomain = 0 < horizonValue < np.inf elif boundaryCondition == NEUMANN: - if kernel.horizon.value == np.inf: + if horizonValue == np.inf: assert False else: tag = NO_BOUNDARY @@ -154,8 +178,8 @@ def build(self, name, kernel, boundaryCondition, noRef=0, useMulti=False, **kwar if not skipMesh: if hasInteractionDomain: - assert 0 < kernel.horizon.value < np.inf - kwargs['horizon'] = kernel.horizon.value + assert 0 < horizonValue < np.inf, horizonValue + kwargs['horizon'] = horizonValue mesh = self.overlappingMeshFactory.build(name, noRef, **kwargs) else: mesh = self.nonOverlappingMeshFactory.build(name, noRef, **kwargs) @@ -260,6 +284,7 @@ def __init__(self, driver): self.addProperty('phiArgs') self.addProperty('admissibleParams') self.admissibleParams = None + self.feFractionalOrder = None def setDriverArgs(self): p = self.driver.addGroup('kernel') @@ -288,6 +313,7 @@ def setDriverArgs(self): self.setDriverFlag('phi', 'const(1.)', argInterpreter=self.argInterpreter(['const', 'twoDomain', 'twoDomainNonSym', 'tempered']), help='kernel coefficient', group=p) self.setDriverFlag('normalized', True, help='kernel normalization', group=p) + self.setDriverFlag('discretizedOrder', False, help='Use a FE function for the fractional order s.') def processCmdline(self, params): dim = nonlocalMeshFactory.getDim(params['domain']) @@ -349,8 +375,13 @@ def processCmdline(self, params): def getDim(self, domain): self.dim = nonlocalMeshFactory.getDim(domain) + @generates('dmAux') + def constructAuxiliarySpace(self): + self.dmAux = None + @generates(['kernel', 'rangedKernel']) - def processKernel(self, dim, kernelType, sType, sArgs, phiType, phiArgs, horizon, interaction, normalized, admissibleParams): + def processKernel(self, dim, kernelType, sType, sArgs, phiType, phiArgs, horizon, interaction, normalized, admissibleParams, + discretizedOrder, dmAux, feFractionalOrder): if kernelType == 'local': self.kernel = None @@ -382,11 +413,26 @@ def processKernel(self, dim, kernelType, sType, sArgs, phiType, phiArgs, horizon self.rangedKernel = None if kType == FRACTIONAL: - try: - sFun = fractionalOrderFactory(sType, *sArgs) - except TypeError: - sArgs = (sArgs, ) - sFun = fractionalOrderFactory(sType, *sArgs) + if feFractionalOrder is None: + if isinstance(sArgs, dict): + if discretizedOrder: + sFun = fractionalOrderFactory(sType, dm=dmAux, **sArgs) + else: + sFun = fractionalOrderFactory(sType, **sArgs) + else: + try: + if discretizedOrder: + sFun = fractionalOrderFactory(sType, *sArgs, dm=dmAux) + else: + sFun = fractionalOrderFactory(sType, *sArgs) + except TypeError: + sArgs = (sArgs, ) + if discretizedOrder: + sFun = fractionalOrderFactory(sType, *sArgs, dm=dmAux) + else: + sFun = fractionalOrderFactory(sType, *sArgs) + else: + sFun = deepcopy(feFractionalOrder) else: sFun = None @@ -491,12 +537,35 @@ def processCmdline(self, params): params['noRef'] = noRef super().processCmdline(params) + @generates('domainParams') + def getDomainParams(self, domain): + meshParams = {} + if domain == 'interval': + radius = 1. + meshParams.update({'a': -radius, 'b': radius}) + elif domain == 'disconnectedInterval': + meshParams['sep'] = 0.1 + elif domain == 'disc': + radius = 1. + meshParams.update({'h': 0.78, 'radius': radius}) + elif domain == 'square': + meshParams.update({'N': 3, 'ax': -1, 'ay': -1, 'bx': 1, 'by': 1}) + elif domain == 'Lshape': + pass + elif domain == 'cutoutCircle': + meshParams.update({'radius': 1., 'cutoutAngle': np.pi/2.}) + elif domain == 'ball': + pass + else: + raise NotImplementedError(domain) + self.domainParams = meshParams + @generates(['analyticSolution', 'exactHsSquared', 'exactL2Squared', 'rhs', 'mesh_domain', 'mesh_params', 'tag', 'boundaryCondition', 'domainIndicator', 'interactionIndicator', 'fluxIndicator', 'zeroExterior', 'rhsData', 'dirichletData', 'fluxData']) - def processProblem(self, kernel, dim, domain, problem, normalized): + def processProblem(self, kernel, dim, domain, domainParams, problem, normalized): s = kernel.s self.analyticSolution = None self.exactHsSquared = None @@ -505,10 +574,8 @@ def processProblem(self, kernel, dim, domain, problem, normalized): assert normalized boundaryCondition = HOMOGENEOUS_DIRICHLET - meshParams = {'kernel': kernel} if domain == 'interval': radius = 1. - meshParams.update({'a': -radius, 'b': radius}) if problem == 'constant': self.rhs = constant(1.) @@ -580,15 +647,12 @@ def fun(x): else: raise NotImplementedError(problem) elif domain == 'disconnectedInterval': - meshParams['sep'] = 0.1 - if problem == 'constant': self.rhs = Lambda(lambda x: 1. if x[0] > 0.5 else 0.) else: raise NotImplementedError(problem) elif domain == 'disc': radius = 1. - meshParams.update({'h': 0.78, 'radius': radius}) if problem == 'constant': self.rhs = constant(1.) @@ -645,8 +709,6 @@ def fun(x): else: raise NotImplementedError(problem) elif domain == 'square': - meshParams.update({'N': 3, 'ax': -1, 'ay': -1, 'bx': 1, 'by': 1}) - if problem == 'constant': self.rhs = constant(1.) elif problem == 'sin': @@ -664,8 +726,6 @@ def fun(x): else: raise NotImplementedError(problem) elif domain == 'cutoutCircle': - meshParams.update({'radius': 1., 'cutoutAngle': np.pi/2.}) - if problem == 'constant': self.rhs = constant(1.) elif problem == 'sin': @@ -687,6 +747,8 @@ def fun(x): raise NotImplementedError(domain) mesh_domain = domain + meshParams = {'kernel': kernel} + meshParams.update(domainParams) self.boundaryCondition = meshParams['boundaryCondition'] = boundaryCondition meshParams['useMulti'] = self.useMulti self.mesh_domain = mesh_domain @@ -695,7 +757,10 @@ def fun(x): self.tag = nI['tag'] self.domainIndicator = nI['domain'] self.interactionIndicator = nI['interaction']+nI['boundary'] - self.fluxIndicator = functionFactory('constant', 0.) + if boundaryCondition in (NEUMANN, HOMOGENEOUS_NEUMANN): + self.fluxIndicator = self.interactionIndicator + else: + self.fluxIndicator = functionFactory('constant', 0.) self.zeroExterior = nI['zeroExterior'] self.dirichletData = None self.fluxData = None @@ -725,6 +790,19 @@ def getApproximationParams(self, dim, kernel, element): def buildMesh(self, mesh_domain, mesh_params): self.mesh, _ = nonlocalMeshFactory.build(mesh_domain, **mesh_params) + @generates('dmAux') + def constructAuxiliarySpace(self, dim, domain, domainParams, kernelType, horizon): + # This is not the actual kernel that we use. + # We just need something to get a mesh to support the fractional order. + kType = getKernelEnum(kernelType) + if kType == FRACTIONAL: + sFun = fractionalOrderFactory('const', 0.75) + kernel = getKernel(dim=dim, kernel=kType, s=sFun, horizon=horizon) + else: + kernel = getKernel(dim=dim, kernel=kType, horizon=horizon) + mesh, _ = nonlocalMeshFactory(domain, kernel=kernel, boundaryCondition=HOMOGENEOUS_DIRICHLET, **domainParams) + self.dmAux = dofmapFactory('P1', mesh, NO_BOUNDARY) + def getIdentifier(self, params): keys = ['domain', 'problem', 's', 'noRef', 'element', 'adaptive'] d = [] @@ -741,10 +819,11 @@ def setDriverArgs(self): super().setDriverArgs() self.setDriverFlag('domain', 'interval', acceptedValues=['gradedInterval', 'square', 'disc', 'gradedDisc', 'discWithIslands'], help='spatial domain') self.addParametrizedArg('indicator', [float, float]) + self.addParametrizedArg('polynomial', [int]) self.setDriverFlag('problem', 'poly-Dirichlet', - argInterpreter=self.argInterpreter(['indicator'], acceptedValues=['poly-Dirichlet', 'poly-Dirichlet2', 'poly-Dirichlet3', - 'poly-Neumann', 'zeroFlux', 'source', 'constant', - 'exact-sin-Dirichlet', 'exact-sin-Neumann', 'discontinuous']), + argInterpreter=self.argInterpreter(['indicator', 'polynomial'], acceptedValues=['poly-Dirichlet', 'polynomial', 'poly-Dirichlet2', 'poly-Dirichlet3', + 'poly-Neumann', 'zeroFlux', 'source', 'constant', + 'exact-sin-Dirichlet', 'exact-sin-Neumann', 'discontinuous']), help="select a problem to solve") self.setDriverFlag('noRef', argInterpreter=int) self.setDriverFlag('element', acceptedValues=['P1', 'P0', 'P2'], help="finite element space") @@ -811,6 +890,44 @@ def processProblem(self, kernel, domain, problem, normalized): self.dirichletData = Lambda(lambda x: 1-x[0]**2) if ((kType == FRACTIONAL and isinstance(sFun, constFractionalOrder)) or kType in (INDICATOR, PERIDYNAMIC, GAUSSIAN)) and phiFun is None and normalized: self.analyticSolution = Lambda(lambda x: 1-x[0]**2) + elif self.parametrizedArg('polynomial').match(problem): + self.domainIndicator = domainIndicator + self.fluxIndicator = constant(0) + self.interactionIndicator = interactionIndicator+boundaryIndicator + self.fluxData = constant(0) + polyOrder = self.parametrizedArg('polynomial').interpret(problem)[0] + knownSolution = (((kType == FRACTIONAL and isinstance(sFun, (constFractionalOrder, variableConstFractionalOrder, singleVariableUnsymmetricFractionalOrder))) or + (kType in (INDICATOR, PERIDYNAMIC, GAUSSIAN))) and + phiFun is None and + normalized and + 0 <= polyOrder <= 3) + if polyOrder == 0: + self.rhsData = functionFactory('constant', 0.) + self.dirichletData = functionFactory('constant', 1.) + if knownSolution: + self.analyticSolution = self.dirichletData + elif polyOrder == 1: + self.rhsData = functionFactory('constant', 0.) + self.dirichletData = functionFactory('x0') + if knownSolution: + self.analyticSolution = self.dirichletData + elif polyOrder == 2: + self.rhsData = functionFactory('constant', -2) + self.dirichletData = functionFactory('x0**2') + if knownSolution: + self.analyticSolution = self.dirichletData + elif polyOrder == 3: + self.rhsData = -6*functionFactory('x0') + self.dirichletData = functionFactory('x0**3') + if knownSolution: + self.analyticSolution = self.dirichletData + else: + self.rhsData = functionFactory('Lambda', lambda x: -polyOrder*(polyOrder-1)*x[0]**(polyOrder-2)) + self.dirichletData = functionFactory('Lambda', lambda x: x[0]**polyOrder) + if knownSolution: + self.analyticSolution = self.dirichletData + + elif problem == 'exact-sin-Dirichlet': assert ((kType == INDICATOR) or (kType == FRACTIONAL)) and phiFun is None and normalized @@ -1142,19 +1259,23 @@ def getApproximationParams(self, dim, kernel, element, target_order): s = kernel.s if dim == 1: - self.eta = 1. if target_order <= 0.: if s is not None: target_order = (1+element-s.min)/dim else: target_order = 2. else: - self.eta = 3. if self.target_order <= 0.: target_order = 1/dim if element == 2: raise NotImplementedError() self.directlySetWithoutChecks('target_order', target_order) + else: + if target_order <= 0.: + target_order = 2/dim + self.directlySetWithoutChecks('target_order', target_order) + if dim == 1: + self.eta = 1. else: self.eta = 3. @@ -1187,8 +1308,8 @@ def setDriverArgs(self): 'rhs', 'rhsData', 'dirichletData', 'analyticSolution', 'exactL2Squared', 'exactHsSquared', 'initial']) - def processProblem(self, kernel, dim, domain, problem, normalized): - super().processProblem(kernel, dim, domain, problem, normalized) + def processProblem(self, kernel, dim, domain, domainParams, problem, normalized): + super().processProblem(kernel, dim, domain, domainParams, problem, normalized) steadyStateRHS = self.rhs steadyStateRHSdata = self.rhsData diff --git a/nl/setup.py b/nl/setup.py index 52d6308..4dc14f1 100644 --- a/nl/setup.py +++ b/nl/setup.py @@ -69,6 +69,8 @@ sources=[p.folder+"twoPointFunctions.pyx"]) p.addExtension("interactionDomains", sources=[p.folder+"interactionDomains.pyx"]) +p.addExtension("kernelNormalization", + sources=[p.folder+"kernelNormalization.pyx"]) p.addExtension("kernelsCy", sources=[p.folder+"kernelsCy.pyx"]) p.addExtension("fractionalOrders", diff --git a/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 b/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 index b60b494..def6aaf 100644 --- a/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.14090504680748486 - L2 error: 0.05317078748170508 - L2 error interpolated: 0.01773511238026114 - Linf error interpolated: 0.040238270520615516 - relative Hs error: 0.09580434396755147 - relative L2 error: 0.04268752259368201 - relative interpolated L2 error: 0.014231452347277154 - relative interpolated Linf error: 0.04675593019703401 + Hs error: 0.1409951440304447 + L2 error: 0.05317302332004089 + L2 error interpolated: 0.017742516045588845 + Linf error interpolated: 0.04024576450108175 + relative Hs error: 0.09586560298938593 + relative L2 error: 0.042689317609403914 + relative interpolated L2 error: 0.014237393381539993 + relative interpolated Linf error: 0.046764638027242504 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense b/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense index dd5f262..9a831c2 100644 --- a/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.14133245719249318 - L2 error: 0.053181189908821994 - L2 error interpolated: 0.017767471098211182 - Linf error interpolated: 0.04023373914321565 - relative Hs error: 0.09609494939630227 - relative L2 error: 0.04269587405627291 - relative interpolated L2 error: 0.014257418438872804 - relative interpolated Linf error: 0.04675066484236883 + Hs error: 0.1413321122331112 + L2 error: 0.05318118151618651 + L2 error interpolated: 0.01776734562313268 + Linf error interpolated: 0.04023229984790638 + relative Hs error: 0.09609471485106749 + relative L2 error: 0.04269586731834719 + relative interpolated L2 error: 0.014257317752022776 + relative interpolated Linf error: 0.0467489924148479 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 b/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 index 4a2e26b..de46d89 100644 --- a/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.06227042611223489 - L2 error: 0.03690130710763057 - L2 error interpolated: 0.055852834714174925 - Linf error interpolated: 0.2093154813930287 - relative Hs error: 0.042338989677306926 - relative L2 error: 0.029625767371517293 - relative interpolated L2 error: 0.04503979955178477 - relative interpolated Linf error: 0.2431971695425788 + Hs error: 0.3961994530890435 + L2 error: 0.05766669454676223 + L2 error interpolated: 0.017395282966069778 + Linf error interpolated: 0.08501026856170302 + relative Hs error: 0.26938445104353914 + relative L2 error: 0.046297007115323835 + relative interpolated L2 error: 0.014027579118370492 + relative interpolated Linf error: 0.09877079592331199 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense b/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense index 44e8e81..84f15c6 100644 --- a/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractional.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.06267848767877376 - L2 error: 0.03690199607863581 - L2 error interpolated: 0.05585332947223387 - Linf error interpolated: 0.20926548500139242 - relative Hs error: 0.04261643943206299 - relative L2 error: 0.02962632050354233 - relative interpolated L2 error: 0.045040198525335806 - relative interpolated Linf error: 0.24313908028490727 + Hs error: 0.39628788712384144 + L2 error: 0.05767235065160996 + L2 error interpolated: 0.01741402710931799 + Linf error interpolated: 0.08492200745915318 + relative Hs error: 0.2694445792283003 + relative L2 error: 0.046301548050580114 + relative interpolated L2 error: 0.014042694420198766 + relative interpolated Linf error: 0.09866824808414609 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 b/tests/cache_runFractional.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 index 9d8bcaa..7e30923 100644 --- a/tests/cache_runFractional.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 +++ b/tests/cache_runFractional.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.06227042587114324 - L2 error: 0.03690130661284065 - L2 error interpolated: 0.05585283301962591 - Linf error interpolated: 0.20931548802063804 - relative Hs error: 0.04233898951338357 - relative L2 error: 0.029625766974281195 - relative interpolated L2 error: 0.045039798185298206 - relative interpolated Linf error: 0.24319717724299267 + Hs error: 0.39619946901297326 + L2 error: 0.057666695197986645 + L2 error interpolated: 0.017395235311283137 + Linf error interpolated: 0.08500979051620516 + relative Hs error: 0.2693844618705583 + relative L2 error: 0.04629700763815147 + relative interpolated L2 error: 0.014027540689487715 + relative interpolated Linf error: 0.09877024049706633 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense b/tests/cache_runFractional.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense index 4d4270a..42865d4 100644 --- a/tests/cache_runFractional.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense +++ b/tests/cache_runFractional.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.06267848744064042 - L2 error: 0.0369019955875316 - L2 error interpolated: 0.05585332779981015 - Linf error interpolated: 0.20926549147124435 - relative Hs error: 0.042616439270151046 - relative L2 error: 0.029626320109265265 - relative interpolated L2 error: 0.04504019717669111 - relative interpolated Linf error: 0.24313908780202767 + Hs error: 0.3962879019940665 + L2 error: 0.05767235132032506 + L2 error interpolated: 0.01741398004575365 + Linf error interpolated: 0.0849215296918851 + relative Hs error: 0.2694445893388832 + relative L2 error: 0.046301548587449934 + relative interpolated L2 error: 0.014042656468078432 + relative interpolated Linf error: 0.09866769298116697 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 index 132d1f6..1f3abfb 100644 --- a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.08624586697655923 - L2 error: 0.013545369494356507 - L2 error interpolated: 0.009431667645849429 - Linf error interpolated: 0.03416022389360812 - relative Hs error: 0.06140946049895989 - relative L2 error: 0.00957802262313688 - relative interpolated L2 error: 0.006668482960460417 - relative interpolated Linf error: 0.03027417215179057 + Hs error: 0.08624507875463457 + L2 error: 0.013545339390748182 + L2 error interpolated: 0.009431622732729125 + Linf error interpolated: 0.034160202861473976 + relative Hs error: 0.061408899262986 + relative L2 error: 0.009578001336671294 + relative interpolated L2 error: 0.0066684512054846585 + relative interpolated Linf error: 0.03027415351226257 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense index e9b8575..4a1fc8f 100644 --- a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.0863471072165334 - L2 error: 0.01354844278444496 - L2 error interpolated: 0.009435986027159142 - Linf error interpolated: 0.034158032671413985 - relative Hs error: 0.061481546370846234 - relative L2 error: 0.009580195767398977 - relative interpolated L2 error: 0.006671536190627329 - relative interpolated Linf error: 0.03027223020205 + Hs error: 0.0863469994893122 + L2 error: 0.013548441805861365 + L2 error interpolated: 0.009435986433937102 + Linf error interpolated: 0.034158068810118125 + relative Hs error: 0.06148146966606296 + relative L2 error: 0.009580195075435881 + relative interpolated L2 error: 0.006671536478232023 + relative interpolated Linf error: 0.030272262229631374 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 index d6f8de4..a939848 100644 --- a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.03176949500712011 - L2 error: 0.005305491586028264 - L2 error interpolated: 0.0246963654473275 - Linf error interpolated: 0.23030513023108914 - relative Hs error: 0.02262076569120575 - relative L2 error: 0.003751549078008755 - relative interpolated L2 error: 0.017474927572792346 - relative interpolated Linf error: 0.2041026074806951 + Hs error: 0.1317754301519303 + L2 error: 0.026655322403489307 + L2 error interpolated: 0.008022646406895418 + Linf error interpolated: 0.046642138999338945 + relative Hs error: 0.09382777814556387 + relative L2 error: 0.018848159226220983 + relative interpolated L2 error: 0.005676752929560787 + relative interpolated Linf error: 0.04133551944192432 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense index de9586c..f36fae6 100644 --- a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.031769401751505486 - L2 error: 0.005305488971312755 - L2 error interpolated: 0.02469636381870763 - Linf error interpolated: 0.23030508944305492 - relative Hs error: 0.0226206992906096 - relative L2 error: 0.0037515472291256876 - relative interpolated L2 error: 0.017474926420395474 - relative interpolated Linf error: 0.20410257133324092 + Hs error: 0.13177539549804498 + L2 error: 0.026655318974538753 + L2 error interpolated: 0.008022651615535436 + Linf error interpolated: 0.04664221602283042 + relative Hs error: 0.0938277534710319 + relative L2 error: 0.018848156801586795 + relative interpolated L2 error: 0.005676756615147934 + relative interpolated Linf error: 0.041335587702216355 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatH2 b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatH2 index bff2f01..f84815d 100644 --- a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatH2 @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.02819524773557293 - L2 error: 0.005129502345279684 - L2 error interpolated: 0.020383923966215654 - Linf error interpolated: 0.21572832544966022 - relative Hs error: 0.020075802038683737 - relative L2 error: 0.0036271058924595627 - relative interpolated L2 error: 0.014416785593362078 - relative interpolated Linf error: 0.19118425059632427 + Hs error: 0.10841094415622979 + L2 error: 0.02292087168437646 + L2 error interpolated: 0.007554723669727594 + Linf error interpolated: 0.05090901291783434 + relative Hs error: 0.07719161307319655 + relative L2 error: 0.016207503798729313 + relative interpolated L2 error: 0.00534317296042096 + relative interpolated Linf error: 0.04511693799600705 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatdense b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatdense index 501186a..0036196 100644 --- a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatdense @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.028195243105532355 - L2 error: 0.005129494634925146 - L2 error interpolated: 0.020383922006257393 - Linf error interpolated: 0.21572822436629102 - relative Hs error: 0.020075798741965835 - relative L2 error: 0.0036271004404155836 - relative interpolated L2 error: 0.014416784207157043 - relative interpolated Linf error: 0.1911841610135208 + Hs error: 0.10841094254111872 + L2 error: 0.022920865169721243 + L2 error interpolated: 0.0075547464500849305 + Linf error interpolated: 0.05090918724915605 + relative Hs error: 0.07719161192319246 + relative L2 error: 0.01620749919217243 + relative interpolated L2 error: 0.005343189072114078 + relative interpolated Linf error: 0.0451170924931183 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatH2 b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatH2 index a50aad6..643be2e 100644 --- a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatH2 @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.020473310606527982 - L2 error: 0.0020890552696959737 - L2 error interpolated: 0.013333440519327819 - Linf error interpolated: 0.18369560001689372 - relative Hs error: 0.01457756763366094 - relative L2 error: 0.0014771851474755144 - relative interpolated L2 error: 0.009430285168640204 - relative interpolated Linf error: 0.16279598682217128 + Hs error: 0.1101929993006766 + L2 error: 0.014277238642413542 + L2 error interpolated: 0.005063067540537936 + Linf error interpolated: 0.022083796893636487 + relative Hs error: 0.0784604859924011 + relative L2 error: 0.01009553226066923 + relative interpolated L2 error: 0.0035809340182038456 + relative interpolated Linf error: 0.019571255423370626 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatdense b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatdense index e4076cf..b35208f 100644 --- a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatdense @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.020472519486592667 - L2 error: 0.0020888387979283995 - L2 error interpolated: 0.013333449492928423 - Linf error interpolated: 0.18369276304121315 - relative Hs error: 0.014577004334222708 - relative L2 error: 0.0014770320788207272 - relative interpolated L2 error: 0.009430291515360135 - relative interpolated Linf error: 0.16279347261793628 + Hs error: 0.11019060819971181 + L2 error: 0.014276977312524999 + L2 error interpolated: 0.005064851277694193 + Linf error interpolated: 0.022093053095892323 + relative Hs error: 0.07845878346188696 + relative L2 error: 0.010095347472532913 + relative interpolated L2 error: 0.003582195594315804 + relative interpolated Linf error: 0.019579458519037182 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemknownSolution--elementP1--solvercg-jacobi--matrixFormatH2 b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemknownSolution--elementP1--solvercg-jacobi--matrixFormatH2 index dcf55eb..b6dd9e8 100644 --- a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemknownSolution--elementP1--solvercg-jacobi--matrixFormatH2 +++ b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemknownSolution--elementP1--solvercg-jacobi--matrixFormatH2 @@ -1,11 +1,11 @@ Timers: {} errors: - L2 error: 0.000607392058826574 - L2 error interpolated: 0.0014406299046700504 - Linf error interpolated: 0.01257048139671692 - relative L2 error: 0.000553353562595389 - relative interpolated L2 error: 0.001312548064580926 - relative interpolated Linf error: 0.01257048139671692 + L2 error: 0.001493004043140445 + L2 error interpolated: 0.000526382451040797 + Linf error interpolated: 0.0036717119596183373 + relative L2 error: 0.0013601743622350762 + relative interpolated L2 error: 0.0004795834551978154 + relative interpolated Linf error: 0.0036717119596183373 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemzeroFlux--elementP1--solverlu--matrixFormatH2 b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemzeroFlux--elementP1--solverlu--matrixFormatH2 index 1c725c6..71e3bc4 100644 --- a/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemzeroFlux--elementP1--solverlu--matrixFormatH2 +++ b/tests/cache_runFractional.py--domaininterval--sconst(0.25)--problemzeroFlux--elementP1--solverlu--matrixFormatH2 @@ -1,11 +1,11 @@ Timers: {} errors: - L2 error: 4.1057062382325376e-05 - L2 error interpolated: 4.105704754817684e-05 - Linf error interpolated: 0.0003000263738740516 - relative L2 error: 5.028442658715754e-05 - relative interpolated L2 error: 5.028440841911015e-05 - relative interpolated Linf error: 0.0003000263738740516 + L2 error: 4.105702722904196e-05 + L2 error interpolated: 4.105701007094699e-05 + Linf error interpolated: 0.00030002633007253365 + relative L2 error: 5.028438353335397e-05 + relative interpolated L2 error: 5.0284362519065095e-05 + relative interpolated Linf error: 0.00030002633007253365 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 b/tests/cache_runFractional.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 index 4891fe0..7826a08 100644 --- a/tests/cache_runFractional.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 +++ b/tests/cache_runFractional.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.031769495153048444 - L2 error: 0.005305492505303933 - L2 error interpolated: 0.024696365064832512 - Linf error interpolated: 0.23030513603944397 - relative Hs error: 0.022620765795110787 - relative L2 error: 0.003751549728034814 - relative interpolated L2 error: 0.01747492730214231 - relative interpolated Linf error: 0.20410261262821552 + Hs error: 0.13177542732204822 + L2 error: 0.02665532198266343 + L2 error interpolated: 0.008022634968889985 + Linf error interpolated: 0.04664209126519353 + relative Hs error: 0.09382777613060897 + relative L2 error: 0.018848158928652152 + relative interpolated L2 error: 0.00567674483613032 + relative interpolated Linf error: 0.04133547713863939 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractional.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense b/tests/cache_runFractional.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense index ec134a6..b790c4d 100644 --- a/tests/cache_runFractional.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense +++ b/tests/cache_runFractional.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense @@ -1,13 +1,13 @@ Timers: {} errors: - Hs error: 0.031769401897318916 - L2 error: 0.005305489889940173 - L2 error interpolated: 0.024696363435431367 - Linf error interpolated: 0.23030509524488546 - relative Hs error: 0.022620699394432828 - relative L2 error: 0.0037515478786933646 - relative interpolated L2 error: 0.017474926149192615 - relative interpolated Linf error: 0.20410257647497937 + Hs error: 0.13177539266817143 + L2 error: 0.026655318553675336 + L2 error interpolated: 0.008022640177556224 + Linf error interpolated: 0.04664216828925877 + relative Hs error: 0.09382775145608307 + relative L2 error: 0.01884815650399142 + relative interpolated L2 error: 0.005676748521736022 + relative interpolated Linf error: 0.041335545399439906 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 b/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 index 76881b5..e8e9dc0 100644 --- a/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.06330842381027034 - L^2(0,T; L^2(Omega)) norm: 1.4944254481530928 - L^2(Omega) error at t=finalTime: 0.02888395949052554 + L^2(0,T; L^2(Omega)) error: 0.06330900779185128 + L^2(0,T; L^2(Omega)) norm: 1.49441738331053 + L^2(Omega) error at t=finalTime: 0.028884938922221277 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense b/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense index 5703507..8037e67 100644 --- a/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.06331101647548046 - L^2(0,T; L^2(Omega)) norm: 1.494387020352447 - L^2(Omega) error at t=finalTime: 0.028888346628916803 + L^2(0,T; L^2(Omega)) error: 0.06331101527900383 + L^2(0,T; L^2(Omega)) norm: 1.49438704354491 + L^2(Omega) error at t=finalTime: 0.02888834024980095 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 b/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 index 62a1368..ded6677 100644 --- a/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.05071531506350525 - L^2(0,T; L^2(Omega)) norm: 1.4949918500501762 - L^2(Omega) error at t=finalTime: 0.020060587191451665 + L^2(0,T; L^2(Omega)) error: 0.07032665300632383 + L^2(0,T; L^2(Omega)) norm: 1.4897692989784312 + L^2(Omega) error at t=finalTime: 0.03169469529581403 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense b/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense index ef40374..93f4136 100644 --- a/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractionalHeat.py--domaindisc--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.050715836753152437 - L^2(0,T; L^2(Omega)) norm: 1.4949755875568616 - L^2(Omega) error at t=finalTime: 0.020060631770579762 + L^2(0,T; L^2(Omega)) error: 0.07032815147555171 + L^2(0,T; L^2(Omega)) norm: 1.4897452496349892 + L^2(Omega) error at t=finalTime: 0.03169764874285769 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 b/tests/cache_runFractionalHeat.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 index 6e58187..65947f1 100644 --- a/tests/cache_runFractionalHeat.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 +++ b/tests/cache_runFractionalHeat.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.05071531544578167 - L^2(0,T; L^2(Omega)) norm: 1.4949918485967646 - L^2(Omega) error at t=finalTime: 0.0200605891890031 + L^2(0,T; L^2(Omega)) error: 0.07032664347646923 + L^2(0,T; L^2(Omega)) norm: 1.4897692989359879 + L^2(Omega) error at t=finalTime: 0.031694695847278256 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense b/tests/cache_runFractionalHeat.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense index 8d36f8e..fad1824 100644 --- a/tests/cache_runFractionalHeat.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense +++ b/tests/cache_runFractionalHeat.py--domaindisc--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.050715837134377995 - L^2(0,T; L^2(Omega)) norm: 1.4949755861023193 - L^2(Omega) error at t=finalTime: 0.020060633763670236 + L^2(0,T; L^2(Omega)) error: 0.07032814196103428 + L^2(0,T; L^2(Omega)) norm: 1.4897452495944319 + L^2(Omega) error at t=finalTime: 0.03169764930080015 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 index 7328612..1f1dad2 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatH2 @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.01494132768387929 - L^2(0,T; L^2(Omega)) norm: 1.7025655539703572 - L^2(Omega) error at t=finalTime: 0.007567188266998093 + L^2(0,T; L^2(Omega)) error: 0.01494132128938275 + L^2(0,T; L^2(Omega)) norm: 1.7025655960568389 + L^2(Omega) error at t=finalTime: 0.0075671787595630235 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense index 8c901b0..45c7090 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP0--solvercg-mg--matrixFormatdense @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.01494131148662652 - L^2(0,T; L^2(Omega)) norm: 1.7025600847367917 - L^2(Omega) error at t=finalTime: 0.0075677600770572865 + L^2(0,T; L^2(Omega)) error: 0.0149413089985309 + L^2(0,T; L^2(Omega)) norm: 1.7025600858867103 + L^2(Omega) error at t=finalTime: 0.007567757829891671 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 index ed58ba5..e397a7b 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatH2 @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.011178762917090485 - L^2(0,T; L^2(Omega)) norm: 1.703480580844389 - L^2(Omega) error at t=finalTime: 0.0026169275065709096 + L^2(0,T; L^2(Omega)) error: 0.03218338879054702 + L^2(0,T; L^2(Omega)) norm: 1.7018299475043261 + L^2(Omega) error at t=finalTime: 0.014558725749186825 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense index f0c135f..448fe76 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP1--solvercg-mg--matrixFormatdense @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.01117876176160322 - L^2(0,T; L^2(Omega)) norm: 1.703480582652017 - L^2(Omega) error at t=finalTime: 0.0026169255320408783 + L^2(0,T; L^2(Omega)) error: 0.03218338586612875 + L^2(0,T; L^2(Omega)) norm: 1.7018299503210628 + L^2(Omega) error at t=finalTime: 0.01455872345929613 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatH2 b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatH2 index 88602a5..0b92f53 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatH2 @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.0114062977758145 - L^2(0,T; L^2(Omega)) norm: 1.7032132520631214 - L^2(Omega) error at t=finalTime: 0.0024190469829229516 + L^2(0,T; L^2(Omega)) error: 0.02749586871028784 + L^2(0,T; L^2(Omega)) norm: 1.701925956510871 + L^2(Omega) error at t=finalTime: 0.012420538758452487 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatdense b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatdense index baacacd..838dbb2 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP2--solvercg-mg--matrixFormatdense @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.011406294738598916 - L^2(0,T; L^2(Omega)) norm: 1.7032132520521739 - L^2(Omega) error at t=finalTime: 0.0024190416441285012 + L^2(0,T; L^2(Omega)) error: 0.027495862469873365 + L^2(0,T; L^2(Omega)) norm: 1.7019259587916384 + L^2(Omega) error at t=finalTime: 0.012420534279834644 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatH2 b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatH2 index eec77bb..297b1c1 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatH2 +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatH2 @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.006802049474540982 - L^2(0,T; L^2(Omega)) norm: 1.703329914967122 - L^2(Omega) error at t=finalTime: 0.00031938257264443675 + L^2(0,T; L^2(Omega)) error: 0.01726746578963863 + L^2(0,T; L^2(Omega)) norm: 1.702633036816579 + L^2(Omega) error at t=finalTime: 0.007746454624948888 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatdense b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatdense index 643f8c5..81341d5 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatdense +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemconstant--elementP3--solvercg-mg--matrixFormatdense @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.006802037231113124 - L^2(0,T; L^2(Omega)) norm: 1.7033299227908947 - L^2(Omega) error at t=finalTime: 0.0003189154238379492 + L^2(0,T; L^2(Omega)) error: 0.017267223086710897 + L^2(0,T; L^2(Omega)) norm: 1.7026331344615124 + L^2(Omega) error at t=finalTime: 0.007746289486904896 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemknownSolution--elementP1--solvercg-jacobi--matrixFormatH2 b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemknownSolution--elementP1--solvercg-jacobi--matrixFormatH2 index 9808fd8..25c054d 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemknownSolution--elementP1--solvercg-jacobi--matrixFormatH2 +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemknownSolution--elementP1--solvercg-jacobi--matrixFormatH2 @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.0009868341060745847 - L^2(0,T; L^2(Omega)) norm: 1.3229163786976783 - L^2(Omega) error at t=finalTime: 0.0005483083851510518 + L^2(0,T; L^2(Omega)) error: 0.0018388585398440504 + L^2(0,T; L^2(Omega)) norm: 1.3228634831094461 + L^2(Omega) error at t=finalTime: 0.0008912208343986159 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemzeroFlux--elementP1--solverlu--matrixFormatH2 b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemzeroFlux--elementP1--solverlu--matrixFormatH2 index a78c7c3..fba0cf9 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemzeroFlux--elementP1--solverlu--matrixFormatH2 +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconst(0.25)--problemzeroFlux--elementP1--solverlu--matrixFormatH2 @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.00027586671804719504 - L^2(0,T; L^2(Omega)) norm: 0.9841382432420539 - L^2(Omega) error at t=finalTime: 0.00041205293671547606 + L^2(0,T; L^2(Omega)) error: 0.00027586673462306533 + L^2(0,T; L^2(Omega)) norm: 0.9841382432557984 + L^2(Omega) error at t=finalTime: 0.0004120529549024724 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 b/tests/cache_runFractionalHeat.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 index c7543aa..62bc5a4 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatH2 @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.011178755913442575 - L^2(0,T; L^2(Omega)) norm: 1.7034805781408084 - L^2(Omega) error at t=finalTime: 0.0026169321223716243 + L^2(0,T; L^2(Omega)) error: 0.03218339459660203 + L^2(0,T; L^2(Omega)) norm: 1.7018299504634995 + L^2(Omega) error at t=finalTime: 0.014558732598655947 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runFractionalHeat.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense b/tests/cache_runFractionalHeat.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense index 78701b5..8938cfc 100644 --- a/tests/cache_runFractionalHeat.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense +++ b/tests/cache_runFractionalHeat.py--domaininterval--sconstantNonSym(0.25)--problemconstant--elementP1--solvergmres-jacobi--matrixFormatdense @@ -1,8 +1,8 @@ Timers: {} errors: - L^2(0,T; L^2(Omega)) error: 0.011178754757950861 - L^2(0,T; L^2(Omega)) norm: 1.7034805799484114 - L^2(Omega) error at t=finalTime: 0.0026169301478875004 + L^2(0,T; L^2(Omega)) error: 0.032183391672112704 + L^2(0,T; L^2(Omega)) norm: 1.7018299532802796 + L^2(Omega) error at t=finalTime: 0.014558730308751077 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runNonlocal.py--domaininterval--kernelTypeconstant--problempoly-Dirichlet--solverlu--matrixFormatH2 b/tests/cache_runNonlocal.py--domaininterval--kernelTypeconstant--problempoly-Dirichlet--solverlu--matrixFormatH2 new file mode 100644 index 0000000..78d8374 --- /dev/null +++ b/tests/cache_runNonlocal.py--domaininterval--kernelTypeconstant--problempoly-Dirichlet--solverlu--matrixFormatH2 @@ -0,0 +1,9 @@ +Timers: {} +errors: + L2 error interpolated: 2.2380604839728094e-13 + Linf error interpolated: 4.313216450668733e-13 + relative interpolated L2 error: 2.142368730105727e-13 + relative interpolated Linf error: 4.313216450668733e-13 +meshes: {} +results: {} +vectors: {} diff --git a/tests/cache_runNonlocal.py--domaininterval--kernelTypeconstant--problempoly-Neumann--solverlu--matrixFormatH2 b/tests/cache_runNonlocal.py--domaininterval--kernelTypeconstant--problempoly-Neumann--solverlu--matrixFormatH2 new file mode 100644 index 0000000..d5c9323 --- /dev/null +++ b/tests/cache_runNonlocal.py--domaininterval--kernelTypeconstant--problempoly-Neumann--solverlu--matrixFormatH2 @@ -0,0 +1,9 @@ +Timers: {} +errors: + L2 error interpolated: 0.4718187556236188 + Linf error interpolated: 0.7335954576992952 + relative interpolated L2 error: 0.4516454115355879 + relative interpolated Linf error: 0.7335954576992952 +meshes: {} +results: {} +vectors: {} diff --git a/tests/cache_runNonlocal.py--domaininterval--kernelTypefractional--problempoly-Dirichlet--solverlu--matrixFormatH2 b/tests/cache_runNonlocal.py--domaininterval--kernelTypefractional--problempoly-Dirichlet--solverlu--matrixFormatH2 new file mode 100644 index 0000000..70f4bce --- /dev/null +++ b/tests/cache_runNonlocal.py--domaininterval--kernelTypefractional--problempoly-Dirichlet--solverlu--matrixFormatH2 @@ -0,0 +1,9 @@ +Timers: {} +errors: + L2 error interpolated: 3.74958407583645e-11 + Linf error interpolated: 3.69601016458887e-11 + relative interpolated L2 error: 3.5892647819396405e-11 + relative interpolated Linf error: 3.69601016458887e-11 +meshes: {} +results: {} +vectors: {} diff --git a/tests/cache_runNonlocal.py--domaininterval--kernelTypefractional--problempoly-Neumann--solverlu--matrixFormatH2 b/tests/cache_runNonlocal.py--domaininterval--kernelTypefractional--problempoly-Neumann--solverlu--matrixFormatH2 new file mode 100644 index 0000000..64e5421 --- /dev/null +++ b/tests/cache_runNonlocal.py--domaininterval--kernelTypefractional--problempoly-Neumann--solverlu--matrixFormatH2 @@ -0,0 +1,9 @@ +Timers: {} +errors: + L2 error interpolated: 0.5476731873817714 + Linf error interpolated: 0.9515340667676401 + relative interpolated L2 error: 0.5242565691885458 + relative interpolated Linf error: 0.9515340667676401 +meshes: {} +results: {} +vectors: {} diff --git a/tests/cache_runNonlocal.py--domaininterval--kernelTypeinverseDistance--problempoly-Dirichlet--solverlu--matrixFormatH2 b/tests/cache_runNonlocal.py--domaininterval--kernelTypeinverseDistance--problempoly-Dirichlet--solverlu--matrixFormatH2 new file mode 100644 index 0000000..597ca61 --- /dev/null +++ b/tests/cache_runNonlocal.py--domaininterval--kernelTypeinverseDistance--problempoly-Dirichlet--solverlu--matrixFormatH2 @@ -0,0 +1,9 @@ +Timers: {} +errors: + L2 error interpolated: 2.0677307416367313e-07 + Linf error interpolated: 2.2039557501241092e-07 + relative interpolated L2 error: 1.9793217005902312e-07 + relative interpolated Linf error: 2.2039557501241092e-07 +meshes: {} +results: {} +vectors: {} diff --git a/tests/cache_runNonlocal.py--domaininterval--kernelTypeinverseDistance--problempoly-Neumann--solverlu--matrixFormatH2 b/tests/cache_runNonlocal.py--domaininterval--kernelTypeinverseDistance--problempoly-Neumann--solverlu--matrixFormatH2 new file mode 100644 index 0000000..854bb88 --- /dev/null +++ b/tests/cache_runNonlocal.py--domaininterval--kernelTypeinverseDistance--problempoly-Neumann--solverlu--matrixFormatH2 @@ -0,0 +1,9 @@ +Timers: {} +errors: + L2 error interpolated: 0.5034968338389819 + Linf error interpolated: 0.8477624694641785 + relative interpolated L2 error: 0.48196904429012694 + relative interpolated Linf error: 0.8477624694641785 +meshes: {} +results: {} +vectors: {} diff --git a/tests/cache_runNonlocal.py--domainsquare--kernelTypeconstant--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 b/tests/cache_runNonlocal.py--domainsquare--kernelTypeconstant--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 new file mode 100644 index 0000000..48a6abc --- /dev/null +++ b/tests/cache_runNonlocal.py--domainsquare--kernelTypeconstant--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 @@ -0,0 +1,9 @@ +Timers: {} +errors: + L2 error interpolated: 0.011552136188779451 + Linf error interpolated: 0.009927687727526946 + relative interpolated L2 error: 0.00714925855055064 + relative interpolated Linf error: 0.009927687727526946 +meshes: {} +results: {} +vectors: {} diff --git a/tests/cache_runNonlocal.py--domainsquare--kernelTypeconstant--problempoly-Dirichlet--solvercg-mg--matrixFormatdense b/tests/cache_runNonlocal.py--domainsquare--kernelTypeconstant--problempoly-Dirichlet--solvercg-mg--matrixFormatdense index 033c18f..41902cf 100644 --- a/tests/cache_runNonlocal.py--domainsquare--kernelTypeconstant--problempoly-Dirichlet--solvercg-mg--matrixFormatdense +++ b/tests/cache_runNonlocal.py--domainsquare--kernelTypeconstant--problempoly-Dirichlet--solvercg-mg--matrixFormatdense @@ -1,9 +1,9 @@ Timers: {} errors: - L2 error interpolated: 0.010620115381058651 - Linf error interpolated: 0.008877474825595977 - relative interpolated L2 error: 0.006572459799220038 - relative interpolated Linf error: 0.008877474825595977 + L2 error interpolated: 0.011660745699545316 + Linf error interpolated: 0.009993667722354216 + relative interpolated L2 error: 0.00721647360591574 + relative interpolated Linf error: 0.009993667722354216 meshes: {} results: {} vectors: {} diff --git a/tests/cache_runNonlocal.py--domainsquare--kernelTypefractional--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 b/tests/cache_runNonlocal.py--domainsquare--kernelTypefractional--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 new file mode 100644 index 0000000..80653b5 --- /dev/null +++ b/tests/cache_runNonlocal.py--domainsquare--kernelTypefractional--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 @@ -0,0 +1,9 @@ +Timers: {} +errors: + L2 error interpolated: 0.0033627264610840775 + Linf error interpolated: 0.002892688108887098 + relative interpolated L2 error: 0.0020810870398514842 + relative interpolated Linf error: 0.002892688108887098 +meshes: {} +results: {} +vectors: {} diff --git a/tests/cache_runNonlocal.py--domainsquare--kernelTypeinverseDistance--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 b/tests/cache_runNonlocal.py--domainsquare--kernelTypeinverseDistance--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 new file mode 100644 index 0000000..854eb60 --- /dev/null +++ b/tests/cache_runNonlocal.py--domainsquare--kernelTypeinverseDistance--problempoly-Dirichlet--solvercg-mg--matrixFormatH2 @@ -0,0 +1,9 @@ +Timers: {} +errors: + L2 error interpolated: 0.008596666193741862 + Linf error interpolated: 0.007407120991754779 + relative interpolated L2 error: 0.005320209897761966 + relative interpolated Linf error: 0.007407120991754779 +meshes: {} +results: {} +vectors: {} diff --git a/tests/cache_runNonlocal.py--domainsquare--kernelTypeinverseDistance--problempoly-Dirichlet--solvercg-mg--matrixFormatdense b/tests/cache_runNonlocal.py--domainsquare--kernelTypeinverseDistance--problempoly-Dirichlet--solvercg-mg--matrixFormatdense index 4c39a23..4ca2967 100644 --- a/tests/cache_runNonlocal.py--domainsquare--kernelTypeinverseDistance--problempoly-Dirichlet--solvercg-mg--matrixFormatdense +++ b/tests/cache_runNonlocal.py--domainsquare--kernelTypeinverseDistance--problempoly-Dirichlet--solvercg-mg--matrixFormatdense @@ -1,9 +1,9 @@ Timers: {} errors: - L2 error interpolated: 0.0078163078883207 - Linf error interpolated: 0.006546634974906151 - relative interpolated L2 error: 0.004837270362046979 - relative interpolated Linf error: 0.006546634974906151 + L2 error interpolated: 0.00867577950640197 + Linf error interpolated: 0.007455365403179481 + relative interpolated L2 error: 0.0053691706715751275 + relative interpolated Linf error: 0.007455365403179481 meshes: {} results: {} vectors: {} diff --git a/tests/test_drivers_intFracLapl.py b/tests/test_drivers_intFracLapl.py index cea9263..6f48360 100644 --- a/tests/test_drivers_intFracLapl.py +++ b/tests/test_drivers_intFracLapl.py @@ -21,18 +21,30 @@ def idfunc(param): @pytest.fixture(scope='module', params=[ - ('interval', 'fractional', 'poly-Dirichlet', 'lu'), - ('interval', 'fractional', 'poly-Neumann', 'lu'), - ('interval', 'constant', 'poly-Dirichlet', 'lu'), - ('interval', 'constant', 'poly-Neumann', 'lu'), - ('interval', 'inverseDistance', 'poly-Dirichlet', 'lu'), - ('interval', 'inverseDistance', 'poly-Neumann', 'lu'), - ('square', 'fractional', 'poly-Dirichlet', 'cg-mg'), - ('square', 'fractional', 'poly-Neumann', 'cg-mg'), - ('square', 'constant', 'poly-Dirichlet', 'cg-mg'), - ('square', 'constant', 'poly-Neumann', 'cg-mg'), - ('square', 'inverseDistance', 'poly-Dirichlet', 'cg-mg'), - ('square', 'inverseDistance', 'poly-Neumann', 'cg-mg'), + ('interval', 'fractional', 'poly-Dirichlet', 'lu', 'dense'), + ('interval', 'fractional', 'poly-Neumann', 'lu', 'dense'), + ('interval', 'constant', 'poly-Dirichlet', 'lu', 'dense'), + ('interval', 'constant', 'poly-Neumann', 'lu', 'dense'), + ('interval', 'inverseDistance', 'poly-Dirichlet', 'lu', 'dense'), + ('interval', 'inverseDistance', 'poly-Neumann', 'lu', 'dense'), + ('square', 'fractional', 'poly-Dirichlet', 'cg-mg', 'dense'), + ('square', 'fractional', 'poly-Neumann', 'cg-mg', 'dense'), + ('square', 'constant', 'poly-Dirichlet', 'cg-mg', 'dense'), + ('square', 'constant', 'poly-Neumann', 'cg-mg', 'dense'), + ('square', 'inverseDistance', 'poly-Dirichlet', 'cg-mg', 'dense'), + ('square', 'inverseDistance', 'poly-Neumann', 'cg-mg', 'dense'), + ('interval', 'fractional', 'poly-Dirichlet', 'lu', 'H2'), + ('interval', 'fractional', 'poly-Neumann', 'lu', 'H2'), + ('interval', 'constant', 'poly-Dirichlet', 'lu', 'H2'), + ('interval', 'constant', 'poly-Neumann', 'lu', 'H2'), + ('interval', 'inverseDistance', 'poly-Dirichlet', 'lu', 'H2'), + ('interval', 'inverseDistance', 'poly-Neumann', 'lu', 'H2'), + ('square', 'fractional', 'poly-Dirichlet', 'cg-mg', 'H2'), + ('square', 'fractional', 'poly-Neumann', 'cg-mg', 'H2'), + ('square', 'constant', 'poly-Dirichlet', 'cg-mg', 'H2'), + ('square', 'constant', 'poly-Neumann', 'cg-mg', 'H2'), + ('square', 'inverseDistance', 'poly-Dirichlet', 'cg-mg', 'H2'), + ('square', 'inverseDistance', 'poly-Neumann', 'cg-mg', 'H2'), ], ids=idfunc) def runNonlocal_params(request): @@ -41,15 +53,15 @@ def runNonlocal_params(request): @pytest.mark.slow def testNonlocal(runNonlocal_params, extra): - domain, kernel, problem, solver = runNonlocal_params + domain, kernel, problem, solver, matrixFormat = runNonlocal_params base = getPath()+'/../' py = ['runNonlocal.py', '--domain', domain, '--kernelType', kernel, '--problem', problem, - '--solver', solver] + '--solver', solver, + '--matrixFormat', matrixFormat] # if kernel != 'fractional': - py += ['--matrixFormat', 'dense'] path = base+'drivers' cacheDir = getPath()+'/' if problem == 'poly-Neumann' and domain == 'square': diff --git a/tests/test_fracLapl.py b/tests/test_fracLapl.py index a7e60df..acefb68 100644 --- a/tests/test_fracLapl.py +++ b/tests/test_fracLapl.py @@ -17,8 +17,8 @@ from PyNucleus.base.myTypes import REAL from scipy.special import gamma from PyNucleus.nl.kernels import getFractionalKernel -from PyNucleus.nl.fractionalOrders import (constFractionalOrder, - variableFractionalLaplacianScaling) +from PyNucleus.nl.fractionalOrders import constFractionalOrder +from PyNucleus.nl.kernelNormalization import variableFractionalLaplacianScaling import pytest import logging LOGGER = logging.getLogger(__name__) diff --git a/tests/test_h2finiteHorizon.py b/tests/test_h2finiteHorizon.py index 44c4693..9581895 100644 --- a/tests/test_h2finiteHorizon.py +++ b/tests/test_h2finiteHorizon.py @@ -13,6 +13,7 @@ leftRightFractionalOrder, variableConstFractionalOrder) from PyNucleus.nl.nonlocalLaplacian import nonlocalBuilder +from PyNucleus.nl.kernelNormalization import variableFractionalLaplacianScaling from PyNucleus.nl.kernels import getFractionalKernel from scipy.linalg import solve import pytest @@ -108,17 +109,23 @@ def test_h2_finite(kernels): # plt.colorbar() # print(err2.max()) - print('\nCORRECTED\n') - builder2 = nonlocalBuilder(mesh2, dm2, kernel1, zeroExterior=False) - A2 = builder2.getH2FiniteHorizon() - A2.setKernel(kernel1) + + buildCorrected = not isinstance(kernel1.scaling, variableFractionalLaplacianScaling) + + if buildCorrected: + print('\nCORRECTED\n') + builder2 = nonlocalBuilder(mesh2, dm2, kernel1, zeroExterior=False) + A2 = builder2.getH2FiniteHorizon() + A2.setKernel(kernel1) A1d = A1.toarray()[np.ix_(idx, idx)] A1_h2d = A1_h2.toarray()[np.ix_(idx, idx)] A1_h2_neard = A1_h2.Anear.toarray()[np.ix_(idx, idx)] A1_neard = A1d.copy() A1_neard[np.where(np.absolute(A1_h2_neard) == 0.)] = 0. - A2d = A2.toarray() + + if buildCorrected: + A2d = A2.toarray() nn = np.absolute(A1d) nn[nn < 1e-16] = 1. @@ -131,13 +138,14 @@ def test_h2_finite(kernels): errDenseH2_near_rel = errDenseH2_near/nn print('errDenseH2_near', errDenseH2_near.max(), errDenseH2_near_rel.max()) - errDenseCor = np.absolute(A1d-A2d) - errDenseCor_rel = errDenseCor/nn - print('errDenseCor', errDenseCor.max(), errDenseCor_rel.max()) + if buildCorrected: + errDenseCor = np.absolute(A1d-A2d) + errDenseCor_rel = errDenseCor/nn + print('errDenseCor', errDenseCor.max(), errDenseCor_rel.max()) - errH2Cor = np.absolute(A1_h2d-A2d) - errH2Cor_rel = errDenseCor/nn - print('errH2Cor', errH2Cor.max(), errH2Cor_rel.max()) + errH2Cor = np.absolute(A1_h2d-A2d) + errH2Cor_rel = errDenseCor/nn + print('errH2Cor', errH2Cor.max(), errH2Cor_rel.max()) # c = dm1.getDoFCoordinates()[idx, 0] # X, Y = np.meshgrid(c, c) @@ -184,24 +192,27 @@ def test_h2_finite(kernels): x1_h2 = np.zeros((A1.shape[0])) x1_h2[idx] = y1_h2 - b2 = dm2.assembleRHS(rhs).toarray() - y2 = solve(A2d, b2) - x2 = np.zeros((A1.shape[0])) - x2[idx] = y2 + if buildCorrected: + b2 = dm2.assembleRHS(rhs).toarray() + y2 = solve(A2d, b2) + x2 = np.zeros((A1.shape[0])) + x2[idx] = y2 # assert np.absolute(A1d[np.ix_(idx, idx)]-A2d).max() < 1e-5 M = dm1.assembleMass() L2_denseH2 = np.sqrt(abs(np.vdot(M*(x1-x1_h2), x1-x1_h2))) - L2_denseCor = np.sqrt(abs(np.vdot(M*(x1-x2), x1-x2))) - L2_H2Cor = np.sqrt(abs(np.vdot(M*(x1_h2-x2), x1_h2-x2))) + if buildCorrected: + L2_denseCor = np.sqrt(abs(np.vdot(M*(x1-x2), x1-x2))) + L2_H2Cor = np.sqrt(abs(np.vdot(M*(x1_h2-x2), x1_h2-x2))) L2_dense = np.sqrt(abs(np.vdot(M*x1, x1))) # L2_cor = np.sqrt(abs(np.vdot(M*x2, x2))) # if not (L2/L2_1 < mesh2.h**(0.5+min(kernel1.s.min, 0.5))): print('L2 errDenseH2', L2_denseH2) - print('L2 errDenseCor', L2_denseCor) - print('L2 errH2Cor', L2_H2Cor) + if buildCorrected: + print('L2 errDenseCor', L2_denseCor) + print('L2 errH2Cor', L2_H2Cor) # mesh1.plotFunction(x1, DoFMap=dm1, label='dense') # mesh1.plotFunction(x1_h2, DoFMap=dm1, label='h2') @@ -209,7 +220,8 @@ def test_h2_finite(kernels): # plt.legend() # plt.show() - assert L2_denseCor/L2_dense < mesh2.h**(0.5+min(kernel1.s.min, 0.5)), (L2_denseCor, L2_dense, L2_denseCor/L2_dense, mesh2.h**(0.5+min(kernel1.s.min, 0.5))) + if buildCorrected: + assert L2_denseCor/L2_dense < mesh2.h**(0.5+min(kernel1.s.min, 0.5)), (L2_denseCor, L2_dense, L2_denseCor/L2_dense, mesh2.h**(0.5+min(kernel1.s.min, 0.5))) mesh3 = meshOverlap(dim, kernel2.horizon.value) dm3 = P1_DoFMap(mesh3) @@ -217,64 +229,65 @@ def test_h2_finite(kernels): ind = dm3.interpolate(Lambda(lambda x: abs(x[0]) < 1-1e-12)) idx = ind.toarray() > 0 - print('\nCORRECTED\n') - A2.setKernel(kernel2) - - print('\nDENSE\n') - builder3 = nonlocalBuilder(mesh3, dm3, kernel2, zeroExterior=False) - A3 = builder3.getDense() - - A2d = A2.toarray() - A3d = A3.toarray()[np.ix_(idx, idx)] - - nn = np.absolute(A3d) - nn[nn < 1e-16] = 1. - - errDenseCor = np.absolute(A3d-A2d) - errDenseCor_rel = errDenseCor/nn - print('errDenseCor', errDenseCor.max(), errDenseCor_rel.max()) - - y2 = solve(A2d, b2) - x2 = np.zeros((A3.shape[0])) - x2[idx] = y2 - - b3 = dm3.assembleRHS(rhs).toarray() - y3 = solve(A3d, b3[idx]) - x3 = np.zeros((A3.shape[0])) - x3[idx] = y3 - - # assert np.absolute(A3d[np.ix_(idx, idx)]-A2d).max() < 1e-5 - - M = dm3.assembleMass() - - L2_denseCor = np.sqrt(abs(np.vdot(M*(x3-x2), x3-x2))) - L2_dense = np.sqrt(abs(np.vdot(M*x3, x3))) - # L2_cor = np.sqrt(abs(np.vdot(M*x2, x2))) - - print('L2 errDenseCor', L2_denseCor) - - # mesh3.plotFunction(x2, DoFMap=dm3, label='corrected') - # mesh3.plotFunction(x3, DoFMap=dm3, label='dense') - # plt.legend() - # plt.show() - - # if not (L2 < mesh2.h**(0.5+min(kernel2.s.value, 0.5))): - # mesh3.plotFunction(x3, DoFMap=dm3) - # mesh3.plotFunction(x2, DoFMap=dm3) - # plt.figure() - # for lvl in A2.Ainf.Pfar: - # for fCP in A2.Ainf.Pfar[lvl]: - # fCP.plot() - # plt.figure() - # diff = np.absolute((A3d[np.ix_(idx, idx)]-A2d)) - # plt.pcolormesh(np.log10(diff)) - # plt.colorbar() - # plt.figure() - # diffRel = np.absolute((A3d[np.ix_(idx, idx)]-A2d)/A2d) - # diffRel[diff < 1e-12] = 0. - # plt.pcolormesh(np.log10(diffRel)) - # print(diffRel[np.isfinite(diffRel)].max(), diffRel[np.isfinite(diffRel)].mean(), np.median(diffRel[np.isfinite(diffRel)])) - # plt.colorbar() - # plt.show() - - assert L2_denseCor/L2_dense < mesh2.h**(0.5+min(kernel2.s.min, 0.5)), (L2_denseCor, L2_dense, L2_denseCor/L2_dense, mesh2.h**(0.5+min(kernel1.s.min, 0.5))) + if buildCorrected: + print('\nCORRECTED\n') + A2.setKernel(kernel2) + + print('\nDENSE\n') + builder3 = nonlocalBuilder(mesh3, dm3, kernel2, zeroExterior=False) + A3 = builder3.getDense() + + A2d = A2.toarray() + A3d = A3.toarray()[np.ix_(idx, idx)] + + nn = np.absolute(A3d) + nn[nn < 1e-16] = 1. + + errDenseCor = np.absolute(A3d-A2d) + errDenseCor_rel = errDenseCor/nn + print('errDenseCor', errDenseCor.max(), errDenseCor_rel.max()) + + y2 = solve(A2d, b2) + x2 = np.zeros((A3.shape[0])) + x2[idx] = y2 + + b3 = dm3.assembleRHS(rhs).toarray() + y3 = solve(A3d, b3[idx]) + x3 = np.zeros((A3.shape[0])) + x3[idx] = y3 + + # assert np.absolute(A3d[np.ix_(idx, idx)]-A2d).max() < 1e-5 + + M = dm3.assembleMass() + + L2_denseCor = np.sqrt(abs(np.vdot(M*(x3-x2), x3-x2))) + L2_dense = np.sqrt(abs(np.vdot(M*x3, x3))) + # L2_cor = np.sqrt(abs(np.vdot(M*x2, x2))) + + print('L2 errDenseCor', L2_denseCor) + + # mesh3.plotFunction(x2, DoFMap=dm3, label='corrected') + # mesh3.plotFunction(x3, DoFMap=dm3, label='dense') + # plt.legend() + # plt.show() + + # if not (L2 < mesh2.h**(0.5+min(kernel2.s.value, 0.5))): + # mesh3.plotFunction(x3, DoFMap=dm3) + # mesh3.plotFunction(x2, DoFMap=dm3) + # plt.figure() + # for lvl in A2.Ainf.Pfar: + # for fCP in A2.Ainf.Pfar[lvl]: + # fCP.plot() + # plt.figure() + # diff = np.absolute((A3d[np.ix_(idx, idx)]-A2d)) + # plt.pcolormesh(np.log10(diff)) + # plt.colorbar() + # plt.figure() + # diffRel = np.absolute((A3d[np.ix_(idx, idx)]-A2d)/A2d) + # diffRel[diff < 1e-12] = 0. + # plt.pcolormesh(np.log10(diffRel)) + # print(diffRel[np.isfinite(diffRel)].max(), diffRel[np.isfinite(diffRel)].mean(), np.median(diffRel[np.isfinite(diffRel)])) + # plt.colorbar() + # plt.show() + + assert L2_denseCor/L2_dense < mesh2.h**(0.5+min(kernel2.s.min, 0.5)), (L2_denseCor, L2_dense, L2_denseCor/L2_dense, mesh2.h**(0.5+min(kernel1.s.min, 0.5))) diff --git a/tests/test_kernels.py b/tests/test_kernels.py index 92525f1..36507dd 100644 --- a/tests/test_kernels.py +++ b/tests/test_kernels.py @@ -7,37 +7,41 @@ import numpy as np from PyNucleus_nl import fractionalOrderFactory, kernelFactory, twoPointFunctionFactory +from PyNucleus_nl.twoPointFunctions import constantTwoPoint from PyNucleus_nl.fractionalOrders import (constFractionalOrder, variableConstFractionalOrder, constantNonSymFractionalOrder, - smoothedLeftRightFractionalOrder) + smoothedLeftRightFractionalOrder, + feFractionalOrder) from PyNucleus_fem.functions import constant -from scipy.special import gamma, erf, digamma +from scipy.special import gamma, erf, digamma, gammaincc from numpy import log from numpy import pi, exp, sqrt from numpy.linalg import norm import pytest -def idfunc(param): - S = [str(p) for p in param] - return '-'.join(S) +def idfuncIntegrable(param): + dim, kernelType, horizon, normalized = param + return f'dim{dim}-kernelType{kernelType}horizon{horizon}-normalized{normalized}' @pytest.fixture(scope='module', params=[ + # 1d kernels (1, 'constant', 0.5, True), (1, 'constant', 0.5, False), - (2, 'constant', 0.5, True), - (2, 'constant', 0.5, False), (1, 'inverseDistance', 0.5, True), (1, 'inverseDistance', 0.5, False), - (2, 'inverseDistance', 0.5, True), - (2, 'inverseDistance', 0.5, False), (1, 'Gaussian', 0.5, True), (1, 'Gaussian', 0.5, False), + # 2d kernels + (2, 'constant', 0.5, True), + (2, 'constant', 0.5, False), + (2, 'inverseDistance', 0.5, True), + (2, 'inverseDistance', 0.5, False), (2, 'Gaussian', 0.5, True), (2, 'Gaussian', 0.5, False), -], ids=idfunc) +], ids=idfuncIntegrable) def integrableKernelParams(request): return request.param @@ -53,6 +57,9 @@ def testIntegrableKernel(integrableKernelParams): else: raise NotImplementedError() kernel = kernelFactory(kernelType, dim=dim, horizon=horizon, normalized=normalized) + infHorizonKernel = kernel.getModifiedKernel(horizon=constant(np.inf)) + boundaryKernelInf = infHorizonKernel.getBoundaryKernel() + horizonValue = kernel.horizon.value if normalized: @@ -84,21 +91,59 @@ def testIntegrableKernel(integrableKernelParams): for x, y in xy_values: + if kernelType == 'constant': + refInf = const + elif kernelType == 'inverseDistance': + refInf = const/norm(x-y) + elif kernelType == 'Gaussian': + invD = (3/horizonValue)**2 + refInf = const*exp(-invD*norm(x-y)**2) + else: + raise NotImplementedError() + if norm(x-y) < horizonValue: - if kernelType == 'constant': - ref = const - elif kernelType == 'inverseDistance': - ref = const/norm(x-y) - elif kernelType == 'Gaussian': - ref = const*exp(-(3.*norm(x-y)/horizonValue)**2) - else: - raise NotImplementedError() + ref = refInf else: ref = 0. - assert np.isclose(kernel(x, y), ref) + if kernelType == 'constant': + refBoundary = refInf*(-1/dim) + elif kernelType == 'inverseDistance': + if dim == 1: + refBoundary = refInf*(-log(norm(x-y))) + else: + refBoundary = refInf*(-1/(dim-1)) + elif kernelType == 'Gaussian': + refBoundary = const*0.5*(invD*norm(x-y)**2)**(-dim/2) * gamma(dim/2) * gammaincc(dim/2, invD*norm(x-y)**2) + else: + raise NotImplementedError() + + # test kernel + assert np.isclose(kernel(x, y), ref), (kernel(x, y), ref) + + # test boundary kernel + assert np.isclose(boundaryKernelInf(x, y), 2*refBoundary*norm(x-y)), (boundaryKernelInf(x, y), 2*refBoundary*norm(x-y)) + + # test that div_y (boundaryKernelInf(x,y) (x-y)/norm(x-y)) == 2*infHorizonKernel(x,y) + eps = 1e-8 + div_fd = 0. + for i in range(dim): + yShifted = y.copy() + yShifted[i] += eps + div_fd += (boundaryKernelInf(x, yShifted) * (x-yShifted)[i]/norm(x-yShifted) - boundaryKernelInf(x, y) * (x-y)[i]/norm(x-y))/eps + assert np.isclose(div_fd, 2*infHorizonKernel(x, y)), (div_fd, 2*infHorizonKernel(x, y)) + + +def idfuncFractional(param): + dim, s, horizon, normalized, phi, derivative = param + return f'dim{dim}-s{s}-horizon{horizon}-normalized{normalized}-phi{phi}-derivative{derivative}' + +from PyNucleus_fem import meshFactory, dofmapFactory +mesh1d = meshFactory('interval', a=-1, b=1, hTarget=1e-2) +dm1d = dofmapFactory('P1', mesh1d, -1) @pytest.fixture(scope='module', params=[ + # 1d kernels (1, fractionalOrderFactory('const', 0.25), np.inf, True, None, 0), (1, fractionalOrderFactory('const', 0.25), np.inf, False, None, 0), (1, fractionalOrderFactory('const', 0.25), 0.5, True, None, 0), @@ -107,35 +152,81 @@ def testIntegrableKernel(integrableKernelParams): (1, fractionalOrderFactory('const', 0.75), np.inf, False, None, 0), (1, fractionalOrderFactory('const', 0.75), 0.5, True, None, 0), (1, fractionalOrderFactory('const', 0.75), 0.5, False, None, 0), + (1, fractionalOrderFactory('varconst', 0.25), np.inf, True, None, 0), + (1, fractionalOrderFactory('varconst', 0.25), np.inf, False, None, 0), + (1, fractionalOrderFactory('varconst', 0.25), 0.5, True, None, 0), + (1, fractionalOrderFactory('varconst', 0.25), 0.5, False, None, 0), (1, fractionalOrderFactory('varconst', 0.75), np.inf, True, None, 0), (1, fractionalOrderFactory('varconst', 0.75), np.inf, False, None, 0), (1, fractionalOrderFactory('varconst', 0.75), 0.5, True, None, 0), (1, fractionalOrderFactory('varconst', 0.75), 0.5, False, None, 0), + (1, fractionalOrderFactory('constantNonSym', 0.25), np.inf, True, None, 0), + (1, fractionalOrderFactory('constantNonSym', 0.25), np.inf, False, None, 0), + (1, fractionalOrderFactory('constantNonSym', 0.25), 0.5, True, None, 0), + (1, fractionalOrderFactory('constantNonSym', 0.25), 0.5, False, None, 0), (1, fractionalOrderFactory('constantNonSym', 0.75), np.inf, True, None, 0), (1, fractionalOrderFactory('constantNonSym', 0.75), np.inf, False, None, 0), (1, fractionalOrderFactory('constantNonSym', 0.75), 0.5, True, None, 0), (1, fractionalOrderFactory('constantNonSym', 0.75), 0.5, False, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), np.inf, True, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), np.inf, False, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), 0.5, True, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), 0.5, False, None, 0), (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, True, None, 0), (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, False, None, 0), (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, True, None, 0), (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, False, None, 0), - (1, fractionalOrderFactory('const', 0.25), np.inf, True, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('const', 0.25), np.inf, False, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('const', 0.25), 0.5, True, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('const', 0.25), 0.5, False, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('const', 0.75), np.inf, True, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('const', 0.75), np.inf, False, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('const', 0.75), 0.5, True, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('const', 0.75), 0.5, False, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('varconst', 0.75), np.inf, True, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('varconst', 0.75), np.inf, False, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('varconst', 0.75), 0.5, True, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('varconst', 0.75), 0.5, False, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, True, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, False, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, True, twoPointFunctionFactory('const', 1.), 0), - (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, False, twoPointFunctionFactory('const', 1.), 0), - # + # discretized fractional order + (1, fractionalOrderFactory('const', 0.25, dm=dm1d), np.inf, True, None, 0), + (1, fractionalOrderFactory('const', 0.25, dm=dm1d), np.inf, False, None, 0), + (1, fractionalOrderFactory('const', 0.25, dm=dm1d), 0.5, True, None, 0), + (1, fractionalOrderFactory('const', 0.25, dm=dm1d), 0.5, False, None, 0), + (1, fractionalOrderFactory('const', 0.75, dm=dm1d), np.inf, True, None, 0), + (1, fractionalOrderFactory('const', 0.75, dm=dm1d), np.inf, False, None, 0), + (1, fractionalOrderFactory('const', 0.75, dm=dm1d), 0.5, True, None, 0), + (1, fractionalOrderFactory('const', 0.75, dm=dm1d), 0.5, False, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25, dm=dm1d), np.inf, True, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25, dm=dm1d), np.inf, False, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25, dm=dm1d), 0.5, True, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25, dm=dm1d), 0.5, False, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75, dm=dm1d), np.inf, True, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75, dm=dm1d), np.inf, False, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75, dm=dm1d), 0.5, True, None, 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75, dm=dm1d), 0.5, False, None, 0), + # now with a trivial phi + (1, fractionalOrderFactory('const', 0.25), np.inf, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('const', 0.25), np.inf, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('const', 0.25), 0.5, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('const', 0.25), 0.5, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('const', 0.75), np.inf, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('const', 0.75), np.inf, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('const', 0.75), 0.5, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('const', 0.75), 0.5, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantSym', 0.25), np.inf, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantSym', 0.25), np.inf, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantSym', 0.25), 0.5, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantSym', 0.25), 0.5, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantSym', 0.75), np.inf, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantSym', 0.75), np.inf, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantSym', 0.75), 0.5, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantSym', 0.75), 0.5, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantNonSym', 0.25), np.inf, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantNonSym', 0.25), np.inf, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantNonSym', 0.25), 0.5, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantNonSym', 0.25), 0.5, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantNonSym', 0.75), np.inf, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantNonSym', 0.75), np.inf, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantNonSym', 0.75), 0.5, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('constantNonSym', 0.75), 0.5, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), np.inf, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), np.inf, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), 0.5, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), 0.5, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, False, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, True, twoPointFunctionFactory('const', 2.), 0), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, False, twoPointFunctionFactory('const', 2.), 0), + # derivative wrt s (1, fractionalOrderFactory('const', 0.25), np.inf, True, None, 1), (1, fractionalOrderFactory('const', 0.25), np.inf, False, None, 1), (1, fractionalOrderFactory('const', 0.25), 0.5, True, None, 1), @@ -144,15 +235,49 @@ def testIntegrableKernel(integrableKernelParams): (1, fractionalOrderFactory('const', 0.75), np.inf, False, None, 1), (1, fractionalOrderFactory('const', 0.75), 0.5, True, None, 1), (1, fractionalOrderFactory('const', 0.75), 0.5, False, None, 1), - (1, fractionalOrderFactory('varconst', 0.75), np.inf, True, None, 1), - (1, fractionalOrderFactory('varconst', 0.75), np.inf, False, None, 1), - (1, fractionalOrderFactory('varconst', 0.75), 0.5, True, None, 1), - (1, fractionalOrderFactory('varconst', 0.75), 0.5, False, None, 1), + (1, fractionalOrderFactory('constantSym', 0.25), np.inf, True, None, 1), + (1, fractionalOrderFactory('constantSym', 0.25), np.inf, False, None, 1), + (1, fractionalOrderFactory('constantSym', 0.25), 0.5, True, None, 1), + (1, fractionalOrderFactory('constantSym', 0.25), 0.5, False, None, 1), + (1, fractionalOrderFactory('constantSym', 0.75), np.inf, True, None, 1), + (1, fractionalOrderFactory('constantSym', 0.75), np.inf, False, None, 1), + (1, fractionalOrderFactory('constantSym', 0.75), 0.5, True, None, 1), + (1, fractionalOrderFactory('constantSym', 0.75), 0.5, False, None, 1), + (1, fractionalOrderFactory('constantNonSym', 0.25), np.inf, True, None, 1), + (1, fractionalOrderFactory('constantNonSym', 0.25), np.inf, False, None, 1), + (1, fractionalOrderFactory('constantNonSym', 0.25), 0.5, True, None, 1), + (1, fractionalOrderFactory('constantNonSym', 0.25), 0.5, False, None, 1), + (1, fractionalOrderFactory('constantNonSym', 0.75), np.inf, True, None, 1), + (1, fractionalOrderFactory('constantNonSym', 0.75), np.inf, False, None, 1), + (1, fractionalOrderFactory('constantNonSym', 0.75), 0.5, True, None, 1), + (1, fractionalOrderFactory('constantNonSym', 0.75), 0.5, False, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), np.inf, True, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), np.inf, False, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), 0.5, True, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), 0.5, False, None, 1), (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, True, None, 1), (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, False, None, 1), (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, True, None, 1), (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, False, None, 1), - # + # discretized fractional order + (1, fractionalOrderFactory('const', 0.25, dm=dm1d), np.inf, True, None, 1), + (1, fractionalOrderFactory('const', 0.25, dm=dm1d), np.inf, False, None, 1), + (1, fractionalOrderFactory('const', 0.25, dm=dm1d), 0.5, True, None, 1), + (1, fractionalOrderFactory('const', 0.25, dm=dm1d), 0.5, False, None, 1), + (1, fractionalOrderFactory('const', 0.75, dm=dm1d), np.inf, True, None, 1), + (1, fractionalOrderFactory('const', 0.75, dm=dm1d), np.inf, False, None, 1), + (1, fractionalOrderFactory('const', 0.75, dm=dm1d), 0.5, True, None, 1), + (1, fractionalOrderFactory('const', 0.75, dm=dm1d), 0.5, False, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25, dm=dm1d), np.inf, True, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25, dm=dm1d), np.inf, False, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25, dm=dm1d), 0.5, True, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25, dm=dm1d), 0.5, False, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75, dm=dm1d), np.inf, True, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75, dm=dm1d), np.inf, False, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75, dm=dm1d), 0.5, True, None, 1), + (1, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75, dm=dm1d), 0.5, False, None, 1), + ################################################## + # 2d kernels (2, fractionalOrderFactory('const', 0.25), np.inf, True, None, 0), (2, fractionalOrderFactory('const', 0.25), np.inf, False, None, 0), (2, fractionalOrderFactory('const', 0.25), 0.5, True, None, 0), @@ -161,15 +286,65 @@ def testIntegrableKernel(integrableKernelParams): (2, fractionalOrderFactory('const', 0.75), np.inf, False, None, 0), (2, fractionalOrderFactory('const', 0.75), 0.5, True, None, 0), (2, fractionalOrderFactory('const', 0.75), 0.5, False, None, 0), - (2, fractionalOrderFactory('varconst', 0.75), np.inf, True, None, 0), - (2, fractionalOrderFactory('varconst', 0.75), np.inf, False, None, 0), - (2, fractionalOrderFactory('varconst', 0.75), 0.5, True, None, 0), - (2, fractionalOrderFactory('varconst', 0.75), 0.5, False, None, 0), + (2, fractionalOrderFactory('constantSym', 0.25), np.inf, True, None, 0), + (2, fractionalOrderFactory('constantSym', 0.25), np.inf, False, None, 0), + (2, fractionalOrderFactory('constantSym', 0.25), 0.5, True, None, 0), + (2, fractionalOrderFactory('constantSym', 0.25), 0.5, False, None, 0), + (2, fractionalOrderFactory('constantSym', 0.75), np.inf, True, None, 0), + (2, fractionalOrderFactory('constantSym', 0.75), np.inf, False, None, 0), + (2, fractionalOrderFactory('constantSym', 0.75), 0.5, True, None, 0), + (2, fractionalOrderFactory('constantSym', 0.75), 0.5, False, None, 0), + (2, fractionalOrderFactory('constantNonSym', 0.25), np.inf, True, None, 0), + (2, fractionalOrderFactory('constantNonSym', 0.25), np.inf, False, None, 0), + (2, fractionalOrderFactory('constantNonSym', 0.25), 0.5, True, None, 0), + (2, fractionalOrderFactory('constantNonSym', 0.25), 0.5, False, None, 0), + (2, fractionalOrderFactory('constantNonSym', 0.75), np.inf, True, None, 0), + (2, fractionalOrderFactory('constantNonSym', 0.75), np.inf, False, None, 0), + (2, fractionalOrderFactory('constantNonSym', 0.75), 0.5, True, None, 0), + (2, fractionalOrderFactory('constantNonSym', 0.75), 0.5, False, None, 0), (2, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, True, None, 0), (2, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, False, None, 0), (2, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, True, None, 0), (2, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, False, None, 0), -], ids=idfunc) + (2, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), np.inf, True, None, 0), + (2, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), np.inf, False, None, 0), + (2, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), 0.5, True, None, 0), + (2, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), 0.5, False, None, 0), + # derivative wrt s + (2, fractionalOrderFactory('const', 0.25), np.inf, True, None, 1), + (2, fractionalOrderFactory('const', 0.25), np.inf, False, None, 1), + (2, fractionalOrderFactory('const', 0.25), 0.5, True, None, 1), + (2, fractionalOrderFactory('const', 0.25), 0.5, False, None, 1), + (2, fractionalOrderFactory('const', 0.75), np.inf, True, None, 1), + (2, fractionalOrderFactory('const', 0.75), np.inf, False, None, 1), + (2, fractionalOrderFactory('const', 0.75), 0.5, True, None, 1), + (2, fractionalOrderFactory('const', 0.75), 0.5, False, None, 1), + (2, fractionalOrderFactory('constantSym', 0.25), np.inf, True, None, 1), + (2, fractionalOrderFactory('constantSym', 0.25), np.inf, False, None, 1), + (2, fractionalOrderFactory('constantSym', 0.25), 0.5, True, None, 1), + (2, fractionalOrderFactory('constantSym', 0.25), 0.5, False, None, 1), + (2, fractionalOrderFactory('constantSym', 0.75), np.inf, True, None, 1), + (2, fractionalOrderFactory('constantSym', 0.75), np.inf, False, None, 1), + (2, fractionalOrderFactory('constantSym', 0.75), 0.5, True, None, 1), + (2, fractionalOrderFactory('constantSym', 0.75), 0.5, False, None, 1), + (2, fractionalOrderFactory('constantNonSym', 0.25), np.inf, True, None, 1), + (2, fractionalOrderFactory('constantNonSym', 0.25), np.inf, False, None, 1), + (2, fractionalOrderFactory('constantNonSym', 0.25), 0.5, True, None, 1), + (2, fractionalOrderFactory('constantNonSym', 0.25), 0.5, False, None, 1), + (2, fractionalOrderFactory('constantNonSym', 0.75), np.inf, True, None, 1), + (2, fractionalOrderFactory('constantNonSym', 0.75), np.inf, False, None, 1), + (2, fractionalOrderFactory('constantNonSym', 0.75), 0.5, True, None, 1), + (2, fractionalOrderFactory('constantNonSym', 0.75), 0.5, False, None, 1), + (2, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), np.inf, True, None, 1), + (2, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), np.inf, False, None, 1), + (2, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), 0.5, True, None, 1), + (2, fractionalOrderFactory('twoDomainNonSym', 0.75, 0.25), 0.5, False, None, 1), + (2, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, True, None, 1), + (2, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), np.inf, False, None, 1), + (2, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, True, None, 1), + (2, fractionalOrderFactory('twoDomainNonSym', 0.25, 0.75), 0.5, False, None, 1), + +], ids=idfuncFractional) def fractionalKernelParams(request): return request.param @@ -186,11 +361,19 @@ def testFractionalKernel(fractionalKernelParams): (np.array([-0.1, 0.1]), np.array([0.5, 0.2]))] else: raise NotImplementedError() - kernel = kernelFactory('fractional', dim=dim, s=s, horizon=horizon, normalized=normalized, phi=phi, derivative=derivative) + + kernel = kernelFactory('fractional', dim=dim, s=s, horizon=horizon, normalized=normalized, phi=phi) + if derivative == 1: + kernel = kernel.getDerivativeKernel() boundaryKernel = kernel.getBoundaryKernel() infHorizonKernel = kernel.getModifiedKernel(horizon=constant(np.inf)) boundaryKernelInf = infHorizonKernel.getBoundaryKernel() + print(kernel) + print(boundaryKernel) + print(infHorizonKernel) + print(boundaryKernelInf) + for x, y in xy_values: sValue = s(x, y) @@ -204,9 +387,11 @@ def testFractionalKernel(fractionalKernelParams): elif isinstance(s, smoothedLeftRightFractionalOrder): assert np.isclose(sValue, s.sFun(x)) if x[0] < 0: - assert np.isclose(sValue, 0.25) + assert np.isclose(sValue, s.sFun.sl) else: - assert np.isclose(sValue, 0.75) + assert np.isclose(sValue, s.sFun.sr) + elif isinstance(s, feFractionalOrder): + pass else: assert False @@ -229,31 +414,86 @@ def testFractionalKernel(fractionalKernelParams): if derivative == 0: refInf = const/norm(x-y)**(dim+2*sValue) * phiValue elif derivative == 1: + refInf = const/norm(x-y)**(dim+2*sValue) * phiValue if normalized: - refInf = const/norm(x-y)**(dim+2*sValue) * phiValue if horizonValue < np.inf: - refInf *= (-2.*log(norm(x-y)/horizonValue) + 1./(sValue-1.)) + refInf *= (-log(norm(x-y)**2/horizonValue**2) - 1./(1.-sValue)) else: - refInf *= (-log(0.25*norm(x-y)**2) + digamma(sValue+0.5*dim) + digamma(-sValue)) + refInf *= (-log(norm(x-y)**2/4) + digamma(sValue+0.5*dim) + digamma(-sValue)) else: - refInf = -2.*0.5/norm(x-y)**(dim+2*sValue) * phiValue * log(norm(x-y)) + refInf *= (-log(norm(x-y)**2)) else: raise NotImplementedError() ref = refInf if (norm(x-y) < horizonValue) else 0. + if derivative == 1: + refInfBoundary = const/norm(x-y)**(dim+2*sValue) * phiValue + if normalized: + if horizonValue == np.inf: + refInfBoundary *= (-log(norm(x-y)**2/4) + digamma(sValue+0.5*dim) + digamma(1.0-sValue)) + else: + refInfBoundary *= (-log(norm(x-y)**2/horizonValue**2) - 1./(1.-sValue) - 1./sValue) + else: + refInfBoundary *= (-log(norm(x-y)**2) - 1.0/sValue) + else: + refInfBoundary = refInf + refBoundary = refInfBoundary if (norm(x-y) < horizonValue) else 0. + # test the kernel with potentially finite horizon assert np.isclose(kernel(x, y), ref) - # test boundary kernel, do not change horizon - np.isclose(boundaryKernel(x, y), ref*norm(x-y)/sValue) - # test kernel with infinite horizon assert np.isclose(infHorizonKernel(x, y), refInf) - if (horizonValue < np.inf) and norm(x-y) > horizonValue: - assert not np.isclose(ref, refInf) - assert np.isclose(ref, 0.) - else: - assert np.isclose(ref, refInf) + if phi is not None and not isinstance(phi, constantTwoPoint): + # Not implemented, as the boundary kernel will depend on phi + return + + # test boundary kernel with potentially finite horizon + assert np.isclose(boundaryKernel(x, y), refBoundary*norm(x-y)/sValue), (boundaryKernel(x, y), refBoundary*norm(x-y)/sValue) + + # test boundary kernel with infinite horizon + assert np.isclose(boundaryKernelInf(x, y), refInfBoundary*norm(x-y)/sValue), (boundaryKernelInf(x, y), refInfBoundary*norm(x-y)/sValue) + + # test that div_y (boundaryKernelInf(x,y) (x-y)/norm(x-y)) == 2*infHorizonKernel(x,y) + eps = 1e-8 + div_fd = 0. + for i in range(dim): + yShifted = y.copy() + yShifted[i] += eps + div_fd += (boundaryKernelInf(x, yShifted) * (x-yShifted)[i]/norm(x-yShifted) - boundaryKernelInf(x, y) * (x-y)[i]/norm(x-y))/eps + assert np.isclose(div_fd, 2*infHorizonKernel(x, y)), (div_fd, 2*infHorizonKernel(x, y)) + + +from PyNucleus import dofmapFactory, fractionalOrderFactory, kernelFactory, meshFactory, nonlocalBuilder, REAL, functionFactory + + +def test_discrete_s_const(): + mesh = meshFactory('interval', a=-1, b=1, hTarget=1e-2) + dmS = dofmapFactory('P1', mesh, -1) + sFun = fractionalOrderFactory('constantNonSym', 0.75) + sFE = fractionalOrderFactory('constantNonSym', 0.75, dm=dmS) + kernelFun = kernelFactory('fractional', dim=mesh.dim, s=sFun) + kernelFE = kernelFactory('fractional', dim=mesh.dim, s=sFE) + dm = dofmapFactory('P1', mesh, 0) + A_fun = dm.assembleNonlocal(kernelFun) + A_fe = dm.assembleNonlocal(kernelFE) + assert np.absolute((A_fun-A_fe).toarray()/A_fun.toarray()).max() == 0. + - # test boundary kernel, infinite horizon - assert np.isclose(boundaryKernelInf(x, y), refInf*norm(x-y)/sValue) +def test_discrete_leftRight(): + mesh = meshFactory('interval', a=-1, b=1, hTarget=1e-2) + dmS = dofmapFactory('P1', mesh, -1) + sFun = fractionalOrderFactory('twoDomainNonSym', sl=0.25, sr=0.75, r=0.3) + sFE = fractionalOrderFactory('twoDomainNonSym', sl=0.25, sr=0.75, r=0.3, dm=dmS) + kernelFun = kernelFactory('fractional', dim=mesh.dim, s=sFun) + kernelFE = kernelFactory('fractional', dim=mesh.dim, s=sFE) + dm = dofmapFactory('P1', mesh, 0) + A_fun = dm.assembleNonlocal(kernelFun) + A_fe = dm.assembleNonlocal(kernelFE) + print('max rel error', np.absolute((A_fun-A_fe).toarray()/A_fun.toarray()).max()) + x = dm.fromArray(np.random.randn(dm.num_dofs)) + a = x.inner(A_fun*x) + b = x.inner(A_fe*x) + print('apply abs error', abs(a-b)) + print('apply rel error', abs(a-b)/a) + assert abs(a-b)/a < 1e-4