Skip to content

Commit e238042

Browse files
cgcgcgcgcgcg
committed
updates
1 parent 5000cc5 commit e238042

24 files changed

+2030
-754
lines changed

‎base/PyNucleus_base/utilsFem.py

Lines changed: 8 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -427,14 +427,14 @@ def getMPIinfo(grp):
427427
if MPI.COMM_WORLD.rank == 0:
428428
hosts = ','.join(set(hosts))
429429
grp.add('MPI library', '{}'.format(MPI.Get_library_version()[:-1]))
430-
for label, value in [('MPI standard supported:', MPI.Get_version()),
431-
('Vendor:', MPI.get_vendor()),
432-
('Level of thread support:', t[MPI.Query_thread()]),
433-
('Is threaded:', MPI.Is_thread_main()),
434-
('Threads requested:', mpi4py.rc.threads),
435-
('Thread level requested:', mpi4py.rc.thread_level),
436-
('Hosts:', hosts),
437-
('Communicator size:', MPI.COMM_WORLD.size)]:
430+
for label, value in [('MPI standard supported', MPI.Get_version()),
431+
('Vendor', MPI.get_vendor()),
432+
('Level of thread support', t[MPI.Query_thread()]),
433+
('Is threaded', MPI.Is_thread_main()),
434+
('Threads requested', mpi4py.rc.threads),
435+
('Thread level requested', mpi4py.rc.thread_level),
436+
('Hosts', hosts),
437+
('Communicator size', MPI.COMM_WORLD.size)]:
438438
grp.add(label, value)
439439

440440

‎drivers/runNonlocal.py

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,7 @@
77
###################################################################################
88

99
from mpi4py import MPI
10-
from PyNucleus import driver, DIRICHLET, NEUMANN
10+
from PyNucleus import driver, NEUMANN
1111
from PyNucleus.nl import (nonlocalPoissonProblem,
1212
discretizedNonlocalProblem)
1313

‎fem/PyNucleus_fem/factories.py

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,8 @@
1212
from . functions import (_rhsFunSin1D, _solSin1D, _rhsFunSin2D, _cos1D, _cos2D, _rhsCos2D, _grad_cos2d_n,
1313
_rhsFunSin3D, _solSin2D, _solSin3D, Lambda, constant,
1414
monomial,
15+
affineFunction,
16+
sqrtAffineFunction,
1517
complexLambda,
1618
real, imag,
1719
_rhsFunSin3D_memoized,
@@ -82,6 +84,8 @@ def rhsFractional2D_nonPeriodic(s):
8284
functionFactory.register('rhsFractional2D', rhsFractional2D)
8385
functionFactory.register('constant', constant)
8486
functionFactory.register('monomial', monomial)
87+
functionFactory.register('affine', affineFunction)
88+
functionFactory.register('sqrt_affine', sqrtAffineFunction)
8589
functionFactory.register('x0', monomial, params={'exponent': np.array([1., 0., 0.])})
8690
functionFactory.register('x1', monomial, params={'exponent': np.array([0., 1., 0.])})
8791
functionFactory.register('x2', monomial, params={'exponent': np.array([0., 0., 1.])})

‎fem/PyNucleus_fem/functions.pxd

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -11,13 +11,26 @@ from numpy cimport uint8_t
1111

1212
cdef class function:
1313
cdef REAL_t eval(self, REAL_t[::1] x)
14+
cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x)
1415

1516

1617
cdef class constant(function):
1718
cdef:
1819
public REAL_t value
1920

2021

22+
cdef class affineFunction(function):
23+
cdef:
24+
public REAL_t[::1] w
25+
public REAL_t c
26+
27+
28+
cdef class sqrtAffineFunction(function):
29+
cdef:
30+
public REAL_t[::1] w
31+
public REAL_t c
32+
33+
2134
ctypedef REAL_t(*volume_t)(REAL_t[:, ::1])
2235

2336

@@ -43,3 +56,9 @@ cdef class matrixFunction:
4356
INDEX_t columns
4457
public BOOL_t symmetric
4558
cdef void eval(self, REAL_t[::1] x, REAL_t[:, ::1] vals)
59+
cdef void evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* vals)
60+
61+
62+
cdef class constantMatrixFunction(matrixFunction):
63+
cdef:
64+
REAL_t[:, ::1] A

‎fem/PyNucleus_fem/functions.pyx

Lines changed: 110 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -26,6 +26,11 @@ cdef class function:
2626
cdef REAL_t eval(self, REAL_t[::1] x):
2727
pass
2828

29+
cdef REAL_t evalPtr(self, INDEX_t dim, REAL_t* x):
30+
cdef:
31+
REAL_t[::1] xA = <REAL_t[:dim]> x
32+
return self.eval(xA)
33+
2934
def __add__(self, function other):
3035
if isinstance(self, mulFunction):
3136
if isinstance(other, mulFunction):
@@ -221,6 +226,75 @@ cdef class monomial(function):
221226
else:
222227
return s
223228

229+
def __reduce__(self):
230+
return monomial, (np.array(self.exponent, copy=True), self.factor)
231+
232+
233+
cdef class affineFunction(function):
234+
def __init__(self, REAL_t[::1] w, REAL_t c):
235+
self.w = w
236+
self.c = c
237+
238+
cdef inline REAL_t eval(self, REAL_t[::1] x):
239+
cdef:
240+
INDEX_t i
241+
REAL_t s = self.c
242+
for i in range(x.shape[0]):
243+
s += x[i]*self.w[i]
244+
return s
245+
246+
def __repr__(self):
247+
cdef:
248+
INDEX_t i
249+
s = ''
250+
for i in range(self.w.shape[0]):
251+
if self.w[i] != 0.:
252+
if len(s) > 0:
253+
s += '+'
254+
if self.w[i] != 1.:
255+
s += '{}*x_{}'.format(self.w[i], i)
256+
else:
257+
s += 'x_{}'.format(i)
258+
if self.c != 0.:
259+
s += '+'+str(self.c)
260+
return s
261+
262+
def __reduce__(self):
263+
return affineFunction, (np.array(self.w, copy=True), self.c)
264+
265+
266+
cdef class sqrtAffineFunction(function):
267+
def __init__(self, REAL_t[::1] w, REAL_t c):
268+
self.w = w
269+
self.c = c
270+
271+
cdef inline REAL_t eval(self, REAL_t[::1] x):
272+
cdef:
273+
INDEX_t i
274+
REAL_t s = self.c
275+
for i in range(x.shape[0]):
276+
s += x[i]*self.w[i]
277+
return sqrt(s)
278+
279+
def __repr__(self):
280+
cdef:
281+
INDEX_t i
282+
s = ''
283+
for i in range(self.w.shape[0]):
284+
if self.w[i] != 0.:
285+
if len(s) > 0:
286+
s += '+'
287+
if self.w[i] != 1.:
288+
s += '{}*x_{}'.format(self.w[i], i)
289+
else:
290+
s += 'x_{}'.format(i)
291+
if self.c != 0.:
292+
s += '+'+str(self.c)
293+
return 'sqrt('+s+')'
294+
295+
def __reduce__(self):
296+
return affineFunction, (np.array(self.w, copy=True), self.c)
297+
224298

225299
cdef class _rhsFunSin1D(function):
226300
cdef inline REAL_t eval(self, REAL_t[::1] x):
@@ -2144,6 +2218,21 @@ cdef class matrixFunction:
21442218
f = self.components[i][j]
21452219
vals[i, j] = f.eval(x)
21462220

2221+
cdef void evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* vals):
2222+
cdef:
2223+
INDEX_t i
2224+
function f
2225+
if self.symmetric:
2226+
for i in range(self.rows):
2227+
for j in range(i, self.columns):
2228+
f = self.components[i][j]
2229+
vals[i*self.columns+j] = vals[j*self.columns+i] = f.evalPtr(dim, x)
2230+
else:
2231+
for i in range(self.rows):
2232+
for j in range(self.columns):
2233+
f = self.components[i][j]
2234+
vals[i*self.columns+j] = f.evalPtr(dim, x)
2235+
21472236
def __repr__(self):
21482237
return '{}({})'.format(self.__class__.__name__, ','.join([f.__repr__() for f in self.components]))
21492238

@@ -2152,6 +2241,27 @@ cdef class matrixFunction:
21522241
return self.components[i][j]
21532242

21542243

2244+
cdef class constantMatrixFunction(matrixFunction):
2245+
def __init__(self, REAL_t[:, ::1] A):
2246+
self.A = A
2247+
self.rows = A.shape[0]
2248+
self.columns = A.shape[1]
2249+
2250+
cdef void eval(self, REAL_t[::1] x, REAL_t[:, ::1] vals):
2251+
cdef:
2252+
INDEX_t i, j
2253+
for i in range(self.rows):
2254+
for j in range(self.columns):
2255+
vals[i, j] = self.A[i, j]
2256+
2257+
cdef void evalPtr(self, INDEX_t dim, REAL_t* x, REAL_t* vals):
2258+
cdef:
2259+
INDEX_t i, j, k = 0
2260+
for i in range(self.rows):
2261+
for j in range(self.columns):
2262+
vals[k] = self.A[i, j]
2263+
k += 1
2264+
21552265
cdef class periodicityFunctor(function):
21562266
cdef:
21572267
function f

‎fem/PyNucleus_fem/mesh.py

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -344,7 +344,7 @@ def squareWithInteractions(ax, ay, bx, by,
344344
(lineHorizontal+bottomLeft) +
345345
(lineHorizontal+(bottomLeft-verticalOffset)))
346346

347-
d2 = (circularSegment(bottomRight, horizon, -0.5*np.pi, 0., numPointsPerUnitLength) +
347+
d2 = (circularSegment(bottomRight, horizon, 1.5*np.pi, 2.*np.pi, numPointsPerUnitLength) +
348348
line(bottomRight, bottomRight+horizontalOffset) +
349349
line(bottomRight, bottomRight-verticalOffset) +
350350
(lineVertical+(bottomRight+horizontalOffset)) +
@@ -395,7 +395,6 @@ def squareWithInteractions(ax, ay, bx, by,
395395
xVals4 = np.sort(mesh.vertices_as_array[idx4, 0])
396396
assert np.allclose(xVals3, xVals4), (xVals3, xVals4)
397397
mesh2 = uniformSquare(ax=ax, ay=ay, bx=bx, by=by, xVals=xVals3, yVals=yVals1)
398-
mesh2.plot()
399398
mesh = snapMeshes(mesh, mesh2)
400399

401400
location = uninitialized((mesh.num_vertices), dtype=INDEX)

‎fem/PyNucleus_fem/vector_{SCALAR}.pxi

Lines changed: 8 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -251,14 +251,12 @@ cdef class {SCALAR_label_lc_}fe_vector:
251251
assign(v.data, self.data)
252252
return v
253253

254-
def __getstate__(self):
255-
return (np.array(self.data, copy=False), self.dm)
256-
257-
def __setstate__(self, state):
258-
self.data = state[0]
259-
self.dm = state[1]
254+
def __reduce__(self):
255+
return {SCALAR_label_lc_}fe_vector, (np.array(self.data), self.dm)
260256

261257
def __getattr__(self, name):
258+
if name == '__deepcopy__':
259+
raise AttributeError()
262260
return getattr(np.array(self.data, copy=False), name)
263261

264262
cpdef REAL_t norm(self, BOOL_t acc=False, BOOL_t asynchronous=False):
@@ -597,14 +595,12 @@ cdef class {SCALAR_label_lc_}multi_fe_vector:
597595
{SCALAR_label_lc_}assign_2d(v.data, self.data)
598596
return v
599597

600-
def __getstate__(self):
601-
return (np.array(self.data, copy=False), self.dm)
602-
603-
def __setstate__(self, state):
604-
self.data = state[0]
605-
self.dm = state[1]
598+
def __reduce__(self):
599+
return {SCALAR_label_lc_}multi_fe_vector, (np.array(self.data), self.dm)
606600

607601
def __getattr__(self, name):
602+
if name == '__deepcopy__':
603+
raise AttributeError()
608604
return getattr(np.array(self.data, copy=False), name)
609605

610606
def linearPart(self):

‎nl/PyNucleus_nl/discretizedProblems.py

Lines changed: 12 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -7,7 +7,6 @@
77

88
import numpy as np
99
from PyNucleus_base import solverFactory
10-
from PyNucleus_base.performanceLogger import PLogger
1110
from PyNucleus_base.utilsFem import (classWithComputedDependencies,
1211
problem,
1312
generates)
@@ -358,7 +357,8 @@ def __init__(self, driver, continuumProblem):
358357

359358
def setDriverArgs(self):
360359
p = self.driver.addGroup('solver')
361-
self.setDriverFlag('solverType', acceptedValues=['cg-mg', 'gmres-mg', 'lu', 'mg', 'cg-jacobi', 'gmres-jacobi'], help='solver for the linear system', group=p)
360+
self.setDriverFlag('solverType', acceptedValues=['cg-mg', 'gmres-mg', 'lu', 'mg',
361+
'cg-jacobi', 'gmres-jacobi'], help='solver for the linear system', group=p)
362362
self.setDriverFlag('maxiter', 100, help='maximum number of iterations', group=p)
363363
self.setDriverFlag('tol', 1e-6, help='solver tolerance', group=p)
364364

@@ -483,7 +483,8 @@ def buildHierarchy(self,
483483

484484
self.hierarchy = hierarchy
485485
if kernel is not None:
486-
assert 2*self.finalMesh.h < kernel.horizon.value, "Please choose horizon bigger than two mesh sizes. h = {}, horizon = {}".format(self.finalMesh.h, kernel.horizon.value)
486+
assert 2*self.finalMesh.h < kernel.max_horizon, ("Please choose horizon bigger than two mesh sizes. " +
487+
"h = {}, horizon = {}").format(self.finalMesh.h, kernel.horizon.value)
487488

488489
@generates('adjointHierarchy')
489490
def buildAdjointHierarchy(self, hierarchy):
@@ -522,7 +523,8 @@ def buildBCoperator(self, dmInterior, dmBC,
522523
assemblyParams['dense'] = matrixFormat == 'dense'
523524
assemblyParams['matrixFormat'] = matrixFormat
524525
assemblyParams['tag'] = tag
525-
self.A_BC = getFracLapl(dmInterior, dm2=dmBC, **assemblyParams)
526+
with self.timer('build BC operator'):
527+
self.A_BC = getFracLapl(dmInterior, dm2=dmBC, **assemblyParams)
526528
else:
527529
self.A_BC = None
528530

@@ -540,8 +542,9 @@ def getOperators(self, hierarchy):
540542

541543
@generates('A_derivative')
542544
def getDerivativeOperator(self, kernel, dmInterior, matrixFormat, eta, target_order):
543-
self.A_derivative = dmInterior.assembleNonlocal(kernel.getDerivativeKernel(derivative=1), matrixFormat=matrixFormat, params={'eta': eta,
544-
'target_order': target_order})
545+
self.A_derivative = dmInterior.assembleNonlocal(kernel.getDerivativeKernel(derivative=1),
546+
matrixFormat=matrixFormat, params={'eta': eta,
547+
'target_order': target_order})
545548

546549
@generates('b')
547550
def buildRHS(self, rhs, dim, A_BC, dmBC, dirichletData, boundaryCondition, solverType, dmInterior, hierarchy):
@@ -669,7 +672,7 @@ def adjointSolve(self, b, dm, dmInterior, P_interior, adjointSolver, tol, maxite
669672
self.adjointModelSolution = stationaryModelSolution(self, u, **data)
670673

671674
def report(self, group):
672-
group.add('kernel',self.continuumProblem.kernel)
675+
group.add('kernel', repr(self.continuumProblem.kernel))
673676
group.add('h', self.finalMesh.h)
674677
group.add('hmin', self.finalMesh.hmin)
675678
if self.continuumProblem.kernel is not None:
@@ -855,6 +858,8 @@ def setInitialCondition(self, dm, initial):
855858

856859
if self.doMovie:
857860
if self._driver.isMaster:
861+
from PyNucleus_base.plot_utils import movieCreator
862+
858863
outputFolder = self.movieFolder
859864
movie_kwargs = {}
860865
if self.continuumProblem.dim == 2:

‎nl/PyNucleus_nl/fractionalLaplacian1D.pyx

Lines changed: 13 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -5,7 +5,7 @@
55
# If you want to use this code, please refer to the README.rst and LICENSE files. #
66
###################################################################################
77

8-
from libc.math cimport (log, ceil, fabs as abs, pow)
8+
from libc.math cimport log, ceil, fabs as abs
99
import numpy as np
1010
cimport numpy as np
1111

@@ -415,6 +415,18 @@ cdef class fractionalLaplacian1D_nonsym(fractionalLaplacian1D):
415415
416416
for the non-symmetric 1D nonlocal Laplacian.
417417
"""
418+
def __init__(self,
419+
Kernel kernel,
420+
meshBase mesh,
421+
DoFMap dm,
422+
target_order=None,
423+
quad_order_diagonal=None,
424+
num_dofs=None,
425+
**kwargs):
426+
super(fractionalLaplacian1D_nonsym, self).__init__(kernel, mesh, dm, num_dofs, **kwargs)
427+
self.symmetricLocalMatrix = False
428+
self.symmetricCells = False
429+
418430
cdef panelType getQuadOrder(self,
419431
const REAL_t h1,
420432
const REAL_t h2,

‎nl/PyNucleus_nl/fractionalLaplacian2D.pyx

Lines changed: 12 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -897,6 +897,18 @@ cdef class fractionalLaplacian2D_nonsym(fractionalLaplacian2D):
897897
898898
for the 2D non-symmetric nonlocal Laplacian.
899899
"""
900+
def __init__(self,
901+
Kernel kernel,
902+
meshBase mesh,
903+
DoFMap dm,
904+
target_order=None,
905+
quad_order_diagonal=None,
906+
num_dofs=None,
907+
**kwargs):
908+
super(fractionalLaplacian2D_nonsym, self).__init__(kernel, mesh, dm, num_dofs, **kwargs)
909+
self.symmetricLocalMatrix = False
910+
self.symmetricCells = False
911+
900912
cdef panelType getQuadOrder(self,
901913
const REAL_t h1,
902914
const REAL_t h2,

0 commit comments

Comments
 (0)