From 1220b188a42e6f255cb7ad388d7deaf3bca3b082 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 11 Nov 2019 20:18:44 +0100 Subject: [PATCH 01/23] Add Said's fixes for boundary assembly. --- psydac/api/ast/basic.py | 5 +- psydac/api/ast/evaluation.py | 34 ++++--- psydac/api/ast/fem.py | 168 ++++++++++++++++++----------------- psydac/api/ast/utilities.py | 68 +++++++------- psydac/api/basic.py | 10 +-- psydac/api/fem.py | 12 +++ psydac/api/grid.py | 44 +++++++-- 7 files changed, 191 insertions(+), 150 deletions(-) diff --git a/psydac/api/ast/basic.py b/psydac/api/ast/basic.py index e614b0bb1..a949ab170 100644 --- a/psydac/api/ast/basic.py +++ b/psydac/api/ast/basic.py @@ -6,7 +6,6 @@ #============================================================================== class SplBasic(Basic): - _discrete_boundary = None def __new__(cls, tag, name=None, prefix=None, debug=False, detailed=False, mapping=None, is_rational_mapping=None): @@ -67,8 +66,8 @@ def is_rational_mapping(self): return self._is_rational_mapping @property - def discrete_boundary(self): - return self._discrete_boundary + def boundary(self): + return self._boundary @property def imports(self): diff --git a/psydac/api/ast/evaluation.py b/psydac/api/ast/evaluation.py index 3a82d20a0..aa5b14a88 100644 --- a/psydac/api/ast/evaluation.py +++ b/psydac/api/ast/evaluation.py @@ -23,7 +23,7 @@ from .basic import SplBasic from .utilities import build_pythran_types_header, variables -from .utilities import filter_loops, filter_product +from .utilities import filter_loops, filter_product, select_loops from .utilities import rationalize_eval_mapping from .utilities import compute_atoms_expr_mapping from .utilities import compute_atoms_expr_field @@ -31,7 +31,7 @@ #============================================================================== class EvalQuadratureMapping(SplBasic): - def __new__(cls, space, mapping, discrete_boundary=None, name=None, + def __new__(cls, space, mapping, boundary=None, name=None, boundary_basis=None, nderiv=1, is_rational_mapping=None, area=None, backend=None): @@ -58,7 +58,7 @@ def __new__(cls, space, mapping, discrete_boundary=None, name=None, is_rational_mapping=is_rational_mapping) obj._space = space - obj._discrete_boundary = discrete_boundary + obj._boundary = boundary obj._boundary_basis = boundary_basis obj._backend = backend @@ -254,9 +254,7 @@ def _initialize(self): # ... # put the body in tests for loops - body = filter_loops(indices_basis, ranges_basis, body, - self.discrete_boundary, - boundary_basis=self.boundary_basis) + body = select_loops(indices_basis, ranges_basis, body, boundary=None) if self.is_rational_mapping: stmts = rationalize_eval_mapping(self.mapping, self.nderiv, @@ -266,7 +264,7 @@ def _initialize(self): # ... if self.area: - weight = filter_product(indices_quad, weights, self.discrete_boundary) + weight = filter_product(indices_quad, weights, self.boundary) stmts = area_eval_mapping(self.mapping, self.area, dim, indices_quad, weight) @@ -275,7 +273,7 @@ def _initialize(self): # put the body in for loops of quadrature points body = filter_loops(indices_quad, ranges_quad, body, - self.discrete_boundary, + self.boundary, boundary_basis=self.boundary_basis) # initialization of the matrix @@ -308,7 +306,7 @@ def _initialize(self): #============================================================================== class EvalQuadratureField(SplBasic): - def __new__(cls, space, fields, discrete_boundary=None, name=None, + def __new__(cls, space, fields, boundary=None, name=None, boundary_basis=None, mapping=None, is_rational_mapping=None,backend=None): if not isinstance(fields, (tuple, list, Tuple)): @@ -320,7 +318,7 @@ def __new__(cls, space, fields, discrete_boundary=None, name=None, obj._space = space obj._fields = Tuple(*fields) - obj._discrete_boundary = discrete_boundary + obj._boundary = boundary obj._boundary_basis = boundary_basis obj._backend = backend obj._func = obj._initialize() @@ -412,13 +410,13 @@ def _initialize(self): # put the body in tests for loops body = filter_loops(indices_basis, ranges_basis, body, - self.discrete_boundary, + self.boundary, boundary_basis=self.boundary_basis) # put the body in for loops of quadrature points body = filter_loops(indices_quad, ranges_quad, body, - self.discrete_boundary, + self.boundary, boundary_basis=self.boundary_basis) @@ -446,7 +444,7 @@ def _initialize(self): #============================================================================== class EvalQuadratureVectorField(SplBasic): - def __new__(cls, space, vector_fields, discrete_boundary=None, name=None, + def __new__(cls, space, vector_fields, boundary=None, name=None, boundary_basis=None, mapping=None, is_rational_mapping=None, backend = None): if not isinstance(vector_fields, (tuple, list, Tuple)): @@ -458,7 +456,7 @@ def __new__(cls, space, vector_fields, discrete_boundary=None, name=None, obj._space = space obj._vector_fields = Tuple(*vector_fields) - obj._discrete_boundary = discrete_boundary + obj._boundary = boundary obj._boundary_basis = boundary_basis obj._backend = backend obj._func = obj._initialize() @@ -548,12 +546,12 @@ def _initialize(self): # put the body in tests for loops body = filter_loops(indices_basis, ranges_basis, body, - self.discrete_boundary, + self.boundary, boundary_basis=self.boundary_basis) # put the body in for loops of quadrature points body = filter_loops(indices_quad, ranges_quad, body, - self.discrete_boundary, + self.boundary, boundary_basis=self.boundary_basis) # initialization of the matrix @@ -597,7 +595,7 @@ def _create_loop(indices, ranges, body): # NOTE: this is used in module 'psydac.api.ast.glt' class EvalArrayField(SplBasic): - def __new__(cls, space, fields, discrete_boundary=None, name=None, + def __new__(cls, space, fields, boundary=None, name=None, boundary_basis=None, mapping=None, is_rational_mapping=None,backend=None): if not isinstance(fields, (tuple, list, Tuple)): @@ -609,7 +607,7 @@ def __new__(cls, space, fields, discrete_boundary=None, name=None, obj._space = space obj._fields = Tuple(*fields) - obj._discrete_boundary = discrete_boundary + obj._boundary = boundary obj._boundary_basis = boundary_basis obj._backend = backend obj._func = obj._initialize() diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py index c0c7521d8..9d8098fa9 100644 --- a/psydac/api/ast/fem.py +++ b/psydac/api/ast/fem.py @@ -64,15 +64,16 @@ FunctionalForms = (BilinearForm, LinearForm, Functional) #============================================================================== -def init_loop_quadrature(indices, ranges, discrete_boundary): +def init_loop_quadrature(indices, ranges, boundary): stmts = [] - if not discrete_boundary: + if not boundary: return stmts - # TODO improve using namedtuple or a specific class ? to avoid the 0 index - # => make it easier to understand - quad_mask = [i[0] for i in discrete_boundary] - quad_ext = [i[1] for i in discrete_boundary] + if isinstance(boundary, Boundary): + quad_mask = [boundary.axis] + quad_ext = [boundary.ext] + else: + raise TypeError(boundary) dim = len(indices) for i in range(dim-1,-1,-1): @@ -88,15 +89,16 @@ def init_loop_quadrature(indices, ranges, discrete_boundary): return stmts #============================================================================== -def init_loop_basis(indices, ranges, discrete_boundary): +def init_loop_basis(indices, ranges, boundary): stmts = [] - if not discrete_boundary: + if not boundary: return stmts - # TODO improve using namedtuple or a specific class ? to avoid the 0 index - # => make it easier to understand - quad_mask = [i[0] for i in discrete_boundary] - quad_ext = [i[1] for i in discrete_boundary] + if isinstance(boundary, Boundary): + quad_mask = [boundary.axis] + quad_ext = [boundary.ext] + else: + raise TypeError(boundary) dim = len(indices) for i in range(dim-1,-1,-1): @@ -124,15 +126,16 @@ def init_loop_support(indices_elm, n_elements, weights_in_elm, weights, test_basis_in_elm, test_basis, trial_basis_in_elm, trial_basis, - is_bilinear, discrete_boundary): + is_bilinear, boundary): stmts = [] - if not discrete_boundary: + if not boundary: return stmts - # TODO improve using namedtuple or a specific class ? to avoid the 0 index - # => make it easier to understand - quad_mask = [i[0] for i in discrete_boundary] - quad_ext = [i[1] for i in discrete_boundary] + if isinstance(boundary, Boundary): + quad_mask = [boundary.axis] + quad_ext = [boundary.ext] + else: + raise TypeError(boundary) dim = len(indices_elm) for i in range(dim-1,-1,-1): @@ -207,7 +210,7 @@ def area_eval_mapping(mapping, area, dim, indices_quad, weight): class Kernel(SplBasic): def __new__(cls, weak_form, kernel_expr, target=None, - discrete_boundary=None, name=None, boundary_basis=None, + boundary=None, name=None, boundary_basis=None, mapping=None, is_rational_mapping=None,symbolic_space=None, backend=None): if not isinstance(weak_form, FunctionalForms): @@ -244,19 +247,15 @@ def __new__(cls, weak_form, kernel_expr, target=None, # ... # ... - if discrete_boundary: - if not isinstance(discrete_boundary, (tuple, list)): - raise TypeError('> Expecting a tuple or list for discrete_boundary') - - discrete_boundary = list(discrete_boundary) - if not isinstance(discrete_boundary[0], (tuple, list)): - discrete_boundary = [discrete_boundary] - # discrete_boundary is now a list of lists + if boundary: + if not isinstance(boundary, Boundary): + print(type(boundary)) + raise TypeError('> Expecting a Boundary for boundary') # ... - # ... discrete_boundary must be given if there are Trace nodes - if on_boundary and not discrete_boundary: - raise ValueError('> discrete_bounary must be provided for a boundary Kernel') + # ... boundary must be given if there are Trace nodes + if on_boundary and not boundary: + raise ValueError('> boundary must be provided for a boundary Kernel') # ... # ... default value for boundary_basis is True if on boundary @@ -272,7 +271,7 @@ def __new__(cls, weak_form, kernel_expr, target=None, obj._weak_form = weak_form obj._kernel_expr = kernel_expr obj._target = target - obj._discrete_boundary = discrete_boundary + obj._boundary = boundary obj._boundary_basis = boundary_basis obj._area = None obj._user_functions = [] @@ -518,9 +517,11 @@ def _initialize(self): space = list(fs)[0].space eval_field = EvalQuadratureField(space, fields_expressions, - discrete_boundary=self.discrete_boundary, - boundary_basis=self.boundary_basis, - mapping=mapping,backend=self.backend) + boundary = self.boundary, + boundary_basis = self.boundary_basis, + mapping = mapping, + backend = self.backend) + fields += list(eval_field.fields) self._eval_fields.append(eval_field) for k,v in eval_field.map_stmts.items(): @@ -555,9 +556,11 @@ def _initialize(self): space = list(fs)[0].space eval_vector_field = EvalQuadratureVectorField(space, vector_fields_expressions, - discrete_boundary=self.discrete_boundary, - boundary_basis=self.boundary_basis, - mapping=mapping,backend=self.backend) + boundary = self.boundary, + boundary_basis = self.boundary_basis, + mapping = mapping, + backend = self.backend) + vector_fields += list(eval_vector_field.vector_fields) self._eval_vector_fields.append(eval_vector_field) for k,v in eval_vector_field.map_stmts.items(): @@ -599,12 +602,13 @@ def _initialize(self): space = self.weak_form.space eval_mapping = EvalQuadratureMapping(space, mapping, - discrete_boundary=self.discrete_boundary, - boundary_basis=self.boundary_basis, - nderiv=nderiv, - is_rational_mapping=self.is_rational_mapping, - area=self.area, - backend=self.backend) + boundary = self.boundary, + boundary_basis = self.boundary_basis, + nderiv = nderiv, + is_rational_mapping = self.is_rational_mapping, + area = self.area, + backend = self.backend) + self._eval_mapping = eval_mapping # update dependencies @@ -821,7 +825,7 @@ def _initialize(self): expr = [e.subs(vector[i], normal_vec[i]) for e in expr] map_stmts, stmts = compute_normal_vector(normal_vec, - self.discrete_boundary, + self.boundary, mapping) elif isinstance(vector, TangentVector): @@ -830,7 +834,7 @@ def _initialize(self): expr = [e.subs(vector[i], tangent_vec[i]) for e in expr] map_stmts, stmts = compute_tangent_vector(tangent_vec, - self.discrete_boundary, + self.boundary, mapping) for stmt in map_stmts: @@ -847,23 +851,24 @@ def _initialize(self): inv_jac = Symbol('inv_jac') det_jac = Symbol('det_jac') - if not isinstance(self.target, Boundary): + if isinstance(self.target, Boundary): + det_jac_bnd = Symbol('det_jac_bnd') - # ... inv jacobian - jac = mapping.det_jacobian - jac = SymbolicExpr(jac) - # ... + # ... inv jacobian + jac = mapping.det_jacobian + jac = SymbolicExpr(jac) + # ... - body += [Assign(det_jac, jac)] - body += [Assign(inv_jac, 1./jac)] + body += [Assign(det_jac, jac)] + body += [Assign(inv_jac, 1./jac)] - # TODO do we use the same inv_jac? - # if not isinstance(self.target, Boundary): - # body += [Assign(inv_jac, 1/jac)] + # TODO do we use the same inv_jac? + # if not isinstance(self.target, Boundary): + # body += [Assign(inv_jac, 1/jac)] - init_map = OrderedDict(sorted(init_map.items())) - for stmt in list(init_map.values()): - body += [stmt.subs(1/jac, inv_jac)] + init_map = OrderedDict(sorted(init_map.items())) + for stmt in list(init_map.values()): + body += [stmt.subs(1/jac, inv_jac)] else: body += [Assign(coordinates[i],positions[i][indices_quad[i]]) @@ -871,7 +876,7 @@ def _initialize(self): # ... # ... - weighted_vol = filter_product(indices_quad, weighted_vols, self.discrete_boundary) + weighted_vol = filter_product(indices_quad, weighted_vols, self.boundary) # ... # ... @@ -909,7 +914,10 @@ def _initialize(self): # TODO use positive mapping all the time? Abs? if mapping: - weighted_vol = weighted_vol * Abs(det_jac) + if isinstance(self.target, Boundary): + weighted_vol = weighted_vol * Abs(det_jac_bnd) + else: + weighted_vol = weighted_vol * Abs(det_jac) body.append(Assign(wvol,weighted_vol)) @@ -926,11 +934,9 @@ def _initialize(self): # ... # put the body in for loops of quadrature points - init_stmts += init_loop_quadrature( indices_quad, ranges_quad, - self.discrete_boundary ) + init_stmts += init_loop_quadrature( indices_quad, ranges_quad, self.boundary ) - body = select_loops( indices_quad, ranges_quad, body, - self.discrete_boundary, + body = select_loops( indices_quad, ranges_quad, body, self.boundary, boundary_basis=self.boundary_basis) # initialization of intermediate vars @@ -957,22 +963,19 @@ def _initialize(self): # ... # put the body in tests and trials for loops if is_bilinear: - init_stmts += init_loop_basis( indices_test, ranges_test, self.discrete_boundary ) - init_stmts += init_loop_basis( indices_trial, ranges_trial, self.discrete_boundary ) + init_stmts += init_loop_basis( indices_test, ranges_test, self.boundary ) + init_stmts += init_loop_basis( indices_trial, ranges_trial, self.boundary ) - body = select_loops(indices_test, ranges_test, body, - self.discrete_boundary, + body = select_loops(indices_test, ranges_test, body, self.boundary, boundary_basis=self.boundary_basis) - body = select_loops(indices_trial, ranges_trial, body, - self.discrete_boundary, + body = select_loops(indices_trial, ranges_trial, body, self.boundary, boundary_basis=self.boundary_basis) if is_linear: - init_stmts += init_loop_basis( indices_test, ranges_test, self.discrete_boundary ) + init_stmts += init_loop_basis( indices_test, ranges_test, self.boundary ) - body = select_loops(indices_test, ranges_test, body, - self.discrete_boundary, + body = select_loops(indices_test, ranges_test, body, self.boundary, boundary_basis=self.boundary_basis) # ... @@ -1016,8 +1019,7 @@ def _initialize(self): # evaluation of the area if the mapping is not used if not mapping: stmts = [AugAssign(self.area, '+', weighted_vol)] - stmts = select_loops( indices_quad, ranges_quad, stmts, - self.discrete_boundary, + stmts = select_loops( indices_quad, ranges_quad, stmts, self.boundary, boundary_basis=self.boundary_basis) body = stmts + body @@ -1078,11 +1080,11 @@ def __new__(cls, kernel, name=None, discrete_space=None, comm=None, prefix='assembly', mapping=mapping, is_rational_mapping=is_rational_mapping) - obj._kernel = kernel + obj._kernel = kernel obj._discrete_space = discrete_space - obj._comm = comm - obj._discrete_boundary = kernel.discrete_boundary - obj._backend = backend + obj._comm = comm + obj._boundary = kernel.boundary + obj._backend = backend # update dependencies obj._dependencies += [kernel] @@ -1163,8 +1165,8 @@ def _initialize(self): n_cols = kernel.n_cols axis_bnd = [] - if self.discrete_boundary: - axis_bnd = [i[0] for i in self.discrete_boundary] + if self.boundary: + axis_bnd = [self.boundary.axis] # ... declarations @@ -1410,10 +1412,10 @@ def _initialize(self): weights_in_elm, weights, test_basis_in_elm, test_basis, trial_basis_in_elm, trial_basis, - is_bilinear, self.discrete_boundary ) + is_bilinear, self.boundary ) body = select_loops(indices_elm, ranges_elm, body, - self.kernel.discrete_boundary, boundary_basis=False) + self.kernel.boundary, boundary_basis=False) body = init_stmts + body diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index e6af6f014..289c7c66c 100644 --- a/psydac/api/ast/utilities.py +++ b/psydac/api/ast/utilities.py @@ -16,6 +16,7 @@ from sympde.topology import ScalarField from sympde.topology import VectorField, IndexedVectorField from sympde.topology import Mapping +from sympde.topology import Boundary from sympde.topology.derivatives import _partial_derivatives from sympde.topology.derivatives import _logical_partial_derivatives from sympde.topology.derivatives import get_atom_derivatives @@ -545,15 +546,17 @@ def rationalize_eval_mapping(mapping, nderiv, space, indices_quad): return stmts #============================================================================== -def filter_product(indices, args, discrete_boundary): +def filter_product(indices, args, boundary): mask = [] ext = [] - if discrete_boundary: - # TODO improve using namedtuple or a specific class ? to avoid the 0 index - # => make it easier to understand - mask = [i[0] for i in discrete_boundary] - ext = [i[1] for i in discrete_boundary] + if boundary: + + if isinstance(boundary, Boundary): + mask = [boundary.axis] + ext = [boundary.ext] + else: + raise TypeError # discrete_boundary gives the perpendicular indices, then we need to # remove them from directions @@ -565,15 +568,17 @@ def filter_product(indices, args, discrete_boundary): #============================================================================== # TODO remove it later -def filter_loops(indices, ranges, body, discrete_boundary, boundary_basis=False): +def filter_loops(indices, ranges, body, boundary, boundary_basis=False): quad_mask = [] quad_ext = [] - if discrete_boundary: - # TODO improve using namedtuple or a specific class ? to avoid the 0 index - # => make it easier to understand - quad_mask = [i[0] for i in discrete_boundary] - quad_ext = [i[1] for i in discrete_boundary] + if boundary: + + if isinstance(boundary, Boundary): + quad_mask = [boundary.axis] + quad_ext = [boundary.ext] + else: + raise TypeError # discrete_boundary gives the perpendicular indices, then we need to # remove them from directions @@ -604,15 +609,17 @@ def filter_loops(indices, ranges, body, discrete_boundary, boundary_basis=False) return body #============================================================================== -def select_loops(indices, ranges, body, discrete_boundary, boundary_basis=False): +def select_loops(indices, ranges, body, boundary, boundary_basis=False): quad_mask = [] quad_ext = [] - if discrete_boundary: - # TODO improve using namedtuple or a specific class ? to avoid the 0 index - # => make it easier to understand - quad_mask = [i[0] for i in discrete_boundary] - quad_ext = [i[1] for i in discrete_boundary] + if boundary: + + if isinstance(boundary, Boundary): + quad_mask = [boundary.axis] + quad_ext = [boundary.ext] + else: + raise TypeError # discrete_boundary gives the perpendicular indices, then we need to # remove them from directions @@ -664,13 +671,14 @@ def fusion_loops(loops): return loops_cp #============================================================================== -def compute_normal_vector(vector, discrete_boundary, mapping): +def compute_normal_vector(vector, boundary, mapping): dim = len(vector) - pdim = dim - len(discrete_boundary) - if len(discrete_boundary) > 1: raise NotImplementedError('TODO') + pdim = dim - 1 - face = discrete_boundary[0] - axis = face[0] ; ext = face[1] + if isinstance(boundary, Boundary): + axis = boundary.axis ; ext = boundary.ext + else: + raise TypeError map_stmts = [] body = [] @@ -682,8 +690,8 @@ def compute_normal_vector(vector, discrete_boundary, mapping): else: M = mapping - inv_jac = Symbol('inv_jac') - det_jac = Symbol('det_jac') + inv_jac_bnd = Symbol('inv_jac_bnd') + det_jac_bnd = Symbol('det_jac_bnd') # ... construct jacobian on manifold lines = [] @@ -711,7 +719,7 @@ def compute_normal_vector(vector, discrete_boundary, mapping): # Fortran between sqrt and Pow(, 1/2)? j = (sum(J[i]**2 for i in range(0, dim)))**(1/2) - values = [inv_jac*J[1], -inv_jac*J[0]] + values = [inv_jac_bnd*J[1], -inv_jac_bnd*J[0]] elif dim == 3: @@ -720,14 +728,14 @@ def compute_normal_vector(vector, discrete_boundary, mapping): values = Cross_3d(x_s, x_t) j = (sum(J[i]**2 for i in range(0, dim)))**(1/2) - values = [inv_jac*v for v in values] + values = [inv_jac_bnd*v for v in values] # change the orientation values = [ext*i for i in values] - map_stmts += [Assign(det_jac, j)] - map_stmts += [Assign(inv_jac, 1./j)] + map_stmts += [Assign(det_jac_bnd, j)] + map_stmts += [Assign(inv_jac_bnd, 1./j)] for i in range(0, dim): body += [Assign(vector[i], values[i])] @@ -739,7 +747,7 @@ def compute_normal_vector(vector, discrete_boundary, mapping): return map_stmts, body #============================================================================== -def compute_tangent_vector(vector, discrete_boundary, mapping): +def compute_tangent_vector(vector, boundary, mapping): raise NotImplementedError('TODO') diff --git a/psydac/api/basic.py b/psydac/api/basic.py index 750669df2..5fe956c9d 100644 --- a/psydac/api/basic.py +++ b/psydac/api/basic.py @@ -521,8 +521,8 @@ class BasicDiscrete(BasicCodeGen): def __init__(self, expr, kernel_expr, **kwargs): # ... - target = kwargs.pop('target', None) - boundary = kwargs.pop('boundary', None) + target = kwargs.pop('target', None) + boundary = kwargs.pop('boundary', None) # ... # ... @@ -543,10 +543,6 @@ def __init__(self, expr, kernel_expr, **kwargs): if not( boundary is target ): raise ValueError('> Unconsistent boundary with symbolic model') - boundary = [boundary.axis, boundary.ext] - boundary = [boundary] - boundary_basis = True # TODO set it to False for Nitch method - # boundary is now a list of boundaries # TODO shall we keep it this way? since this is the simplest # interface to be able to compute Functional on a curve in 3d @@ -611,7 +607,7 @@ def _create_ast(self, expr, tag, **kwargs): target = target, mapping = mapping, is_rational_mapping = is_rational_mapping, - discrete_boundary = boundary, + boundary = boundary, boundary_basis = boundary_basis, symbolic_space = symbolic_space, backend = backend ) diff --git a/psydac/api/fem.py b/psydac/api/fem.py index e51ab217c..5816536aa 100644 --- a/psydac/api/fem.py +++ b/psydac/api/fem.py @@ -73,6 +73,12 @@ def __init__(self, expr, kernel_expr, *args, **kwargs): kwargs['is_rational_mapping'] = is_rational_mapping kwargs['comm'] = domain_h.comm + boundary = kwargs.pop('boundary', []) + if boundary and isinstance(boundary, list): + kwargs['boundary'] = boundary[0] + elif boundary: + kwargs['boundary'] = boundary + BasicDiscrete.__init__(self, expr, kernel_expr, **kwargs) # ... @@ -158,6 +164,12 @@ def __init__(self, expr, kernel_expr, *args, **kwargs): kwargs['is_rational_mapping'] = is_rational_mapping kwargs['comm'] = domain_h.comm + boundary = kwargs.pop('boundary', []) + if boundary and isinstance(boundary, list): + kwargs['boundary'] = boundary[0] + elif boundary: + kwargs['boundary'] = boundary + BasicDiscrete.__init__(self, expr, kernel_expr, **kwargs) # ... diff --git a/psydac/api/grid.py b/psydac/api/grid.py index 963f0208b..db41e8d5a 100644 --- a/psydac/api/grid.py +++ b/psydac/api/grid.py @@ -6,6 +6,7 @@ from psydac.core.bsplines import quadrature_grid from psydac.core.bsplines import find_span from psydac.core.bsplines import basis_funs_all_ders +from psydac.core.bsplines import basis_ders_on_quad_grid from psydac.fem.splines import SplineSpace from psydac.fem.tensor import TensorFemSpace from psydac.fem.vector import ProductFemSpace @@ -162,8 +163,8 @@ def __init__( self, V, axis, ext, quad_order=None ): QuadratureGrid.__init__( self, V, quad_order=quad_order ) - points = self.points - weights = self.weights + points = self.points + weights = self.weights # ... if V.ldim == 1: @@ -173,8 +174,8 @@ def __init__( self, V, axis, ext, quad_order=None ): bounds[-1] = V.domain[0] bounds[1] = V.domain[1] - points[axis] = np.asarray([[bounds[ext]]]) - weights[axis] = np.asarray([[1.]]) + points[axis] = np.asarray([[bounds[ext]]]) + weights[axis] = np.asarray([[1.]]) elif V.ldim in [2, 3]: assert( isinstance( V, TensorFemSpace ) ) @@ -183,13 +184,22 @@ def __init__( self, V, axis, ext, quad_order=None ): bounds[-1] = V.spaces[axis].domain[0] bounds[1] = V.spaces[axis].domain[1] - points[axis] = np.asarray([[bounds[ext]]]) - weights[axis] = np.asarray([[1.]]) + points[axis] = np.asarray([[bounds[ext]]]) + weights[axis] = np.asarray([[1.]]) # ... - self._points = points - self._weights = weights + self._axis = axis + self._ext = ext + self._points = points + self._weights = weights + @property + def axis(self): + return self._axis + + @property + def ext(self): + return self._ext #============================================================================== def create_fem_assembly_grid(V, quad_order=None, nderiv=1): @@ -243,7 +253,23 @@ def __init__( self, V, grid, nderiv ): else: self._spans = [g.spans for g in global_quad_grid] self._basis = [g.basis for g in global_quad_grid] - + + # Modify data for boundary grid + if isinstance(grid, BoundaryQuadratureGrid): + axis = grid.axis + ext = grid.ext + if isinstance(V, ProductFemSpace): + for i in range(len(V.spaces)): + sp_space = V.spaces[i].spaces[axis] + points = grid.points[axis] + boundary_basis = basis_ders_on_quad_grid(sp_space.knots, sp_space.degree, points, nderiv) + self._basis[i][axis][0:1, :, :, 0:1] = boundary_basis + else: + sp_space = V.spaces[axis] + points = grid.points[axis] + boundary_basis = basis_ders_on_quad_grid(sp_space.knots, sp_space.degree, points, nderiv) + self._basis[axis][0:1, :, :, 0:1] = boundary_basis + @property def basis(self): return self._basis From 49d5b2dc790e9fad35693ffcd7f1fc8461156bc3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 12 Nov 2019 18:28:16 +0100 Subject: [PATCH 02/23] Fix import statements for system of equations. --- psydac/api/ast/fem.py | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py index 9d8098fa9..1a2207cbb 100644 --- a/psydac/api/ast/fem.py +++ b/psydac/api/ast/fem.py @@ -249,7 +249,6 @@ def __new__(cls, weak_form, kernel_expr, target=None, # ... if boundary: if not isinstance(boundary, Boundary): - print(type(boundary)) raise TypeError('> Expecting a Boundary for boundary') # ... @@ -749,6 +748,7 @@ def _initialize(self): self._vector_fields_logical = vector_fields_logical self._vector_fields_coeffs = vector_fields_coeffs self._mapping_coeffs = mapping_coeffs + self._imports = set() # ... # ranges @@ -1032,11 +1032,27 @@ def _initialize(self): body = len_quads + body # get math functions and constants - math_elements = math_atoms_as_str(self.kernel_expr, 'numpy') + def get_math_elements(function_body, lib): + math_elements = set() + for i in function_body: + if isinstance(i, For): + new = get_math_elements(i.body, lib) + math_elements = math_elements.union(new) + elif isinstance(i, (Assign, AugAssign)): + new = math_atoms_as_str(i.rhs, lib) + math_elements = math_elements.union(new) + elif isinstance(i, FunctionCall): + pass + else: + raise TypeError(i) + return math_elements + + math_elements = get_math_elements(body, 'numpy') math_imports = [Import(e, 'numpy') for e in math_elements] imports += math_imports - self._imports = imports + self._imports = self._imports.union(imports) + # function args mats_args = tuple([mats[i] for i in range(start, end) if not( i in zero_terms )]) func_args = fields_coeffs + vector_fields_coeffs + mapping_coeffs + mats_args @@ -1054,7 +1070,7 @@ def _initialize(self): funcs[i_row][i_col] = FunctionDef(self.name+'_'+str(i_row)+str(i_col), list(func_args), [], body, decorators=decorators,header=header) - + return funcs #============================================================================== From 519625a19c189685748b044fc368b9a21570737a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 12 Nov 2019 18:29:30 +0100 Subject: [PATCH 03/23] Add function 'compute_boundary_jacobian' to api/ast/utilities: - This function computes the correct metric determinant to be used for integration over a boundary manifold; - With this function we can correctly assemble boundary integrals also in the case where the expression does not contain a NormalVector or TangentVector. TODO: fix calculation of metric determinant in function ComputeNormalVector, using the same formula. --- psydac/api/ast/fem.py | 30 ++++++++++++++++-------------- psydac/api/ast/utilities.py | 26 ++++++++++++++++++++++++++ 2 files changed, 42 insertions(+), 14 deletions(-) diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py index 1a2207cbb..93da0e2c1 100644 --- a/psydac/api/ast/fem.py +++ b/psydac/api/ast/fem.py @@ -54,6 +54,7 @@ from .evaluation import EvalQuadratureMapping, EvalQuadratureField, EvalQuadratureVectorField from .utilities import random_string from .utilities import build_pythran_types_header, variables +from .utilities import compute_boundary_jacobian from .utilities import compute_normal_vector, compute_tangent_vector from .utilities import select_loops, filter_product from .utilities import compute_atoms_expr @@ -848,27 +849,28 @@ def _initialize(self): # ... if mapping: - inv_jac = Symbol('inv_jac') - det_jac = Symbol('det_jac') if isinstance(self.target, Boundary): det_jac_bnd = Symbol('det_jac_bnd') + if not vectors: + body += compute_boundary_jacobian(self.target, mapping) - # ... inv jacobian - jac = mapping.det_jacobian - jac = SymbolicExpr(jac) - # ... + else: - body += [Assign(det_jac, jac)] - body += [Assign(inv_jac, 1./jac)] + inv_jac = Symbol('inv_jac') + det_jac = Symbol('det_jac') - # TODO do we use the same inv_jac? - # if not isinstance(self.target, Boundary): - # body += [Assign(inv_jac, 1/jac)] + # ... inv jacobian + jac = mapping.det_jacobian + jac = SymbolicExpr(jac) + # ... - init_map = OrderedDict(sorted(init_map.items())) - for stmt in list(init_map.values()): - body += [stmt.subs(1/jac, inv_jac)] + body += [Assign(det_jac, jac)] + body += [Assign(inv_jac, 1/det_jac)] + + init_map = OrderedDict(sorted(init_map.items())) + for stmt in list(init_map.values()): + body += [stmt.subs(1/jac, inv_jac)] else: body += [Assign(coordinates[i],positions[i][indices_quad[i]]) diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index 289c7c66c..8af118b35 100644 --- a/psydac/api/ast/utilities.py +++ b/psydac/api/ast/utilities.py @@ -670,6 +670,32 @@ def fusion_loops(loops): else: return loops_cp +#============================================================================== +def compute_boundary_jacobian(boundary, mapping): + + # Sanity check on arguments + if isinstance(boundary, Boundary): + axis = boundary.axis + ext = boundary.ext + else: + raise TypeError(boundary) + + if not isinstance(mapping, Mapping): + raise TypeError(mapping) + + # Compute metric determinant g on manifold + J = SymbolicExpr(mapping.jacobian) + Jm = J[:, [i for i in range(J.shape[1]) if i != axis]] + g = (Jm.T * Jm).det() + + # Create statements for computing sqrt(g) and 1/sqrt(g) + inv_jac_bnd = Symbol('inv_jac_bnd') + det_jac_bnd = Symbol('det_jac_bnd') + map_stmts = [Assign(det_jac_bnd, sympy_sqrt(g)), + Assign(inv_jac_bnd, 1/det_jac_bnd)] + + return map_stmts + #============================================================================== def compute_normal_vector(vector, boundary, mapping): dim = len(vector) From 7955baecb74d74821e3bb78a7ff1f162779b5abf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 14 Nov 2019 18:59:39 +0100 Subject: [PATCH 04/23] Add more tests for biharmonic eqn. --- mesh/{quart_circle.h5 => quarter_annulus.h5} | Bin .../api/tests/test_api_2d_scalar_mapping.py | 64 ++++++++++++++++-- .../tests/test_api_glt_2d_scalar_mapping.py | 6 +- 3 files changed, 60 insertions(+), 10 deletions(-) rename mesh/{quart_circle.h5 => quarter_annulus.h5} (100%) diff --git a/mesh/quart_circle.h5 b/mesh/quarter_annulus.h5 similarity index 100% rename from mesh/quart_circle.h5 rename to mesh/quarter_annulus.h5 diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py index 2706f7959..f1cf0b950 100644 --- a/psydac/api/tests/test_api_2d_scalar_mapping.py +++ b/psydac/api/tests/test_api_2d_scalar_mapping.py @@ -1,4 +1,15 @@ # -*- coding: UTF-8 -*- +# +# A note on the mappings used in these tests: +# +# - 'identity_2d.h5' is the identity mapping on the unit square [0, 1] X [0, 1] +# +# - 'collela_2d.h5' is a NURBS mapping from the unit square [0, 1]^2 to the +# larger square [-1, 1]^2, with deformations going as sin(pi x) * sin(pi y) +# +# - 'quarter_annulus.h5' is a NURBS transformation from the unit square [0, 1]^2 +# to the quarter annulus in the lower-left quadrant of the Cartesian place +# (hence both x and y are negative), with r_min = 0.5 and r_max = 0.5 from mpi4py import MPI from sympy import pi, cos, sin @@ -360,8 +371,8 @@ def test_api_poisson_2d_dir_collela(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dir_quart_circle(): - filename = os.path.join(mesh_dir, 'quart_circle.h5') +def test_api_poisson_2d_dir_quarter_annulus(): + filename = os.path.join(mesh_dir, 'quarter_annulus.h5') from sympy.abc import x,y @@ -597,7 +608,7 @@ def test_api_laplace_2d_neu_identity(): #============================================================================== def test_api_biharmonic_2d_dir_identity(): - filename = os.path.join(mesh_dir, 'quart_circle.h5') + filename = os.path.join(mesh_dir, 'identity_2d.h5') from sympy.abc import x,y @@ -606,15 +617,54 @@ def test_api_biharmonic_2d_dir_identity(): l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f) - expected_l2_error = 0.3092792236008727 - expected_h1_error = 1.3320441589030887 - expected_h2_error = 6.826223014197719 + expected_l2_error = 0.015086415626060034 + expected_h1_error = 0.08773346232941553 + expected_h2_error = 1.9368842415954024 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) assert( abs(h2_error - expected_h2_error) < 1.e-7) -# TODO: add biharmonic on Collela and quart_circle mappings +#============================================================================== +def test_api_biharmonic_2d_dir_collela(): + filename = os.path.join(mesh_dir, 'collela_2d.h5') + + from sympy.abc import x,y + + solution = (cos(pi*x/2)*cos(pi*y/2))**2 + f = laplace(laplace(solution)) + + l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f) + + expected_l2_error = 0.10977627980052021 + expected_h1_error = 0.32254511059711766 + expected_h2_error = 1.87205519824758 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + assert( abs(h2_error - expected_h2_error) < 1.e-7) + +#============================================================================== +def test_api_biharmonic_2d_dir_quarter_annulus(): + filename = os.path.join(mesh_dir, 'quarter_annulus.h5') + + from sympy.abc import x,y + + r_in = 0.5 + r_out = 1 + kappa = 1 / 0.00643911127175763 + solution = kappa * (x * y * (x**2 + y**2 - r_in**2) * (x**2 + y**2 - r_out**2))**2 + f = laplace(laplace(solution)) + + l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f) + + expected_l2_error = 0.016730298635551484 + expected_h1_error = 0.21243295522291714 + expected_h2_error = 7.572921831391894 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + assert( abs(h2_error - expected_h2_error) < 1.e-7) ##============================================================================== ## TODO DEBUG, not working since merge with devel diff --git a/psydac/api/tests/test_api_glt_2d_scalar_mapping.py b/psydac/api/tests/test_api_glt_2d_scalar_mapping.py index f66c6a2b3..eaf418d3f 100644 --- a/psydac/api/tests/test_api_glt_2d_scalar_mapping.py +++ b/psydac/api/tests/test_api_glt_2d_scalar_mapping.py @@ -96,8 +96,8 @@ def test_api_glt_poisson_2d_dir_collela(): assert(np.allclose([error], [0.04655602895206486])) #============================================================================== -def test_api_glt_poisson_2d_dir_quart_circle(): - filename = os.path.join(mesh_dir, 'quart_circle.h5') +def test_api_glt_poisson_2d_dir_quarter_annulus(): + filename = os.path.join(mesh_dir, 'quarter_annulus.h5') error = run_poisson_2d_dir(filename) assert(np.allclose([error], [0.04139096668630673])) @@ -116,4 +116,4 @@ def teardown_function(): #test_api_glt_poisson_2d_dir_identity() #test_api_glt_poisson_2d_dir_collela() -#test_api_glt_poisson_2d_dir_quart_circle() +#test_api_glt_poisson_2d_dir_quarter_annulus() From 7beb6ce55fb003803ba63867f20cfe7d2e9c4327 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 14 Nov 2019 19:48:28 +0100 Subject: [PATCH 05/23] Fix broken unit tests w/ Neumann BCs. --- .../api/tests/test_api_2d_scalar_mapping.py | 62 +++++++++---------- 1 file changed, 30 insertions(+), 32 deletions(-) diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py index f1cf0b950..1ac8ab539 100644 --- a/psydac/api/tests/test_api_2d_scalar_mapping.py +++ b/psydac/api/tests/test_api_2d_scalar_mapping.py @@ -520,22 +520,21 @@ def test_api_poisson_2d_dirneu_identity_123(): ##============================================================================== -# TODO DEBUG, not working since merge with devel -#def test_api_poisson_2d_dirneu_collela_1(): -# filename = os.path.join(mesh_dir, 'collela_2d.h5') -# -# from sympy.abc import x,y -# -# solution = cos(0.25*pi*x)*sin(pi*y) -# f = (17./16.)*pi**2*solution -# -# l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, [{'axis': 0, 'ext': -1}]) -# -# expected_l2_error = 0.013540717397796734 -# expected_h1_error = 0.19789463571596025 -# -# assert( abs(l2_error - expected_l2_error) < 1.e-7) -# assert( abs(h1_error - expected_h1_error) < 1.e-7) +def test_api_poisson_2d_dirneu_collela_1(): + filename = os.path.join(mesh_dir, 'collela_2d.h5') + + from sympy.abc import x,y + + solution = sin(0.25*pi*(x-1))*sin(pi*y) + f = (17./16.)*pi**2*solution + + l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, [{'axis': 0, 'ext': -1}]) + + expected_l2_error = 0.013540717397796734 + expected_h1_error = 0.19789463571596025 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== def test_api_poisson_2d_dirneu_collela_2(): @@ -555,22 +554,21 @@ def test_api_poisson_2d_dirneu_collela_2(): assert( abs(h1_error - expected_h1_error) < 1.e-7) ##============================================================================== -## TODO DEBUG, not working since merge with devel -#def test_api_poisson_2d_dirneu_collela_3(): -# filename = os.path.join(mesh_dir, 'collela_2d.h5') -# -# from sympy.abc import x,y -# -# solution = cos(0.25*pi*y)*sin(pi*x) -# f = (17./16.)*pi**2*solution -# -# l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, [{'axis': 1, 'ext': -1}]) -# -# expected_l2_error = 0.013540717397817427 -# expected_h1_error = 0.19789463571595994 -# -# assert( abs(l2_error - expected_l2_error) < 1.e-7) -# assert( abs(h1_error - expected_h1_error) < 1.e-7) +def test_api_poisson_2d_dirneu_collela_3(): + filename = os.path.join(mesh_dir, 'collela_2d.h5') + + from sympy.abc import x,y + + solution = sin(0.25*pi*(y-1))*sin(pi*x) + f = (17./16.)*pi**2*solution + + l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, [{'axis': 1, 'ext': -1}]) + + expected_l2_error = 0.013540717397817427 + expected_h1_error = 0.19789463571595994 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== def test_api_poisson_2d_dirneu_collela_4(): From 3c1dcdc79dd2eb6ed62ceac0b21fe1912a44292e Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 14 Nov 2019 20:11:31 +0100 Subject: [PATCH 06/23] Add unit tests with non-zero Neumann BCs. --- psydac/api/tests/test_api_2d_scalar.py | 67 ++++++++++++++++++++++++++ 1 file changed, 67 insertions(+) diff --git a/psydac/api/tests/test_api_2d_scalar.py b/psydac/api/tests/test_api_2d_scalar.py index 01ea0e0bb..eff7c8b19 100644 --- a/psydac/api/tests/test_api_2d_scalar.py +++ b/psydac/api/tests/test_api_2d_scalar.py @@ -512,6 +512,73 @@ def test_api_poisson_2d_dirneu_123(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) +#============================================================================== +def test_api_poisson_2d_dir_zero_neu_nonzero_1(): + + from sympy.abc import x,y + + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution + + l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 0, 'ext': -1}], + ncells=[2**3, 2**3], degree=[2, 2]) + + expected_l2_error = 0.00021786960672322118 + expected_h1_error = 0.01302350067761091 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#============================================================================== +def test_api_poisson_2d_dir_zero_neu_nonzero_2(): + + from sympy.abc import x,y + + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution + + l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 0, 'ext': 1}], + ncells=[2**3, 2**3], degree=[2, 2]) + + expected_l2_error = 0.00021786960672322118 + expected_h1_error = 0.01302350067761091 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#============================================================================== +def test_api_poisson_2d_dir_zero_neu_nonzero_3(): + + from sympy.abc import x,y + + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution + + l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 1, 'ext': -1}], + ncells=[2**3, 2**3], degree=[2, 2]) + + expected_l2_error = 0.00021786960672322118 + expected_h1_error = 0.01302350067761091 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#============================================================================== +def test_api_poisson_2d_dir_zero_neu_nonzero_4(): + + from sympy.abc import x,y + + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution + + l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 1, 'ext': 1}], + ncells=[2**3, 2**3], degree=[2, 2]) + + expected_l2_error = 0.00021786960672322118 + expected_h1_error = 0.01302350067761091 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== def test_api_laplace_2d_neu(): From ab60adb5e6e6515de3df5fab5a372b5d9ddf8f00 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 14 Nov 2019 20:26:01 +0100 Subject: [PATCH 07/23] Add unit tests w/ non-zero Neumann BCs & mapping: So far only 2D Poisson equation on square domain (identity mapping). Tests only pass if inhomogeneous Neumann BCs are along direction 0, and fail otherwise. Identified bug in 'compute_normal_vector'. --- .../api/tests/test_api_2d_scalar_mapping.py | 72 +++++++++++++++++++ 1 file changed, 72 insertions(+) diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py index 1ac8ab539..de5d23ca0 100644 --- a/psydac/api/tests/test_api_2d_scalar_mapping.py +++ b/psydac/api/tests/test_api_2d_scalar_mapping.py @@ -587,6 +587,78 @@ def test_api_poisson_2d_dirneu_collela_4(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) +#============================================================================== +def test_api_poisson_2d_dir_zero_neu_nonzero_identity_1(): + filename = os.path.join(mesh_dir, 'identity_2d.h5') + + from sympy.abc import x,y + + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution + + l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, + [{'axis': 0, 'ext': -1}]) + + expected_l2_error = 0.00021786960671761908 + expected_h1_error = 0.01302350067761177 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#============================================================================== +def test_api_poisson_2d_dir_zero_neu_nonzero_identity_2(): + filename = os.path.join(mesh_dir, 'identity_2d.h5') + + from sympy.abc import x,y + + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution + + l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, + [{'axis': 0, 'ext': 1}]) + + expected_l2_error = 0.00021786960671761908 + expected_h1_error = 0.01302350067761177 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#============================================================================== +def test_api_poisson_2d_dir_zero_neu_nonzero_identity_3(): + filename = os.path.join(mesh_dir, 'identity_2d.h5') + + from sympy.abc import x,y + + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution + + l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, + [{'axis': 1, 'ext': -1}]) + + expected_l2_error = 0.00021786960671761908 + expected_h1_error = 0.01302350067761177 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#============================================================================== +def test_api_poisson_2d_dir_zero_neu_nonzero_identity_4(): + filename = os.path.join(mesh_dir, 'identity_2d.h5') + + from sympy.abc import x,y + + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution + + l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, + [{'axis': 1, 'ext': 1}]) + + expected_l2_error = 0.00021786960671761908 + expected_h1_error = 0.01302350067761177 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + #============================================================================== def test_api_laplace_2d_neu_identity(): filename = os.path.join(mesh_dir, 'identity_2d.h5') From 1f5008d3283303a88b8e6e5f1fc1713bb1bcadee Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 14 Nov 2019 20:35:31 +0100 Subject: [PATCH 08/23] Fix bug in function 'compute_normal_vector' of api/ast/utilities: The normal vector is identified with a dual basis vector of the coordinate transformation (reversed if needed). Hence it is computed as a row of the inverse of the Jacobian matrix of the mapping. This procedure is simple and general, and it applies to domains with any number of dimensions. TODO: verify if a normalization is required. --- psydac/api/ast/fem.py | 3 +- psydac/api/ast/utilities.py | 85 +++++++++---------------------------- 2 files changed, 22 insertions(+), 66 deletions(-) diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py index 93da0e2c1..772dedded 100644 --- a/psydac/api/ast/fem.py +++ b/psydac/api/ast/fem.py @@ -852,8 +852,7 @@ def _initialize(self): if isinstance(self.target, Boundary): det_jac_bnd = Symbol('det_jac_bnd') - if not vectors: - body += compute_boundary_jacobian(self.target, mapping) + body += compute_boundary_jacobian(self.target, mapping) else: diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index 8af118b35..78e61a43b 100644 --- a/psydac/api/ast/utilities.py +++ b/psydac/api/ast/utilities.py @@ -698,77 +698,34 @@ def compute_boundary_jacobian(boundary, mapping): #============================================================================== def compute_normal_vector(vector, boundary, mapping): - dim = len(vector) - pdim = dim - 1 + # Sanity check on arguments if isinstance(boundary, Boundary): - axis = boundary.axis ; ext = boundary.ext + axis = boundary.axis + ext = boundary.ext else: - raise TypeError - - map_stmts = [] - body = [] + raise TypeError(boundary) - if not mapping: - values = np.zeros(dim) - values[axis] = ext + if mapping is None: + map_stmts = [] + values = [(ext if i==axis else 0) for i, v in enumerate(vector)] else: - M = mapping - inv_jac_bnd = Symbol('inv_jac_bnd') - det_jac_bnd = Symbol('det_jac_bnd') - - # ... construct jacobian on manifold - lines = [] - n_row,n_col = M.jacobian.shape - range_row = [i for i in range(0,n_row) if not(i == axis)] - range_col = range(0,n_col) - for i_row in range_row: - line = [] - for i_col in range_col: - line.append(M.jacobian[i_col, i_row]) - - lines.append(line) - - J = Matrix(lines) - # ... - - J = SymbolicExpr(J) - - if dim == 1: - raise NotImplementedError('TODO') - - elif dim == 2: - J = J[0,:] - # TODO shall we use sympy_sqrt here? is there any difference in - # Fortran between sqrt and Pow(, 1/2)? - j = (sum(J[i]**2 for i in range(0, dim)))**(1/2) - - values = [inv_jac_bnd*J[1], -inv_jac_bnd*J[0]] - - elif dim == 3: - - x_s = J[0,:] - x_t = J[1,:] - - values = Cross_3d(x_s, x_t) - j = (sum(J[i]**2 for i in range(0, dim)))**(1/2) - values = [inv_jac_bnd*v for v in values] - - - # change the orientation - values = [ext*i for i in values] - - map_stmts += [Assign(det_jac_bnd, j)] - map_stmts += [Assign(inv_jac_bnd, 1./j)] - - for i in range(0, dim): - body += [Assign(vector[i], values[i])] - -# print(map_stmts) -# print(body) -# import sys; sys.exit(0) + # Given the Jacobian matrix J, we need to extract the (i=axis) row + # of J^(-1). For efficiency we separately compute det(J) and the + # cofactors C[i=0...dim] of the (j=axis) column of J. + # NOTE: we also change the orientation according to 'ext' + J = SymbolicExpr(mapping.jacobian) + det = J.det() + cof = [J.cofactor(i, j=axis) for i in range(J.shape[0])] + + inv_jac = Symbol('inv_jac') + map_stmts = [Assign(inv_jac, 1 / det)] + values = [ext * inv_jac * cof[i] for i in range(J.shape[0])] + + # Create statements for computing normal vector components + body = [Assign(lhs, rhs) for lhs, rhs in zip(vector, values)] return map_stmts, body From e90947c8f88295190c56a60049d02c14b6b1d30d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Fri, 15 Nov 2019 19:04:45 +0100 Subject: [PATCH 09/23] Clean up Kernel class in api/ast/fem, and add extensive comments. --- psydac/api/ast/fem.py | 306 ++++++++++++++++++------------------ psydac/api/ast/utilities.py | 42 +++-- 2 files changed, 171 insertions(+), 177 deletions(-) diff --git a/psydac/api/ast/fem.py b/psydac/api/ast/fem.py index 772dedded..6752cb3c5 100644 --- a/psydac/api/ast/fem.py +++ b/psydac/api/ast/fem.py @@ -209,7 +209,15 @@ def area_eval_mapping(mapping, area, dim, indices_quad, weight): #============================================================================== # target is used when there are multiple expression (domain/boundaries) class Kernel(SplBasic): + """ + Generate the AST of a function for computing an integral form over a + single domain element, or boundary element. + For a bilinear form, such a function will compute an 'element matrix'. + For a linear form, it will compute an 'element vector'. + For a functional, it will compute an 'element value' (scalar). + + """ def __new__(cls, weak_form, kernel_expr, target=None, boundary=None, name=None, boundary_basis=None, mapping=None, is_rational_mapping=None,symbolic_space=None, backend=None): @@ -635,35 +643,20 @@ def _initialize(self): rank = 1 if isinstance(expr, Matrix): - sh = expr.shape - - # ... - mats = [] - for i_row in range(0, sh[0]): - for i_col in range(0, sh[1]): - mats.append('mat_{}{}'.format(i_row, i_col)) + nr, nc = expr.shape + mats = ['mat_{}{}'.format(i, j) for j in range(nc) for i in range(nr)] + vals = ['val_{}{}'.format(i, j) for j in range(nc) for i in range(nr)] mats = variables(mats, dtype='real', rank=rank, cls=IndexedVariable) - # ... - - # ... - v = [] - for i_row in range(0, sh[0]): - for i_col in range(0, sh[1]): - v.append('v_{}{}'.format(i_row, i_col)) - - v = variables(v, 'real') - # ... + vals = variables(vals, 'real') expr = expr[:] ln = len(expr) else: mats = (IndexedVariable('mat_00', dtype='real', rank=rank),) - - v = (Variable('real', 'v_00'),) - ln = 1 - + vals = (Variable('real', 'val_00'),) expr = [expr] + ln = 1 # ... looking for 0 terms zero_terms = [i for i,e in enumerate(expr) if e == 0] @@ -710,6 +703,16 @@ def _initialize(self): positions = variables('quad_u1:%s'%(dim+1), dtype='real', rank=1, cls=IndexedVariable) + # Used only if there is a mapping + inv_jac = Symbol('inv_jac') + det_jac = Symbol('det_jac') + + # Used only in the case of a boundary assembly + vectors = self.kernel_expr.atoms(BoundaryVector) + normal_vec = symbols('normal_1:%d'%(dim+1)) + tangent_vec = symbols('tangent_1:%d'%(dim+1)) + det_jac_bnd = symbols('det_jac_bnd') + # ... # ... @@ -759,10 +762,6 @@ def _initialize(self): # ... # body of kernel - - init_basis = OrderedDict() - init_map = OrderedDict() - init_stmts, map_stmts = compute_atoms_expr(atomic_expr, indices_quad, indices_test, @@ -774,176 +773,157 @@ def _initialize(self): is_linear, mapping) - for stmt in init_stmts: - init_basis[str(stmt.lhs)] = stmt + # Sort statements according to left-hand-side name + init_basis = sorted(init_stmts, key=lambda s: str(s.lhs)) + init_map = sorted( map_stmts, key=lambda s: str(s.lhs)) - for stmt in map_stmts: - init_map[str(stmt.lhs)] = stmt - + # Prepare matrix of functions for assembling each element of a block + # matrix or block vector. if unique_scalar_space: - ln = 1 + ln = 1 funcs = [[None]] - else: - funcs = [[None]*self._n_cols for i in range(self._n_rows)] + funcs = [[None] * self._n_cols] * self._n_rows + + # Compute jacobian determinant expression + if mapping: + jac = SymbolicExpr(mapping.det_jacobian) + #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + # FOR LOOP OVER THE EQUATIONS IN A SYSTEM + #++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ for indx in range(ln): if not unique_scalar_space and indx in zero_terms: continue - + elif not unique_scalar_space: start = indx end = indx + 1 - i_row = indx//self._n_cols - i_col = indx -i_row*self._n_cols - + i_row = indx // self._n_cols + i_col = indx - i_row * self._n_cols + else: i_row = 0 i_col = 0 start = 0 end = len(expr) - + + #------------------------------------------------------------------ + # Body of inner quadrature loop: all quantities are scalars + #------------------------------------------------------------------ + + # New body body = [] - init_basis = OrderedDict(sorted(init_basis.items())) - body += list(init_basis.values()) + # Evaluation of basis functions (and their derivatives) + # along each direction in the parametric domain + body += init_basis + + # Evaluation of physical coordinates (x, y, z). + # In the case of a mapping, also compute the components + # of the Jacobian matrix, as well as its determinant. if mapping: - body += [Assign(lhs, rhs[indices_quad]) for lhs, rhs in zip(mapping_elements, - mapping_values)] + body += [Assign(lhs, rhs[indices_quad]) + for lhs, rhs in zip(mapping_elements, mapping_values)] - # ... normal/tangent vectors - init_map_bnd = OrderedDict() + body += [Assign(det_jac, jac), + Assign(inv_jac, 1/det_jac)] + else: + body += [Assign(coordinates[i], positions[i][indices_quad[i]]) for i in range(dim)] + + # If assembling a boundary integral, compute any normal/tangent vectors, + # as well as the metric determinant for integration on the manifold. if isinstance(self.target, Boundary): - vectors = self.kernel_expr.atoms(BoundaryVector) - normal_vec = symbols('normal_1:%d'%(dim+1)) - tangent_vec = symbols('tangent_1:%d'%(dim+1)) for vector in vectors: + if isinstance(vector, NormalVector): # replace n[i] by its scalar components - for i in range(0, dim): - expr = [e.subs(vector[i], normal_vec[i]) for e in expr] - - map_stmts, stmts = compute_normal_vector(normal_vec, - self.boundary, - mapping) + expr = [e.subs(zip(vector, normal_vec)) for e in expr] + body += compute_normal_vector(locals(), normal_vec, self.target, mapping) elif isinstance(vector, TangentVector): # replace t[i] by its scalar components - for i in range(0, dim): - expr = [e.subs(vector[i], tangent_vec[i]) for e in expr] + expr = [e.subs(zip(vector, tangent_vec)) for e in expr] + body += compute_tangent_vector(locals(), tangent_vec, self.target, mapping) - map_stmts, stmts = compute_tangent_vector(tangent_vec, - self.boundary, - mapping) - - for stmt in map_stmts: - init_map_bnd[str(stmt.lhs)] = stmt + else: + raise TypeError(vector) - init_map_bnd = OrderedDict(sorted(init_map_bnd.items())) - for stmt in list(init_map_bnd.values()): - body += [stmt] + body += compute_boundary_jacobian(locals(), self.target, mapping) - body += stmts - # ... + # Evaluation of physical derivatives of basis functions. If there + # is no mapping, 'init_map' is an empty list and nothing is done. + body += [stmt.subs(1/jac, inv_jac) for stmt in init_map] + # Evaluation of any scalar and vector fields (and their derivatives) + # that appear in the kernel expression. if mapping: + body += [Assign(f, v[indices_quad]) for f, v in zip(fields_logical, fields_val)] + body += [Assign(f, v[indices_quad]) + for f, v in zip(vector_fields_logical, vector_fields_val)] - if isinstance(self.target, Boundary): - det_jac_bnd = Symbol('det_jac_bnd') - body += compute_boundary_jacobian(self.target, mapping) - - else: - - inv_jac = Symbol('inv_jac') - det_jac = Symbol('det_jac') - - # ... inv jacobian - jac = mapping.det_jacobian - jac = SymbolicExpr(jac) - # ... - - body += [Assign(det_jac, jac)] - body += [Assign(inv_jac, 1/det_jac)] - - init_map = OrderedDict(sorted(init_map.items())) - for stmt in list(init_map.values()): - body += [stmt.subs(1/jac, inv_jac)] + # Convert logical derivatives to physical derivatives + body += [stmt.subs(1/jac, inv_jac) for stmt in self._map_stmts_fields.values()] else: - body += [Assign(coordinates[i],positions[i][indices_quad[i]]) - for i in range(dim)] - # ... - - # ... + body += [Assign(f, v[indices_quad]) for f, v in zip(fields, fields_val)] + body += [Assign(f, v[indices_quad]) + for f, v in zip(vector_fields, vector_fields_val)] + + # Compute the scalar coefficient V_{ij..} by which the point value + # of the kernel expression E must be multiplied in order to + # approximate the integral by a high-order quadrature rule: + # + # integral(domain, E) ~= sum_{ij..} (E_{ij..} * V_{ij..} * J_{ij..}), + # + # where E_{ij..} = E(x1_i, x2_j, ...), J_{ij..} is the Jacobian + # determinant, and V_{ij..} = w1_i * w2_j * ... is simply the + # product of the rescaled 1D Gaussian weights along each direction. + # + # On a boundary with xk=const, we must not multiply by the weight + # along xk; hence the need for a 'filter product'. weighted_vol = filter_product(indices_quad, weighted_vols, self.boundary) - # ... - - # ... - # add fields and vector fields - if not mapping: - # ... fields - for i in range(len(fields_val)): - body.append(Assign(fields[i],fields_val[i][indices_quad])) - # ... - - # ... vector_fields - for i in range(len(vector_fields_val)): - body.append(Assign(vector_fields[i],vector_fields_val[i][indices_quad])) - # ... - - else: - # ... fields - for i in range(len(fields_val)): - body.append(Assign(fields_logical[i],fields_val[i][indices_quad])) - # ... - - # ... vector_fields - # if vector_fields_val: - # print(vector_fields_logical) - # print(vector_fields_val) - # import sys; sys.exit(0) - for i in range(len(vector_fields_val)): - body.append(Assign(vector_fields_logical[i],vector_fields_val[i][indices_quad])) - # ... - - # ... substitute expression of inv_jac - for k,stmt in self._map_stmts_fields.items(): - body += [stmt.subs(1/jac, inv_jac)] - # ... + # Multiply by the correct metric determinant # TODO use positive mapping all the time? Abs? if mapping: if isinstance(self.target, Boundary): - weighted_vol = weighted_vol * Abs(det_jac_bnd) + weighted_vol *= Abs(det_jac_bnd) else: - weighted_vol = weighted_vol * Abs(det_jac) - - body.append(Assign(wvol,weighted_vol)) + weighted_vol *= Abs(det_jac) + # Update value of integrals: + # 1. Assign 'weighted_vol' expression to variable 'wvol', + # 2. Multiply each kernel expression by 'wvol', + # 3. Increment each integral variables by the expressions above. + body.append(Assign(wvol, weighted_vol)) for i in range(start, end): if not( i in zero_terms ): - e = SymbolicExpr(Mul(expr[i],wvol)) + e = SymbolicExpr(Mul(expr[i], wvol)) + body.append(AugAssign(vals[i], '+', e)) - body.append(AugAssign(v[i],'+', e)) - # ... + #------------------------------------------------------------------ + # Body of external loop over basis functions + #------------------------------------------------------------------ - # ... stmts for initializtion: only when boundary is present - init_stmts = [] - # ... + # Collect initialization statements in this list. These statements + # will be placed at the beginning of the function, outside any loops. + init_stmts = [] - # ... - # put the body in for loops of quadrature points + # If we are assemblying on a boundary, some indices are set to zero init_stmts += init_loop_quadrature( indices_quad, ranges_quad, self.boundary ) + # Put the body inside for loops of quadrature points body = select_loops( indices_quad, ranges_quad, body, self.boundary, boundary_basis=self.boundary_basis) - # initialization of intermediate vars - init_vars = [Assign(v[i],0.0) for i in range(start, end) if not( i in zero_terms )] + # Initialize intermediate variables 'vals': they are set to zero + # just before the inner quadrature loop, and they are incremented + # inside that loop. + init_vars = [Assign(vals[i], 0.0) for i in range(start, end) if not( i in zero_terms )] body = init_vars + body - # ... if dim_trial: trial_idxs = tuple([indices_trial[i]+trial_pads[i]-indices_test[i] for i in range(dim)]) @@ -954,15 +934,18 @@ def _initialize(self): if is_bilinear or is_linear: for i in range(start, end): if not( i in zero_terms ): - body.append(Assign(mats[i][idxs],v[i])) + body.append(Assign(mats[i][idxs], vals[i])) elif is_function: for i in range(start, end): if not( i in zero_terms ): - body.append(Assign(mats[i][0],v[i])) + body.append(Assign(mats[i][0], vals[i])) - # ... - # put the body in tests and trials for loops + #------------------------------------------------------------------ + # Body of kernel + #------------------------------------------------------------------ + + # Put the body inside for loops of test and trial functions if is_bilinear: init_stmts += init_loop_basis( indices_test, ranges_test, self.boundary ) init_stmts += init_loop_basis( indices_trial, ranges_trial, self.boundary ) @@ -1001,8 +984,6 @@ def _initialize(self): body = [FunctionCall(eval_field.func, args)] + body - imports = [] - # call eval vector_field for eval_vector_field in self.eval_vector_fields: args = test_degrees + basis_test + vector_fields_coeffs + vector_fields_val @@ -1032,16 +1013,20 @@ def _initialize(self): len_quads = [Assign(k, Len(u)) for k,u in zip(qds_dim, positions)] body = len_quads + body - # get math functions and constants + #------------------------------------------------------------------ + # Compute module-wise import statements + #------------------------------------------------------------------ + + # Search recursively for math functions and constants def get_math_elements(function_body, lib): math_elements = set() for i in function_body: if isinstance(i, For): new = get_math_elements(i.body, lib) - math_elements = math_elements.union(new) + math_elements.update(new) elif isinstance(i, (Assign, AugAssign)): new = math_atoms_as_str(i.rhs, lib) - math_elements = math_elements.union(new) + math_elements.update(new) elif isinstance(i, FunctionCall): pass else: @@ -1050,16 +1035,20 @@ def get_math_elements(function_body, lib): math_elements = get_math_elements(body, 'numpy') math_imports = [Import(e, 'numpy') for e in math_elements] + self._imports.update(math_imports) - imports += math_imports - self._imports = self._imports.union(imports) + #------------------------------------------------------------------ + # Identify function arguments + #------------------------------------------------------------------ - # function args mats_args = tuple([mats[i] for i in range(start, end) if not( i in zero_terms )]) func_args = fields_coeffs + vector_fields_coeffs + mapping_coeffs + mats_args - func_args = self.build_arguments(func_args) + #------------------------------------------------------------------ + # Add header and decorators, if needed + #------------------------------------------------------------------ + decorators = {} header = None if self.backend['name'] == 'pyccel': @@ -1070,13 +1059,22 @@ def get_math_elements(function_body, lib): header = build_pythran_types_header(self.name, func_args) funcs[i_row][i_col] = FunctionDef(self.name+'_'+str(i_row)+str(i_col), list(func_args), [], body, - decorators=decorators,header=header) + decorators=decorators, header=header) return funcs #============================================================================== class Assembly(SplBasic): + """ + Generate the AST of a function for computing an integral form over the + whole domain, or the whole boundary. This is obtained by 'assembling' the + single contributions over each element (see class Kernel). + + For a bilinear form, such a function will assemble a matrix. + For a linear form, it will assemble a vector. + For a functional, it will assemble a scalar value. + """ def __new__(cls, kernel, name=None, discrete_space=None, comm=None, mapping=None, is_rational_mapping=None, backend=None): diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index 78e61a43b..461f8ee7b 100644 --- a/psydac/api/ast/utilities.py +++ b/psydac/api/ast/utilities.py @@ -671,7 +671,7 @@ def fusion_loops(loops): return loops_cp #============================================================================== -def compute_boundary_jacobian(boundary, mapping): +def compute_boundary_jacobian(parent_namespace, boundary, mapping=None): # Sanity check on arguments if isinstance(boundary, Boundary): @@ -680,24 +680,23 @@ def compute_boundary_jacobian(boundary, mapping): else: raise TypeError(boundary) - if not isinstance(mapping, Mapping): - raise TypeError(mapping) + if mapping is None: + stmts = [] - # Compute metric determinant g on manifold - J = SymbolicExpr(mapping.jacobian) - Jm = J[:, [i for i in range(J.shape[1]) if i != axis]] - g = (Jm.T * Jm).det() + else: + # Compute metric determinant g on manifold + J = SymbolicExpr(mapping.jacobian) + Jm = J[:, [i for i in range(J.shape[1]) if i != axis]] + g = (Jm.T * Jm).det() - # Create statements for computing sqrt(g) and 1/sqrt(g) - inv_jac_bnd = Symbol('inv_jac_bnd') - det_jac_bnd = Symbol('det_jac_bnd') - map_stmts = [Assign(det_jac_bnd, sympy_sqrt(g)), - Assign(inv_jac_bnd, 1/det_jac_bnd)] + # Create statements for computing sqrt(g) + det_jac_bnd = parent_namespace['det_jac_bnd'] + stmts = [Assign(det_jac_bnd, sympy_sqrt(g))] - return map_stmts + return stmts #============================================================================== -def compute_normal_vector(vector, boundary, mapping): +def compute_normal_vector(parent_namespace, vector, boundary, mapping=None): # Sanity check on arguments if isinstance(boundary, Boundary): @@ -708,8 +707,7 @@ def compute_normal_vector(vector, boundary, mapping): if mapping is None: - map_stmts = [] - values = [(ext if i==axis else 0) for i, v in enumerate(vector)] + values = [(ext if i==axis else 0) for i, v in enumerate(vector)] else: # Given the Jacobian matrix J, we need to extract the (i=axis) row @@ -720,20 +718,18 @@ def compute_normal_vector(vector, boundary, mapping): det = J.det() cof = [J.cofactor(i, j=axis) for i in range(J.shape[0])] - inv_jac = Symbol('inv_jac') - map_stmts = [Assign(inv_jac, 1 / det)] - values = [ext * inv_jac * cof[i] for i in range(J.shape[0])] + inv_jac = parent_namespace['inv_jac'] + values = [ext * inv_jac * cof[i] for i in range(J.shape[0])] # Create statements for computing normal vector components - body = [Assign(lhs, rhs) for lhs, rhs in zip(vector, values)] + stmts = [Assign(lhs, rhs) for lhs, rhs in zip(vector, values)] - return map_stmts, body + return stmts #============================================================================== -def compute_tangent_vector(vector, boundary, mapping): +def compute_tangent_vector(parent_namespace, vector, boundary, mapping): raise NotImplementedError('TODO') - #============================================================================== _range = re.compile('([0-9]*:[0-9]+|[a-zA-Z]?:[a-zA-Z])') From 71d4590656db6edac6025dcd44b1466129562db4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 20 Nov 2019 21:18:08 +0100 Subject: [PATCH 10/23] Fix order of (trial, test) function in DiscreteEquation. --- psydac/api/discretization.py | 9 +++++---- psydac/api/fem.py | 4 ++-- 2 files changed, 7 insertions(+), 6 deletions(-) diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index fc110c5a8..f15e505c5 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -101,9 +101,10 @@ def __init__(self, expr, *args, **kwargs): # since lhs and rhs are calls, we need to take their expr # ... - test_trial = args[1] - test_space = test_trial[0] - trial_space = test_trial[1] + domain = args[0] + trial_test = args[1] + trial_space = trial_test[0] + test_space = trial_test[1] # ... # ... @@ -121,7 +122,7 @@ def __init__(self, expr, *args, **kwargs): kwargs['boundary'] = boundaries_lhs newargs = list(args) - newargs[1] = test_trial + newargs[1] = trial_test self._lhs = discretize(expr.lhs, *newargs, **kwargs) # ... diff --git a/psydac/api/fem.py b/psydac/api/fem.py index 5816536aa..0fb384c0e 100644 --- a/psydac/api/fem.py +++ b/psydac/api/fem.py @@ -82,8 +82,8 @@ def __init__(self, expr, kernel_expr, *args, **kwargs): BasicDiscrete.__init__(self, expr, kernel_expr, **kwargs) # ... - test_space = self.spaces[0] - trial_space = self.spaces[1] + trial_space = self.spaces[0] + test_space = self.spaces[1] # ... # ... From ff97dc2446459757b770db419c4a95f80b376aa8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 26 Nov 2019 19:07:37 +0100 Subject: [PATCH 11/23] New driver for 2D Poisson tests w/out mapping. A unique function is now used to define unit tests for the 2D Poisson equation in a unit square. It handles any combination of homogeneous and inhomogeneous boundary conditions of Dirichlet or Neumann type. --- psydac/api/tests/test_api_2d_scalar.py | 264 +++++++++++-------------- 1 file changed, 120 insertions(+), 144 deletions(-) diff --git a/psydac/api/tests/test_api_2d_scalar.py b/psydac/api/tests/test_api_2d_scalar.py index eff7c8b19..46d56f94a 100644 --- a/psydac/api/tests/test_api_2d_scalar.py +++ b/psydac/api/tests/test_api_2d_scalar.py @@ -20,115 +20,63 @@ from psydac.api.discretization import discretize #============================================================================== -def run_poisson_2d_dir(solution, f, ncells, degree, comm=None): +def get_boundaries(*args): - # ... abstract model - domain = Square() - - V = ScalarFunctionSpace('V', domain) - - F = element_of(V, name='F') - - v = element_of(V, name='v') - u = element_of(V, name='u') - - int_0 = lambda expr: integral(domain , expr) - - expr = dot(grad(v), grad(u)) - a = BilinearForm((v,u), int_0(expr)) - - expr = f*v - l = LinearForm(v, int_0(expr)) - - error = F - solution - l2norm = Norm(error, domain, kind='l2') - h1norm = Norm(error, domain, kind='h1') - - bc = EssentialBC(u, 0, domain.boundary) - equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc) - # ... - - # ... create the computational domain from a topological domain - domain_h = discretize(domain, ncells=ncells, comm=comm) - # ... - - # ... discrete spaces - Vh = discretize(V, domain_h, degree=degree) - # ... - - # ... dsicretize the equation using Dirichlet bc - equation_h = discretize(equation, domain_h, [Vh, Vh]) - # ... - - # ... discretize norms - l2norm_h = discretize(l2norm, domain_h, Vh) - h1norm_h = discretize(h1norm, domain_h, Vh) - # ... - - # ... solve the discrete equation - x = equation_h.solve() - # ... - - # ... - phi = FemField( Vh, x ) - # ... + if not args: + return () + else: + assert all(1 <= a <= 4 for a in args) + assert len(set(args)) == len(args) - # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) - # ... + boundaries = {1: {'axis': 0, 'ext': -1}, + 2: {'axis': 0, 'ext': 1}, + 3: {'axis': 1, 'ext': -1}, + 4: {'axis': 1, 'ext': 1}} - return l2_error, h1_error + return tuple(boundaries[i] for i in args) #============================================================================== -def run_poisson_2d_dirneu(solution, f, boundary, ncells, degree, comm=None): +def run_poisson_2d(solution, f, dir_zero_boundary, dir_nonzero_boundary, ncells, degree, comm=None): - assert( isinstance(boundary, (list, tuple)) ) + assert( isinstance(dir_zero_boundary , (list, tuple)) ) + assert( isinstance(dir_nonzero_boundary, (list, tuple)) ) # ... abstract model domain = Square() - V = ScalarFunctionSpace('V', domain) - - B_neumann = [domain.get_boundary(**kw) for kw in boundary] - if len(B_neumann) == 1: - B_neumann = B_neumann[0] - - else: - B_neumann = Union(*B_neumann) - - x,y = domain.coordinates - - F = element_of(V, name='F') - - v = element_of(V, name='v') - u = element_of(V, name='u') + B_dirichlet_0 = Union(*[domain.get_boundary(**kw) for kw in dir_zero_boundary]) + B_dirichlet_i = Union(*[domain.get_boundary(**kw) for kw in dir_nonzero_boundary]) + B_dirichlet = Union(B_dirichlet_0, B_dirichlet_i) + B_neumann = domain.boundary.complement(B_dirichlet) + V = ScalarFunctionSpace('V', domain) + u = element_of(V, name='u') + v = element_of(V, name='v') nn = NormalVector('nn') - int_0 = lambda expr: integral(domain , expr) - int_1 = lambda expr: integral(B_neumann , expr) - - expr = dot(grad(v), grad(u)) - a = BilinearForm((v,u), int_0(expr)) + # Bilinear Form: a(u, v) + a = BilinearForm((u, v), integral(domain, dot(grad(u), grad(v)))) - expr = f*v - l0 = LinearForm(v, int_0(expr)) + # Linear form: l(v) + l0 = LinearForm(v, integral(domain, f * v)) + if B_neumann: + l1 = LinearForm(v, integral(B_neumann, v * dot(grad(solution), nn))) + l = LinearForm(v, l0(v) + l1(v)) + else: + l = l0 - expr = v*dot(grad(solution), nn) - l_B_neumann = LinearForm(v, int_1(expr)) + # Dirichlet boundary conditions + bc = [] + if B_dirichlet_0: bc += [EssentialBC(u, 0, B_dirichlet_0)] + if B_dirichlet_i: bc += [EssentialBC(u, solution, B_dirichlet_i)] - expr = l0(v) + l_B_neumann(v) - l = LinearForm(v, expr) + # Variational model + equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc) - error = F-solution + # Error norms + error = u - solution l2norm = Norm(error, domain, kind='l2') h1norm = Norm(error, domain, kind='h1') - - B_dirichlet = domain.boundary.complement(B_neumann) - bc = EssentialBC(u, 0, B_dirichlet) - - equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc) # ... # ... create the computational domain from a topological domain @@ -153,12 +101,12 @@ def run_poisson_2d_dirneu(solution, f, boundary, ncells, degree, comm=None): # ... # ... - phi = FemField( Vh, x ) + uh = FemField( Vh, x ) # ... # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) + l2_error = l2norm_h.assemble(u=uh) + h1_error = h1norm_h.assemble(u=uh) # ... return l2_error, h1_error @@ -299,7 +247,6 @@ def run_biharmonic_2d_dir(solution, f, ncells, degree, comm=None): return l2_error, h1_error - #============================================================================== def run_poisson_user_function_2d_dir(f, solution, ncells, degree, comm=None): @@ -369,16 +316,18 @@ def run_poisson_user_function_2d_dir(f, solution, ncells, degree, comm=None): # SERIAL TESTS ############################################################################### -#============================================================================== -def test_api_poisson_2d_dir_1(): +def test_api_poisson_2d_dir0_1234(): from sympy.abc import x,y solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*sin(pi*x)*sin(pi*y) - l2_error, h1_error = run_poisson_2d_dir(solution, f, - ncells=[2**3,2**3], degree=[2,2]) + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.00021808678604760232 expected_h1_error = 0.013023570720360362 @@ -387,15 +336,18 @@ def test_api_poisson_2d_dir_1(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dirneu_1(): +def test_api_poisson_2d_dir0_234_neu0_1(): from sympy.abc import x,y solution = cos(0.5*pi*x)*sin(pi*y) f = (5./4.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 0, 'ext': -1}], - ncells=[2**3,2**3], degree=[2,2]) + dir_zero_boundary = get_boundaries(2, 3, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.00015546057796452772 expected_h1_error = 0.00926930278452745 @@ -404,15 +356,18 @@ def test_api_poisson_2d_dirneu_1(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dirneu_2(): +def test_api_poisson_2d_dir0_134_neu0_2(): from sympy.abc import x,y solution = sin(0.5*pi*x)*sin(pi*y) f = (5./4.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 0, 'ext': 1}], - ncells=[2**3,2**3], degree=[2,2]) + dir_zero_boundary = get_boundaries(1, 3, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.0001554605779481901 expected_h1_error = 0.009269302784527256 @@ -421,15 +376,18 @@ def test_api_poisson_2d_dirneu_2(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dirneu_3(): +def test_api_poisson_2d_dir0_124_neu0_3(): from sympy.abc import x,y solution = sin(pi*x)*cos(0.5*pi*y) f = (5./4.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 1, 'ext': -1}], - ncells=[2**3,2**3], degree=[2,2]) + dir_zero_boundary = get_boundaries(1, 2, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.0001554605779681901 expected_h1_error = 0.009269302784528678 @@ -438,15 +396,18 @@ def test_api_poisson_2d_dirneu_3(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dirneu_4(): +def test_api_poisson_2d_dir0_123_neu0_4(): from sympy.abc import x,y solution = sin(pi*x)*sin(0.5*pi*y) f = (5./4.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 1, 'ext': 1}], - ncells=[2**3,2**3], degree=[2,2]) + dir_zero_boundary = get_boundaries(1, 2, 3) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.00015546057796339546 expected_h1_error = 0.009269302784526841 @@ -455,17 +416,18 @@ def test_api_poisson_2d_dirneu_4(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dirneu_13(): +def test_api_poisson_2d_dir0_24_neu0_13(): from sympy.abc import x,y solution = cos(0.5*pi*x)*cos(0.5*pi*y) f = (1./2.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, - [{'axis': 0, 'ext': -1}, - {'axis': 1, 'ext': -1}], - ncells=[2**3,2**3], degree=[2,2]) + dir_zero_boundary = get_boundaries(2, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 2.6119892736036942e-05 expected_h1_error = 0.0016032430287934746 @@ -474,17 +436,18 @@ def test_api_poisson_2d_dirneu_13(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dirneu_24(): +def test_api_poisson_2d_dir0_13_neu0_24(): from sympy.abc import x,y solution = sin(0.5*pi*x)*sin(0.5*pi*y) f = (1./2.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, - [{'axis': 0, 'ext': 1}, - {'axis': 1, 'ext': 1}], - ncells=[2**3,2**3], degree=[2,2]) + dir_zero_boundary = get_boundaries(1, 3) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 2.611989253883369e-05 expected_h1_error = 0.0016032430287973409 @@ -493,18 +456,18 @@ def test_api_poisson_2d_dirneu_24(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dirneu_123(): +def test_api_poisson_2d_dir0_4_neu0_123(): from sympy.abc import x,y solution = cos(pi*x)*cos(0.5*pi*y) f = 5./4.*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, - [{'axis': 0, 'ext': -1}, - {'axis': 0, 'ext': 1}, - {'axis': 1, 'ext': -1}], - ncells=[2**3,2**3], degree=[2,2]) + dir_zero_boundary = get_boundaries(4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.00015494478505412876 expected_h1_error = 0.009242166414700994 @@ -513,15 +476,18 @@ def test_api_poisson_2d_dirneu_123(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dir_zero_neu_nonzero_1(): +def test_api_poisson_2d_dir0_234_neui_1(): from sympy.abc import x,y solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 0, 'ext': -1}], - ncells=[2**3, 2**3], degree=[2, 2]) + dir_zero_boundary = get_boundaries(2, 3, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.00021786960672322118 expected_h1_error = 0.01302350067761091 @@ -530,15 +496,18 @@ def test_api_poisson_2d_dir_zero_neu_nonzero_1(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dir_zero_neu_nonzero_2(): +def test_api_poisson_2d_dir0_134_neui_2(): from sympy.abc import x,y solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 0, 'ext': 1}], - ncells=[2**3, 2**3], degree=[2, 2]) + dir_zero_boundary = get_boundaries(1, 3, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.00021786960672322118 expected_h1_error = 0.01302350067761091 @@ -547,15 +516,18 @@ def test_api_poisson_2d_dir_zero_neu_nonzero_2(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dir_zero_neu_nonzero_3(): +def test_api_poisson_2d_dir0_124_neui_3(): from sympy.abc import x,y solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 1, 'ext': -1}], - ncells=[2**3, 2**3], degree=[2, 2]) + dir_zero_boundary = get_boundaries(1, 2, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.00021786960672322118 expected_h1_error = 0.01302350067761091 @@ -564,15 +536,18 @@ def test_api_poisson_2d_dir_zero_neu_nonzero_3(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dir_zero_neu_nonzero_4(): +def test_api_poisson_2d_dir0_123_neui_4(): from sympy.abc import x,y solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(solution, f, [{'axis': 1, 'ext': 1}], - ncells=[2**3, 2**3], degree=[2, 2]) + dir_zero_boundary = get_boundaries(1, 2, 3) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.00021786960672322118 expected_h1_error = 0.01302350067761091 @@ -617,8 +592,6 @@ def test_api_biharmonic_2d_dir_1(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) - - #============================================================================== def test_api_poisson_user_function_2d_dir_1(): @@ -651,16 +624,19 @@ def f(x,y): #============================================================================== @pytest.mark.parallel -def test_api_poisson_2d_dir_1_parallel(): +def test_api_poisson_2d_dir0_1234_parallel(): from sympy.abc import x,y solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*sin(pi*x)*sin(pi*y) - l2_error, h1_error = run_poisson_2d_dir(solution, f, - ncells=[2**3,2**3], degree=[2,2], - comm=MPI.COMM_WORLD) + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2], + comm=MPI.COMM_WORLD) expected_l2_error = 0.00021808678604760232 expected_h1_error = 0.013023570720360362 From 8bf3a2937d3123da362450d4ace22128339abb67 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 20 Nov 2019 21:19:30 +0100 Subject: [PATCH 12/23] Implement inhomogeneous Dirichlet BCs. --- psydac/api/discretization.py | 68 ++++++++++++++++++++++++-- psydac/api/tests/test_api_2d_scalar.py | 41 ++++++++++++++++ 2 files changed, 105 insertions(+), 4 deletions(-) diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index f15e505c5..239c38fe0 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -139,8 +139,9 @@ def __init__(self, expr, *args, **kwargs): self._bc = bc self._linear_system = None + self._domain = domain self._trial_space = trial_space - self._test_space = test_space + self._test_space = test_space @property def expr(self): @@ -155,13 +156,17 @@ def rhs(self): return self._rhs @property - def test_space(self): - return self._test_space + def domain(self): + return self._domain @property def trial_space(self): return self._trial_space + @property + def test_space(self): + return self._test_space + @property def bc(self): return self._bc @@ -195,7 +200,7 @@ def assemble(self, **kwargs): self._linear_system = LinearSystem(M, rhs) def solve(self, **kwargs): - settings = kwargs.pop('settings', _default_solver) + settings = kwargs.pop('settings', _default_solver.copy()) rhs = kwargs.pop('rhs', None) if rhs: kwargs['assemble_rhs'] = False @@ -207,6 +212,61 @@ def solve(self, **kwargs): L = LinearSystem(L.lhs, rhs) self._linear_system = L + #---------------------------------------------------------------------- + # [YG, 18/11/2019] + # + # Impose inhomogeneous Dirichlet boundary conditions through + # L2 projection on the boundary. This requires setting up a + # new variational formulation and solving the resulting linear + # system to obtain a solution that does not live in the space + # of homogeneous solutions. Such a solution is then used as + # initial guess when the model equation is to be solved by an + # iterative method. Our current method of solution does not + # modify the initial guess at the boundary. + # + if self.bc: + + # Inhomogeneous Dirichlet boundary conditions + idbcs = [i for i in self.bc if i.rhs != 0] + + if idbcs: + + from sympde.expr import integral + from sympde.expr import find + from sympde.topology import element_of, ScalarTestFunction + + # Extract trial functions from model equation + u = self.expr.trial_functions + + # Create test functions in same space of trial functions + # TODO: check if we should generate random names + V = ProductSpace(*[ui.space for ui in u]) + v = element_of(V, name='v:{}'.format(len(u))) + + # In a system, each essential boundary condition is applied to + # only one component (bc.variable) of the state vector. Hence + # we will select the correct test function using a dictionary. + test_dict = dict(zip(u, v)) + + # Construct variational formulation that performs L2 projection + # of boundary conditions onto the correct space + product = lambda f, g: (f * g if isinstance(g, ScalarTestFunction) else dot(f, g)) + factor = lambda bc : bc.lhs.xreplace(test_dict) + lhs_expr = sum(integral(i.boundary, product(i.lhs, factor(i))) for i in idbcs) + rhs_expr = sum(integral(i.boundary, product(i.rhs, factor(i))) for i in idbcs) + equation = find(u, forall=v, lhs=lhs_expr, rhs=rhs_expr) + + # Discretize weak form and find inhomogeneous solution + domain_h = self.domain + Vh = self.trial_space + equation_h = discretize(equation, domain_h, [Vh, Vh]) + X = equation_h.solve() + + # Use inhomogeneous solution as initial guess to solver + settings['x0'] = X + + #---------------------------------------------------------------------- + return driver_solve(self.linear_system, **settings) #============================================================================== diff --git a/psydac/api/tests/test_api_2d_scalar.py b/psydac/api/tests/test_api_2d_scalar.py index 46d56f94a..452a87595 100644 --- a/psydac/api/tests/test_api_2d_scalar.py +++ b/psydac/api/tests/test_api_2d_scalar.py @@ -556,6 +556,47 @@ def test_api_poisson_2d_dir0_123_neui_4(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== +def test_api_poisson_2d_dir0_123_diri_4(): + + from sympy.abc import x,y + + solution = sin(pi * x) * sin(0.5*pi * y) + f = 5/4*pi**2 * solution + + dir_zero_boundary = get_boundaries(1, 2, 3) + dir_nonzero_boundary = get_boundaries(4) + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) + + expected_l2_error = 0.00015292215711784052 + expected_h1_error = 0.009293161646614652 + + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + +#============================================================================== +def test_api_poisson_2d_dir0_13_diri_24(): + + from sympy.abc import x,y + + solution = sin(3*pi/2 * x) * sin(3*pi/2 * y) + f = 9/2*pi**2 * solution + + dir_zero_boundary = get_boundaries(1, 3) + dir_nonzero_boundary = get_boundaries(2, 4) + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) + + expected_l2_error = 0.0007786454571731944 + expected_h1_error = 0.0449669071240554 + + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + +#============================================================================== + def test_api_laplace_2d_neu(): from sympy.abc import x,y From 565cc38141a44ed6c5061bb6ce2e980d85098291 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 27 Nov 2019 19:26:43 +0100 Subject: [PATCH 13/23] Cleanup 2D Poisson test-case w/out mapping. --- psydac/api/tests/test_api_2d_scalar.py | 40 +++++++++++++------------- 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/psydac/api/tests/test_api_2d_scalar.py b/psydac/api/tests/test_api_2d_scalar.py index 452a87595..19d407785 100644 --- a/psydac/api/tests/test_api_2d_scalar.py +++ b/psydac/api/tests/test_api_2d_scalar.py @@ -38,10 +38,12 @@ def get_boundaries(*args): #============================================================================== def run_poisson_2d(solution, f, dir_zero_boundary, dir_nonzero_boundary, ncells, degree, comm=None): - assert( isinstance(dir_zero_boundary , (list, tuple)) ) - assert( isinstance(dir_nonzero_boundary, (list, tuple)) ) + assert isinstance(dir_zero_boundary , (list, tuple)) + assert isinstance(dir_nonzero_boundary, (list, tuple)) - # ... abstract model + #+++++++++++++++++++++++++++++++ + # 1. Abstract model + #+++++++++++++++++++++++++++++++ domain = Square() B_dirichlet_0 = Union(*[domain.get_boundary(**kw) for kw in dir_zero_boundary]) @@ -54,7 +56,7 @@ def run_poisson_2d(solution, f, dir_zero_boundary, dir_nonzero_boundary, ncells, v = element_of(V, name='v') nn = NormalVector('nn') - # Bilinear Form: a(u, v) + # Bilinear form: a(u, v) a = BilinearForm((u, v), integral(domain, dot(grad(u), grad(v)))) # Linear form: l(v) @@ -77,37 +79,35 @@ def run_poisson_2d(solution, f, dir_zero_boundary, dir_nonzero_boundary, ncells, error = u - solution l2norm = Norm(error, domain, kind='l2') h1norm = Norm(error, domain, kind='h1') - # ... - # ... create the computational domain from a topological domain + #+++++++++++++++++++++++++++++++ + # 2. Discretization + #+++++++++++++++++++++++++++++++ + + # Create computational domain from topological domain domain_h = discretize(domain, ncells=ncells, comm=comm) - # ... - # ... discrete spaces + # Discrete spaces Vh = discretize(V, domain_h, degree=degree) - # ... - # ... dsicretize the equation using Dirichlet bc + # Discretize equation using Dirichlet bc equation_h = discretize(equation, domain_h, [Vh, Vh]) - # ... - # ... discretize norms + # Discretize error norms l2norm_h = discretize(l2norm, domain_h, Vh) h1norm_h = discretize(h1norm, domain_h, Vh) - # ... - # ... solve the discrete equation - x = equation_h.solve() - # ... + #+++++++++++++++++++++++++++++++ + # 3. Solution + #+++++++++++++++++++++++++++++++ - # ... + # Solve linear system + x = equation_h.solve() uh = FemField( Vh, x ) - # ... - # ... compute norms + # Compute error norms l2_error = l2norm_h.assemble(u=uh) h1_error = h1norm_h.assemble(u=uh) - # ... return l2_error, h1_error From d13fd9131f908d4a1e31b9fd34ab4b0355dabcfe Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Wed, 27 Nov 2019 19:34:18 +0100 Subject: [PATCH 14/23] Extend 2D biharmonic test-cases w/out mapping. --- psydac/api/discretization.py | 5 +- psydac/api/tests/test_api_2d_scalar.py | 150 +++++++++++++++++-------- 2 files changed, 108 insertions(+), 47 deletions(-) diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index 239c38fe0..ec6d81298 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -248,9 +248,12 @@ def solve(self, **kwargs): # we will select the correct test function using a dictionary. test_dict = dict(zip(u, v)) + # TODO: use dot product for vector quantities +# product = lambda f, g: (f * g if isinstance(g, ScalarTestFunction) else dot(f, g)) + product = lambda f, g: f * g + # Construct variational formulation that performs L2 projection # of boundary conditions onto the correct space - product = lambda f, g: (f * g if isinstance(g, ScalarTestFunction) else dot(f, g)) factor = lambda bc : bc.lhs.xreplace(test_dict) lhs_expr = sum(integral(i.boundary, product(i.lhs, factor(i))) for i in idbcs) rhs_expr = sum(integral(i.boundary, product(i.rhs, factor(i))) for i in idbcs) diff --git a/psydac/api/tests/test_api_2d_scalar.py b/psydac/api/tests/test_api_2d_scalar.py index 19d407785..c2db441d4 100644 --- a/psydac/api/tests/test_api_2d_scalar.py +++ b/psydac/api/tests/test_api_2d_scalar.py @@ -185,67 +185,81 @@ def run_laplace_2d_neu(solution, f, ncells, degree, comm=None): return l2_error, h1_error #============================================================================== -def run_biharmonic_2d_dir(solution, f, ncells, degree, comm=None): +def run_biharmonic_2d_dir(solution, f, dir_zero_boundary, ncells, degree, comm=None): - # ... abstract model + assert isinstance(dir_zero_boundary, (list, tuple)) + + #+++++++++++++++++++++++++++++++ + # 1. Abstract model + #+++++++++++++++++++++++++++++++ domain = Square() - V = ScalarFunctionSpace('V', domain) + B_dirichlet_0 = Union(*[domain.get_boundary(**kw) for kw in dir_zero_boundary]) + B_dirichlet_i = domain.boundary.complement(B_dirichlet_0) - F = element_of(V, name='F') + V = ScalarFunctionSpace('V', domain) + v = element_of(V, name='v') + u = element_of(V, name='u') + nn = NormalVector('nn') - v = element_of(V, name='v') - u = element_of(V, name='u') + # Bilinear form a: V x V --> R + a = BilinearForm((u, v), integral(domain, laplace(u) * laplace(v))) - int_0 = lambda expr: integral(domain , expr) + # Linear form l: V --> R + l = LinearForm(v, integral(domain, f * v)) - expr = laplace(v) * laplace(u) - a = BilinearForm((v,u),int_0(expr)) - expr = f*v - l = LinearForm(v, int_0(expr)) + # Essential boundary conditions + dn = lambda a: dot(grad(a), nn) + bc = [] + if B_dirichlet_0: + bc += [EssentialBC( u , 0, B_dirichlet_0)] + bc += [EssentialBC(dn(u), 0, B_dirichlet_0)] + if B_dirichlet_i: + bc += [EssentialBC( u , solution , B_dirichlet_i)] + bc += [EssentialBC(dn(u), dn(solution), B_dirichlet_i)] - error = F - solution + # Variational model + equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc) + + # Error norms + error = u - solution l2norm = Norm(error, domain, kind='l2') h1norm = Norm(error, domain, kind='h1') + h2norm = Norm(error, domain, kind='h2') - nn = NormalVector('nn') - bc = [EssentialBC(u, 0, domain.boundary)] - bc += [EssentialBC(dot(grad(u), nn), 0, domain.boundary)] - equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc) - # ... + #+++++++++++++++++++++++++++++++ + # 2. Discretization + #+++++++++++++++++++++++++++++++ - # ... create the computational domain from a topological domain + # Create computational domain from topological domain domain_h = discretize(domain, ncells=ncells, comm=comm) - # ... - # ... discrete spaces + # Discrete spaces Vh = discretize(V, domain_h, degree=degree) - # ... - # ... dsicretize the equation using Dirichlet bc + # Discretize equation using Dirichlet bc equation_h = discretize(equation, domain_h, [Vh, Vh]) - # ... - # ... discretize norms + # Discretize error norms l2norm_h = discretize(l2norm, domain_h, Vh) h1norm_h = discretize(h1norm, domain_h, Vh) - # ... + h2norm_h = discretize(h2norm, domain_h, Vh) - # ... solve the discrete equation - x = equation_h.solve() - # ... + #+++++++++++++++++++++++++++++++ + # 3. Solution + #+++++++++++++++++++++++++++++++ - # ... - phi = FemField( Vh, x ) - # ... + # Solve linear system + x = equation_h.solve() + uh = FemField( Vh, x ) - # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) - # ... + # Compute error norms + l2_error = l2norm_h.assemble(u=uh) + h1_error = h1norm_h.assemble(u=uh) + h2_error = h2norm_h.assemble(u=uh) - return l2_error, h1_error + return l2_error, h1_error, h2_error #============================================================================== def run_poisson_user_function_2d_dir(f, solution, ncells, degree, comm=None): @@ -613,25 +627,69 @@ def test_api_laplace_2d_neu(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_biharmonic_2d_dir_1(): +def test_api_biharmonic_2d_dir0_1234(): from sympy.abc import x,y - from sympde.expr import TerminalExpr - solution = (sin(pi*x)*sin(pi*y))**2 + solution = sin(pi * x)**2 * sin(pi * y)**2 + f = laplace(laplace(solution)) - # compute the analytical solution - f = laplace(laplace(solution)) - f = TerminalExpr(f, dim=2) + dir_zero_boundary = get_boundaries(1, 2, 3, 4) - l2_error, h1_error = run_biharmonic_2d_dir(solution, f, - ncells=[2**3,2**3], degree=[2,2]) + l2_error, h1_error, h2_error = run_biharmonic_2d_dir(solution, f, + dir_zero_boundary, ncells=[2**3, 2**3], degree=[3, 3]) - expected_l2_error = 0.015086415626061608 - expected_h1_error = 0.08773346232942228 + expected_l2_error = 0.00019981371108040476 + expected_h1_error = 0.0063205179028178295 + expected_h2_error = 0.2123929568623994 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) + assert( abs(h2_error - expected_h2_error) < 1.e-7) + +#============================================================================== +@pytest.mark.xfail +def test_api_biharmonic_2d_dir0_123_diri_4(): + + from sympy.abc import x,y + + solution = sin(pi * x)**2 * sin(0.5*pi * y)**2 + f = laplace(laplace(solution)) + + dir_zero_boundary = get_boundaries(1, 2, 3) + + l2_error, h1_error, h2_error = run_biharmonic_2d_dir(solution, f, + dir_zero_boundary, ncells=[2**3, 2**3], degree=[3, 3]) + + print() + print(l2_error) + print(h1_error) + print(h2_error) + print() + + assert False + +#============================================================================== +@pytest.mark.xfail +def test_api_biharmonic_2d_dir0_13_diri_24(): + + from sympy.abc import x,y + + solution = sin(3*pi/2 * x)**2 * sin(3*pi/2 * y)**2 + f = laplace(laplace(solution)) + + dir_zero_boundary = get_boundaries(1, 3) + + l2_error, h1_error, h2_error = run_biharmonic_2d_dir(solution, f, + dir_zero_boundary, ncells=[2**3, 2**3], degree=[3, 3]) + + print() + print(l2_error) + print(h1_error) + print(h2_error) + print() + + assert False #============================================================================== def test_api_poisson_user_function_2d_dir_1(): From c41db0ad4f0f3775c6c012b5f7e2009c701a4d47 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 28 Nov 2019 15:22:05 +0100 Subject: [PATCH 15/23] Reorganize 2D scalar tests w/out mapping --- psydac/api/tests/test_api_2d_scalar.py | 340 +++++++++---------------- 1 file changed, 127 insertions(+), 213 deletions(-) diff --git a/psydac/api/tests/test_api_2d_scalar.py b/psydac/api/tests/test_api_2d_scalar.py index c2db441d4..6a9839ca7 100644 --- a/psydac/api/tests/test_api_2d_scalar.py +++ b/psydac/api/tests/test_api_2d_scalar.py @@ -2,6 +2,7 @@ from mpi4py import MPI from sympy import pi, cos, sin +from sympy.abc import x, y from sympy.utilities.lambdify import implemented_function import pytest @@ -36,7 +37,8 @@ def get_boundaries(*args): return tuple(boundaries[i] for i in args) #============================================================================== -def run_poisson_2d(solution, f, dir_zero_boundary, dir_nonzero_boundary, ncells, degree, comm=None): +def run_poisson_2d(solution, f, dir_zero_boundary, dir_nonzero_boundary, + ncells, degree, comm=None): assert isinstance(dir_zero_boundary , (list, tuple)) assert isinstance(dir_nonzero_boundary, (list, tuple)) @@ -56,10 +58,10 @@ def run_poisson_2d(solution, f, dir_zero_boundary, dir_nonzero_boundary, ncells, v = element_of(V, name='v') nn = NormalVector('nn') - # Bilinear form: a(u, v) + # Bilinear form a: V x V --> R a = BilinearForm((u, v), integral(domain, dot(grad(u), grad(v)))) - # Linear form: l(v) + # Linear form l: V --> R l0 = LinearForm(v, integral(domain, f * v)) if B_neumann: l1 = LinearForm(v, integral(B_neumann, v * dot(grad(solution), nn))) @@ -112,75 +114,79 @@ def run_poisson_2d(solution, f, dir_zero_boundary, dir_nonzero_boundary, ncells, return l2_error, h1_error #============================================================================== -def run_laplace_2d_neu(solution, f, ncells, degree, comm=None): - - # ... abstract model - domain = Square() - - V = ScalarFunctionSpace('V', domain) - - B_neumann = domain.boundary +def run_laplace_2d(solution, f, dir_zero_boundary, dir_nonzero_boundary, + ncells, degree, comm=None): - x,y = domain.coordinates + assert isinstance(dir_zero_boundary , (list, tuple)) + assert isinstance(dir_nonzero_boundary, (list, tuple)) - F = element_of(V, name='F') + #+++++++++++++++++++++++++++++++ + # 1. Abstract model + #+++++++++++++++++++++++++++++++ + domain = Square() - v = element_of(V, name='v') - u = element_of(V, name='u') + B_dirichlet_0 = Union(*[domain.get_boundary(**kw) for kw in dir_zero_boundary]) + B_dirichlet_i = Union(*[domain.get_boundary(**kw) for kw in dir_nonzero_boundary]) + B_dirichlet = Union(B_dirichlet_0, B_dirichlet_i) + B_neumann = domain.boundary.complement(B_dirichlet) + V = ScalarFunctionSpace('V', domain) + u = element_of(V, name='u') + v = element_of(V, name='v') nn = NormalVector('nn') - int_0 = lambda expr: integral(domain , expr) - int_1 = lambda expr: integral(B_neumann , expr) - - expr = dot(grad(v), grad(u)) + v*u - a = BilinearForm((v,u), int_0(expr)) + # Bilinear form a: V x V --> R + a = BilinearForm((u, v), integral(domain, dot(grad(u), grad(v)) + u * v)) - expr = f*v - l0 = LinearForm(v, int_0(expr)) + # Linear form l: V --> R + l0 = LinearForm(v, integral(domain, f * v)) + if B_neumann: + l1 = LinearForm(v, integral(B_neumann, v * dot(grad(solution), nn))) + l = LinearForm(v, l0(v) + l1(v)) + else: + l = l0 - expr = v*dot(grad(solution), nn) - l_B_neumann = LinearForm(v, int_1(expr)) + # Dirichlet boundary conditions + bc = [] + if B_dirichlet_0: bc += [EssentialBC(u, 0, B_dirichlet_0)] + if B_dirichlet_i: bc += [EssentialBC(u, solution, B_dirichlet_i)] - expr = l0(v) + l_B_neumann(v) - l = LinearForm(v, expr) + # Variational model + equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc) - error = F-solution + # Error norms + error = u - solution l2norm = Norm(error, domain, kind='l2') h1norm = Norm(error, domain, kind='h1') - equation = find(u, forall=v, lhs=a(u,v), rhs=l(v)) - # ... + #+++++++++++++++++++++++++++++++ + # 2. Discretization + #+++++++++++++++++++++++++++++++ - # ... create the computational domain from a topological domain + # Create computational domain from topological domain domain_h = discretize(domain, ncells=ncells, comm=comm) - # ... - # ... discrete spaces + # Discrete spaces Vh = discretize(V, domain_h, degree=degree) - # ... - # ... dsicretize the equation using Dirichlet bc + # Discretize equation using Dirichlet bc equation_h = discretize(equation, domain_h, [Vh, Vh]) - # ... - # ... discretize norms + # Discretize error norms l2norm_h = discretize(l2norm, domain_h, Vh) h1norm_h = discretize(h1norm, domain_h, Vh) - # ... - # ... solve the discrete equation - x = equation_h.solve() - # ... + #+++++++++++++++++++++++++++++++ + # 3. Solution + #+++++++++++++++++++++++++++++++ - # ... - phi = FemField( Vh, x ) - # ... + # Solve linear system + x = equation_h.solve() + uh = FemField( Vh, x ) - # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) - # ... + # Compute error norms + l2_error = l2norm_h.assemble(u=uh) + h1_error = h1norm_h.assemble(u=uh) return l2_error, h1_error @@ -198,8 +204,8 @@ def run_biharmonic_2d_dir(solution, f, dir_zero_boundary, ncells, degree, comm=N B_dirichlet_i = domain.boundary.complement(B_dirichlet_0) V = ScalarFunctionSpace('V', domain) - v = element_of(V, name='v') u = element_of(V, name='u') + v = element_of(V, name='v') nn = NormalVector('nn') # Bilinear form a: V x V --> R @@ -208,7 +214,6 @@ def run_biharmonic_2d_dir(solution, f, dir_zero_boundary, ncells, degree, comm=N # Linear form l: V --> R l = LinearForm(v, integral(domain, f * v)) - # Essential boundary conditions dn = lambda a: dot(grad(a), nn) bc = [] @@ -261,78 +266,14 @@ def run_biharmonic_2d_dir(solution, f, dir_zero_boundary, ncells, degree, comm=N return l2_error, h1_error, h2_error -#============================================================================== -def run_poisson_user_function_2d_dir(f, solution, ncells, degree, comm=None): - - # ... abstract model - domain = Square() - x,y = domain.coordinates - - f = implemented_function('f', f) - - V = ScalarFunctionSpace('V', domain) - - F = element_of(V, name='F') - - v = element_of(V, name='v') - u = element_of(V, name='u') - - int_0 = lambda expr: integral(domain , expr) - - expr = dot(grad(v), grad(u)) - a = BilinearForm((v,u), int_0(expr)) - - expr = f(x,y)*v - l = LinearForm(v, int_0(expr)) - - error = F - solution - l2norm = Norm(error, domain, kind='l2') - h1norm = Norm(error, domain, kind='h1') - - bc = EssentialBC(u, 0, domain.boundary) - equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc) - # ... - - # ... create the computational domain from a topological domain - domain_h = discretize(domain, ncells=ncells, comm=comm) - # ... - - # ... discrete spaces - Vh = discretize(V, domain_h, degree=degree) - # ... - - # ... dsicretize the equation using Dirichlet bc - equation_h = discretize(equation, domain_h, [Vh, Vh]) - # ... - - # ... discretize norms - l2norm_h = discretize(l2norm, domain_h, Vh) - h1norm_h = discretize(h1norm, domain_h, Vh) - # ... - - # ... solve the discrete equation - x = equation_h.solve() - # ... - - # ... - phi = FemField( Vh, x ) - # ... - - # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) - # ... - - return l2_error, h1_error - - ############################################################################### # SERIAL TESTS ############################################################################### -def test_api_poisson_2d_dir0_1234(): - - from sympy.abc import x,y +#============================================================================== +# 2D Poisson's equation +#============================================================================== +def test_poisson_2d_dir0_1234(): solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*sin(pi*x)*sin(pi*y) @@ -349,10 +290,8 @@ def test_api_poisson_2d_dir0_1234(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_234_neu0_1(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_234_neu0_1(): solution = cos(0.5*pi*x)*sin(pi*y) f = (5./4.)*pi**2*solution @@ -369,10 +308,8 @@ def test_api_poisson_2d_dir0_234_neu0_1(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_134_neu0_2(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_134_neu0_2(): solution = sin(0.5*pi*x)*sin(pi*y) f = (5./4.)*pi**2*solution @@ -389,10 +326,8 @@ def test_api_poisson_2d_dir0_134_neu0_2(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_124_neu0_3(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_124_neu0_3(): solution = sin(pi*x)*cos(0.5*pi*y) f = (5./4.)*pi**2*solution @@ -409,10 +344,8 @@ def test_api_poisson_2d_dir0_124_neu0_3(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_123_neu0_4(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_123_neu0_4(): solution = sin(pi*x)*sin(0.5*pi*y) f = (5./4.)*pi**2*solution @@ -429,10 +362,8 @@ def test_api_poisson_2d_dir0_123_neu0_4(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_24_neu0_13(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_24_neu0_13(): solution = cos(0.5*pi*x)*cos(0.5*pi*y) f = (1./2.)*pi**2*solution @@ -449,10 +380,8 @@ def test_api_poisson_2d_dir0_24_neu0_13(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_13_neu0_24(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_13_neu0_24(): solution = sin(0.5*pi*x)*sin(0.5*pi*y) f = (1./2.)*pi**2*solution @@ -469,10 +398,8 @@ def test_api_poisson_2d_dir0_13_neu0_24(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_4_neu0_123(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_4_neu0_123(): solution = cos(pi*x)*cos(0.5*pi*y) f = 5./4.*pi**2*solution @@ -489,10 +416,8 @@ def test_api_poisson_2d_dir0_4_neu0_123(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_234_neui_1(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_234_neui_1(): solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*solution @@ -509,10 +434,8 @@ def test_api_poisson_2d_dir0_234_neui_1(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_134_neui_2(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_134_neui_2(): solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*solution @@ -529,10 +452,8 @@ def test_api_poisson_2d_dir0_134_neui_2(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_124_neui_3(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_124_neui_3(): solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*solution @@ -549,10 +470,8 @@ def test_api_poisson_2d_dir0_124_neui_3(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_123_neui_4(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_123_neui_4(): solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*solution @@ -569,10 +488,8 @@ def test_api_poisson_2d_dir0_123_neui_4(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir0_123_diri_4(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_123_diri_4(): solution = sin(pi * x) * sin(0.5*pi * y) f = 5/4*pi**2 * solution @@ -589,10 +506,8 @@ def test_api_poisson_2d_dir0_123_diri_4(): assert abs(l2_error - expected_l2_error) < 1.e-7 assert abs(h1_error - expected_h1_error) < 1.e-7 -#============================================================================== -def test_api_poisson_2d_dir0_13_diri_24(): - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_13_diri_24(): solution = sin(3*pi/2 * x) * sin(3*pi/2 * y) f = 9/2*pi**2 * solution @@ -609,16 +524,47 @@ def test_api_poisson_2d_dir0_13_diri_24(): assert abs(l2_error - expected_l2_error) < 1.e-7 assert abs(h1_error - expected_h1_error) < 1.e-7 -#============================================================================== +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_1234_user_function(): -def test_api_laplace_2d_neu(): + solution = sin(pi*x)*sin(pi*y) + + # ... + # User provides right-hand side in the form of a callable Python function: + def f(x, y): + from numpy import pi, sin + return 2*pi**2*sin(pi*x)*sin(pi*y) + + # Python function is converted to Sympy's "implemented function" and then + # called with symbolic arguments (x, y): + f = implemented_function('f', f)(x, y) + # ... + + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) + + expected_l2_error = 0.00021808678604760232 + expected_h1_error = 0.013023570720360362 - from sympy.abc import x,y + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#============================================================================== +# 2D "Laplace-like" equation +#============================================================================== +def test_laplace_2d_neu0_1234(): solution = cos(pi*x)*cos(pi*y) f = (2.*pi**2 + 1.)*solution - l2_error, h1_error = run_laplace_2d_neu(solution, f, ncells=[2**3,2**3], degree=[2,2]) + dir_zero_boundary = get_boundaries() + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error = run_laplace_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) expected_l2_error = 0.0002172846538950129 expected_h1_error = 0.012984852988125026 @@ -627,9 +573,9 @@ def test_api_laplace_2d_neu(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_biharmonic_2d_dir0_1234(): - - from sympy.abc import x,y +# 2D biharmonic equation +#============================================================================== +def test_biharmonic_2d_dir0_1234(): solution = sin(pi * x)**2 * sin(pi * y)**2 f = laplace(laplace(solution)) @@ -647,11 +593,9 @@ def test_api_biharmonic_2d_dir0_1234(): assert( abs(h1_error - expected_h1_error) < 1.e-7) assert( abs(h2_error - expected_h2_error) < 1.e-7) -#============================================================================== +#------------------------------------------------------------------------------ @pytest.mark.xfail -def test_api_biharmonic_2d_dir0_123_diri_4(): - - from sympy.abc import x,y +def test_biharmonic_2d_dir0_123_diri_4(): solution = sin(pi * x)**2 * sin(0.5*pi * y)**2 f = laplace(laplace(solution)) @@ -669,11 +613,9 @@ def test_api_biharmonic_2d_dir0_123_diri_4(): assert False -#============================================================================== +#------------------------------------------------------------------------------ @pytest.mark.xfail -def test_api_biharmonic_2d_dir0_13_diri_24(): - - from sympy.abc import x,y +def test_biharmonic_2d_dir0_13_diri_24(): solution = sin(3*pi/2 * x)**2 * sin(3*pi/2 * y)**2 f = laplace(laplace(solution)) @@ -691,41 +633,13 @@ def test_api_biharmonic_2d_dir0_13_diri_24(): assert False -#============================================================================== -def test_api_poisson_user_function_2d_dir_1(): - - from sympy.abc import x,y - - solution = sin(pi*x)*sin(pi*y) - - # ... - def f(x,y): - from numpy import pi - from numpy import cos - from numpy import sin - - value = 2*pi**2*sin(pi*x)*sin(pi*y) - return value - # ... - - l2_error, h1_error = run_poisson_user_function_2d_dir(f, solution, - ncells=[2**3,2**3], degree=[2,2]) - - expected_l2_error = 0.00021808678604760232 - expected_h1_error = 0.013023570720360362 - - assert( abs(l2_error - expected_l2_error) < 1.e-7) - assert( abs(h1_error - expected_h1_error) < 1.e-7) - ############################################################################### # PARALLEL TESTS ############################################################################### #============================================================================== @pytest.mark.parallel -def test_api_poisson_2d_dir0_1234_parallel(): - - from sympy.abc import x,y +def test_poisson_2d_dir0_1234_parallel(): solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*sin(pi*x)*sin(pi*y) From 4e5ce522aa957986dd11e9d57e0779860db0df74 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 28 Nov 2019 15:27:32 +0100 Subject: [PATCH 16/23] Reorganize 2D scalar tests w/ mapping --- .../api/tests/test_api_2d_scalar_mapping.py | 847 +++++++++--------- 1 file changed, 434 insertions(+), 413 deletions(-) diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py index de5d23ca0..51b98b4fb 100644 --- a/psydac/api/tests/test_api_2d_scalar_mapping.py +++ b/psydac/api/tests/test_api_2d_scalar_mapping.py @@ -13,6 +13,7 @@ from mpi4py import MPI from sympy import pi, cos, sin +from sympy.abc import x, y import pytest import os @@ -41,294 +42,250 @@ # ... #============================================================================== -def run_poisson_2d_dir(filename, solution, f, dir_boundary=None, comm=None): +def get_boundaries(*args): - # ... abstract model - domain = Domain.from_file(filename) - - V = ScalarFunctionSpace('V', domain) - - if dir_boundary is None: - B_dirichlet = domain.boundary - elif len(dir_boundary) == 1: - B_dirichlet = domain.get_boundary(**dir_boundary[0]) + if not args: + return () else: - B_dirichlet = Union(*[domain.get_boundary(**kw) for kw in dir_boundary]) - - x,y = domain.coordinates - - F = element_of(V, name='F') + assert all(1 <= a <= 4 for a in args) + assert len(set(args)) == len(args) - v = element_of(V, name='v') - u = element_of(V, name='u') + boundaries = {1: {'axis': 0, 'ext': -1}, + 2: {'axis': 0, 'ext': 1}, + 3: {'axis': 1, 'ext': -1}, + 4: {'axis': 1, 'ext': 1}} - int_0 = lambda expr: integral(domain , expr) - - expr = dot(grad(v), grad(u)) - a = BilinearForm((v,u), int_0(expr)) - - expr = f*v - l = LinearForm(v, int_0(expr)) - - error = F - solution - l2norm = Norm(error, domain, kind='l2') - h1norm = Norm(error, domain, kind='h1') - - bc = EssentialBC(u, 0, B_dirichlet) - equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc) - # ... - - # ... create the computational domain from a topological domain - domain_h = discretize(domain, filename=filename, comm=comm) - # ... - - # ... discrete spaces - Vh = discretize(V, domain_h) - # ... - - # ... dsicretize the equation using Dirichlet bc - equation_h = discretize(equation, domain_h, [Vh, Vh]) - # ... - - # ... discretize norms - l2norm_h = discretize(l2norm, domain_h, Vh) - h1norm_h = discretize(h1norm, domain_h, Vh) - # ... - - # ... solve the discrete equation - x = equation_h.solve() - # ... - - # ... - phi = FemField( Vh, x ) - # ... - - # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) - # ... - - return l2_error, h1_error + return tuple(boundaries[i] for i in args) #============================================================================== -def run_poisson_2d_dirneu(filename, solution, f, boundary, comm=None): +def run_poisson_2d(filename, solution, f, dir_zero_boundary, + dir_nonzero_boundary, neumann_boundary, comm=None): - assert( isinstance(boundary, (list, tuple)) ) + assert isinstance( dir_zero_boundary, (list, tuple)) + assert isinstance(dir_nonzero_boundary, (list, tuple)) + assert isinstance( neumann_boundary, (list, tuple)) - # ... abstract model + #+++++++++++++++++++++++++++++++ + # 1. Abstract model + #+++++++++++++++++++++++++++++++ domain = Domain.from_file(filename) - V = ScalarFunctionSpace('V', domain) - - B_neumann = [domain.get_boundary(**kw) for kw in boundary] - if len(B_neumann) == 1: - B_neumann = B_neumann[0] - - else: - B_neumann = Union(*B_neumann) - - x,y = domain.coordinates - - int_0 = lambda expr: integral(domain , expr) - int_1 = lambda expr: integral(B_neumann , expr) - - F = element_of(V, name='F') - - v = element_of(V, name='v') - u = element_of(V, name='u') + B_dirichlet_0 = Union(*[domain.get_boundary(**kw) for kw in dir_zero_boundary]) + B_dirichlet_i = Union(*[domain.get_boundary(**kw) for kw in dir_nonzero_boundary]) + B_neumann = Union(*[domain.get_boundary(**kw) for kw in neumann_boundary]) + V = ScalarFunctionSpace('V', domain) + u = element_of(V, name='u') + v = element_of(V, name='v') nn = NormalVector('nn') - expr = dot(grad(v), grad(u)) - a = BilinearForm((v,u), int_0(expr)) + # Bilinear form a: V x V --> R + a = BilinearForm((u, v), integral(domain, dot(grad(u), grad(v)))) - expr = f*v - l0 = LinearForm(v, int_0(expr)) + # Linear form l: V --> R + l0 = LinearForm(v, integral(domain, f * v)) + if B_neumann: + l1 = LinearForm(v, integral(B_neumann, v * dot(grad(solution), nn))) + l = LinearForm(v, l0(v) + l1(v)) + else: + l = l0 - expr = v*dot(grad(solution), nn) - l_B_neumann = LinearForm(v, int_1(expr)) + # Dirichlet boundary conditions + bc = [] + if B_dirichlet_0: bc += [EssentialBC(u, 0, B_dirichlet_0)] + if B_dirichlet_i: bc += [EssentialBC(u, solution, B_dirichlet_i)] - expr = l0(v) + l_B_neumann(v) - l = LinearForm(v, expr) + # Variational model + equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc) - error = F-solution + # Error norms + error = u - solution l2norm = Norm(error, domain, kind='l2') h1norm = Norm(error, domain, kind='h1') - B_dirichlet = domain.boundary.complement(B_neumann) - bc = EssentialBC(u, 0, B_dirichlet) + #+++++++++++++++++++++++++++++++ + # 2. Discretization + #+++++++++++++++++++++++++++++++ - equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc) - # ... - - # ... create the computational domain from a topological domain + # Create computational domain from topological domain domain_h = discretize(domain, filename=filename, comm=comm) - # ... - # ... discrete spaces + # Discrete spaces Vh = discretize(V, domain_h) - # ... - # ... dsicretize the equation using Dirichlet bc + # Discretize equation using Dirichlet bc equation_h = discretize(equation, domain_h, [Vh, Vh]) - # ... - # ... discretize norms + # Discretize error norms l2norm_h = discretize(l2norm, domain_h, Vh) h1norm_h = discretize(h1norm, domain_h, Vh) - # ... - # ... solve the discrete equation - x = equation_h.solve() - # ... + #+++++++++++++++++++++++++++++++ + # 3. Solution + #+++++++++++++++++++++++++++++++ - # ... - phi = FemField( Vh, x ) - # ... + # Solve linear system + X = equation_h.solve() + uh = FemField(Vh, X) - # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) - # ... + # Compute error norms + l2_error = l2norm_h.assemble(u=uh) + h1_error = h1norm_h.assemble(u=uh) return l2_error, h1_error #============================================================================== -def run_laplace_2d_neu(filename, solution, f, comm=None): - - # ... abstract model - domain = Domain.from_file(filename) - - V = ScalarFunctionSpace('V', domain) - - B_neumann = domain.boundary +def run_laplace_2d(filename, solution, f, dir_zero_boundary, + dir_nonzero_boundary, neumann_boundary, comm=None): - x,y = domain.coordinates + assert isinstance( dir_zero_boundary, (list, tuple)) + assert isinstance(dir_nonzero_boundary, (list, tuple)) + assert isinstance( neumann_boundary, (list, tuple)) - F = element_of(V, name='F') + #+++++++++++++++++++++++++++++++ + # 1. Abstract model + #+++++++++++++++++++++++++++++++ + domain = Domain.from_file(filename) - v = element_of(V, name='v') - u = element_of(V, name='u') + B_dirichlet_0 = Union(*[domain.get_boundary(**kw) for kw in dir_zero_boundary]) + B_dirichlet_i = Union(*[domain.get_boundary(**kw) for kw in dir_nonzero_boundary]) + B_neumann = Union(*[domain.get_boundary(**kw) for kw in neumann_boundary]) + V = ScalarFunctionSpace('V', domain) + u = element_of(V, name='u') + v = element_of(V, name='v') nn = NormalVector('nn') - int_0 = lambda expr: integral(domain , expr) - int_1 = lambda expr: integral(B_neumann , expr) - - expr = dot(grad(v), grad(u)) + v*u - a = BilinearForm((v,u), int_0(expr)) + # Bilinear form a: V x V --> R + a = BilinearForm((u, v), integral(domain, dot(grad(u), grad(v)) + u * v)) - expr = f*v - l0 = LinearForm(v, int_0(expr)) + # Linear form l: V --> R + l0 = LinearForm(v, integral(domain, f * v)) + if B_neumann: + l1 = LinearForm(v, integral(B_neumann, v * dot(grad(solution), nn))) + l = LinearForm(v, l0(v) + l1(v)) + else: + l = l0 - expr = v*dot(grad(solution), nn) - l_B_neumann = LinearForm(v, int_1(expr)) + # Dirichlet boundary conditions + bc = [] + if B_dirichlet_0: bc += [EssentialBC(u, 0, B_dirichlet_0)] + if B_dirichlet_i: bc += [EssentialBC(u, solution, B_dirichlet_i)] - expr = l0(v) + l_B_neumann(v) - l = LinearForm(v, expr) + # Variational model + equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc) - error = F-solution + # Error norms + error = u - solution l2norm = Norm(error, domain, kind='l2') h1norm = Norm(error, domain, kind='h1') - equation = find(u, forall=v, lhs=a(u,v), rhs=l(v)) - # ... + #+++++++++++++++++++++++++++++++ + # 2. Discretization + #+++++++++++++++++++++++++++++++ - # ... create the computational domain from a topological domain + # Create computational domain from topological domain domain_h = discretize(domain, filename=filename, comm=comm) - # ... - # ... discrete spaces + # Discrete spaces Vh = discretize(V, domain_h) - # ... - # ... dsicretize the equation using Dirichlet bc + # Discretize equation using Dirichlet bc equation_h = discretize(equation, domain_h, [Vh, Vh]) - # ... - # ... discretize norms + # Discretize error norms l2norm_h = discretize(l2norm, domain_h, Vh) h1norm_h = discretize(h1norm, domain_h, Vh) - # ... - # ... solve the discrete equation - x = equation_h.solve() - # ... + #+++++++++++++++++++++++++++++++ + # 3. Solution + #+++++++++++++++++++++++++++++++ - # ... - phi = FemField( Vh, x ) - # ... + # Solve linear system + X = equation_h.solve() + uh = FemField(Vh, X) - # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) - # ... + # Compute error norms + l2_error = l2norm_h.assemble(u=uh) + h1_error = h1norm_h.assemble(u=uh) return l2_error, h1_error #============================================================================== -def run_biharmonic_2d_dir(filename, solution, f, comm=None): +def run_biharmonic_2d_dir(filename, solution, f, dir_zero_boundary, + dir_nonzero_boundary, comm=None): + + assert isinstance( dir_zero_boundary, (list, tuple)) + assert isinstance(dir_nonzero_boundary, (list, tuple)) - # ... abstract model + #+++++++++++++++++++++++++++++++ + # 1. Abstract model + #+++++++++++++++++++++++++++++++ domain = Domain.from_file(filename) - V = ScalarFunctionSpace('V', domain) + B_dirichlet_0 = Union(*[domain.get_boundary(**kw) for kw in dir_zero_boundary]) + B_dirichlet_i = Union(*[domain.get_boundary(**kw) for kw in dir_nonzero_boundary]) - F = element_of(V, name='F') + V = ScalarFunctionSpace('V', domain) + u = element_of(V, name='u') + v = element_of(V, name='v') + nn = NormalVector('nn') - v = element_of(V, name='v') - u = element_of(V, name='u') + # Bilinear form a: V x V --> R + a = BilinearForm((u, v), integral(domain, laplace(u) * laplace(v))) - int_0 = lambda expr: integral(domain , expr) + # Linear form l: V --> R + l = LinearForm(v, integral(domain, f * v)) - expr = laplace(v) * laplace(u) - a = BilinearForm((v,u),int_0(expr)) + # Essential boundary conditions + dn = lambda a: dot(grad(a), nn) + bc = [] + if B_dirichlet_0: + bc += [EssentialBC( u , 0, B_dirichlet_0)] + bc += [EssentialBC(dn(u), 0, B_dirichlet_0)] + if B_dirichlet_i: + bc += [EssentialBC( u , solution , B_dirichlet_i)] + bc += [EssentialBC(dn(u), dn(solution), B_dirichlet_i)] - expr = f*v - l = LinearForm(v, int_0(expr)) + # Variational model + equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc) - error = F - solution + # Error norms + error = u - solution l2norm = Norm(error, domain, kind='l2') h1norm = Norm(error, domain, kind='h1') h2norm = Norm(error, domain, kind='h2') - nn = NormalVector('nn') - bc = [EssentialBC(u, 0, domain.boundary)] - bc += [EssentialBC(dot(grad(u), nn), 0, domain.boundary)] - equation = find(u, forall=v, lhs=a(u,v), rhs=l(v), bc=bc) - # ... + #+++++++++++++++++++++++++++++++ + # 2. Discretization + #+++++++++++++++++++++++++++++++ - # ... create the computational domain from a topological domain + # Create computational domain from topological domain domain_h = discretize(domain, filename=filename, comm=comm) - # ... - # ... discrete spaces + # Discrete spaces Vh = discretize(V, domain_h) - # ... - # ... dsicretize the equation using Dirichlet bc + # Discretize equation using Dirichlet bc equation_h = discretize(equation, domain_h, [Vh, Vh]) - # ... - # ... discretize norms + # Discretize error norms l2norm_h = discretize(l2norm, domain_h, Vh) h1norm_h = discretize(h1norm, domain_h, Vh) h2norm_h = discretize(h2norm, domain_h, Vh) - # ... - # ... solve the discrete equation - x = equation_h.solve() - # ... + #+++++++++++++++++++++++++++++++ + # 3. Solution + #+++++++++++++++++++++++++++++++ - # ... - phi = FemField( Vh, x ) - # ... + # Solve linear system + X = equation_h.solve() + uh = FemField(Vh, X) - # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) - h2_error = h2norm_h.assemble(F=phi) - # ... + # Compute error norms + l2_error = l2norm_h.assemble(u=uh) + h1_error = h1norm_h.assemble(u=uh) + h2_error = h2norm_h.assemble(u=uh) return l2_error, h1_error, h2_error @@ -337,15 +294,20 @@ def run_biharmonic_2d_dir(filename, solution, f, comm=None): ############################################################################### #============================================================================== -def test_api_poisson_2d_dir_identity(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') - - from sympy.abc import x,y +# 2D Poisson's equation with identity map +#============================================================================== +def test_poisson_2d_identity_dir0_1234(): + filename = os.path.join(mesh_dir, 'identity_2d.h5') solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*sin(pi*x)*sin(pi*y) - l2_error, h1_error = run_poisson_2d_dir(filename, solution, f) + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) expected_l2_error = 0.00021808678604159413 expected_h1_error = 0.013023570720357957 @@ -353,182 +315,241 @@ def test_api_poisson_2d_dir_identity(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir_collela(): - filename = os.path.join(mesh_dir, 'collela_2d.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_234_neu0_1(): - from sympy.abc import x,y + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = cos(0.5*pi*x)*sin(pi*y) + f = (5./4.)*pi**2*solution - solution = sin(pi*x)*sin(pi*y) - f = 2*pi**2*sin(pi*x)*sin(pi*y) + dir_zero_boundary = get_boundaries(2, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1) - l2_error, h1_error = run_poisson_2d_dir(filename, solution, f) + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.03032933682661518 - expected_h1_error = 0.41225081526853247 + expected_l2_error = 0.00015546057795986509 + expected_h1_error = 0.009269302784527035 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir_quarter_annulus(): - filename = os.path.join(mesh_dir, 'quarter_annulus.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_134_neu0_2(): - from sympy.abc import x,y + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(0.5*pi*x)*sin(pi*y) + f = (5./4.)*pi**2*solution - c = pi / (1. - 0.5**2) - r2 = 1. - x**2 - y**2 - solution = x*y*sin(c * r2) - f = 4.*c**2*x*y*(x**2 + y**2)*sin(c * r2) + 12.*c*x*y*cos(c * r2) + dir_zero_boundary = get_boundaries(1, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(2) - l2_error, h1_error = run_poisson_2d_dir(filename, solution, f) + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.00010289930281268989 - expected_h1_error = 0.009473407914765117 + expected_l2_error = 0.00015546057795095866 + expected_h1_error = 0.009269302784528054 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dir_circle(): - filename = os.path.join(mesh_dir, 'circle.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_124_neu0_3(): - from sympy.abc import x,y + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(pi*x)*cos(0.5*pi*y) + f = (5./4.)*pi**2*solution - solution = (1 - (x**2 + y**2)) * cos(2*pi*x) * cos(2*pi*y) - f = -laplace(solution) + dir_zero_boundary = get_boundaries(1, 2, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(3) - dir_boundary = [{'axis': 0, 'ext': 1}] - l2_error, h1_error = run_poisson_2d_dir(filename, solution, f, dir_boundary) + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.0015245737751297718 - expected_h1_error = 0.06653900724243668 + expected_l2_error = 0.00015546057796188848 + expected_h1_error = 0.009269302784527448 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dirneu_identity_1(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_123_neu0_4(): - from sympy.abc import x,y - - solution = cos(0.5*pi*x)*sin(pi*y) + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(pi*x)*sin(0.5*pi*y) f = (5./4.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, - [{'axis': 0, 'ext': -1}]) + dir_zero_boundary = get_boundaries(1, 2, 3) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(4) - expected_l2_error = 0.00015546057795986509 - expected_h1_error = 0.009269302784527035 + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.00015546057795073548 + expected_h1_error = 0.009269302784522822 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dirneu_identity_2(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_24_neu0_13(): - from sympy.abc import x,y + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = cos(0.5*pi*x)*cos(0.5*pi*y) + f = (1./2.)*pi**2*solution - solution = sin(0.5*pi*x)*sin(pi*y) - f = (5./4.)*pi**2*solution + dir_zero_boundary = get_boundaries(2, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1, 3) - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, - [{'axis': 0, 'ext': 1}]) + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.00015546057795095866 - expected_h1_error = 0.009269302784528054 + expected_l2_error = 2.6119892693464717e-05 + expected_h1_error = 0.0016032430287989195 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dirneu_identity_3(): +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_4_neu0_123(): + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = cos(pi*x)*cos(0.5*pi*y) + f = 5./4.*pi**2*solution - from sympy.abc import x,y + dir_zero_boundary = get_boundaries(4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1, 2, 3) - solution = sin(pi*x)*cos(0.5*pi*y) - f = (5./4.)*pi**2*solution - - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, - [{'axis': 1, 'ext': -1}]) + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.00015546057796188848 - expected_h1_error = 0.009269302784527448 + expected_l2_error = 0.00015492540684276186 + expected_h1_error = 0.009242166615517364 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dirneu_identity_4(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_234_neui_1(): - from sympy.abc import x,y + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution - solution = sin(pi*x)*sin(0.5*pi*y) - f = (5./4.)*pi**2*solution + dir_zero_boundary = get_boundaries(2, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1) - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, - [{'axis': 1, 'ext': 1}]) + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.00015546057795073548 - expected_h1_error = 0.009269302784522822 + expected_l2_error = 0.00021786960671761908 + expected_h1_error = 0.01302350067761177 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dirneu_identity_13(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_134_neui_2(): - from sympy.abc import x,y + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution - solution = cos(0.5*pi*x)*cos(0.5*pi*y) - f = (1./2.)*pi**2*solution + dir_zero_boundary = get_boundaries(1, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(2) - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, - [{'axis': 0, 'ext': -1}, - {'axis': 1, 'ext': -1}]) + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 2.6119892693464717e-05 - expected_h1_error = 0.0016032430287989195 + expected_l2_error = 0.00021786960671761908 + expected_h1_error = 0.01302350067761177 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dirneu_identity_123(): +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_124_neui_3(): + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution - from sympy.abc import x,y + dir_zero_boundary = get_boundaries(1, 2, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(3) - solution = cos(pi*x)*cos(0.5*pi*y) - f = 5./4.*pi**2*solution + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, - [{'axis': 0, 'ext': -1}, - {'axis': 0, 'ext': 1}, - {'axis': 1, 'ext': -1}]) + expected_l2_error = 0.00021786960671761908 + expected_h1_error = 0.01302350067761177 - expected_l2_error = 0.00015492540684276186 - expected_h1_error = 0.009242166615517364 + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_123_neui_4(): + + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution + + dir_zero_boundary = get_boundaries(1, 2, 3) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(4) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.00021786960671761908 + expected_h1_error = 0.01302350067761177 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) +#============================================================================== +# 2D Poisson's equation with "Collela" map +#============================================================================== +def test_poisson_2d_collela_dir0_1234(): -##============================================================================== -def test_api_poisson_2d_dirneu_collela_1(): filename = os.path.join(mesh_dir, 'collela_2d.h5') + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*sin(pi*x)*sin(pi*y) + + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.03032933682661518 + expected_h1_error = 0.41225081526853247 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_234_neu0_1(): + filename = os.path.join(mesh_dir, 'collela_2d.h5') solution = sin(0.25*pi*(x-1))*sin(pi*y) f = (17./16.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, [{'axis': 0, 'ext': -1}]) + dir_zero_boundary = get_boundaries(2, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) expected_l2_error = 0.013540717397796734 expected_h1_error = 0.19789463571596025 @@ -536,16 +557,19 @@ def test_api_poisson_2d_dirneu_collela_1(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dirneu_collela_2(): - filename = os.path.join(mesh_dir, 'collela_2d.h5') - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_134_neu0_2(): + filename = os.path.join(mesh_dir, 'collela_2d.h5') solution = sin(0.25*pi*(x+1.))*sin(pi*y) f = (17./16.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, [{'axis': 0, 'ext': 1}]) + dir_zero_boundary = get_boundaries(1, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(2) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) expected_l2_error = 0.012890849094111699 expected_h1_error = 0.19553563279728328 @@ -553,16 +577,19 @@ def test_api_poisson_2d_dirneu_collela_2(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -##============================================================================== -def test_api_poisson_2d_dirneu_collela_3(): - filename = os.path.join(mesh_dir, 'collela_2d.h5') - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_124_neu0_3(): + filename = os.path.join(mesh_dir, 'collela_2d.h5') solution = sin(0.25*pi*(y-1))*sin(pi*x) f = (17./16.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, [{'axis': 1, 'ext': -1}]) + dir_zero_boundary = get_boundaries(1, 2, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(3) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) expected_l2_error = 0.013540717397817427 expected_h1_error = 0.19789463571595994 @@ -570,16 +597,19 @@ def test_api_poisson_2d_dirneu_collela_3(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_poisson_2d_dirneu_collela_4(): - filename = os.path.join(mesh_dir, 'collela_2d.h5') - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_123_neu0_4(): + filename = os.path.join(mesh_dir, 'collela_2d.h5') solution = sin(0.25*pi*(y+1.))*sin(pi*x) f = (17./16.)*pi**2*solution - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, [{'axis': 1, 'ext': 1}]) + dir_zero_boundary = get_boundaries(1, 2, 3) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(4) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) expected_l2_error = 0.012890849094111942 expected_h1_error = 0.19553563279728325 @@ -588,104 +618,107 @@ def test_api_poisson_2d_dirneu_collela_4(): assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dir_zero_neu_nonzero_identity_1(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') - - from sympy.abc import x,y - - solution = sin(pi*x)*sin(pi*y) - f = 2*pi**2*solution - - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, - [{'axis': 0, 'ext': -1}]) - - expected_l2_error = 0.00021786960671761908 - expected_h1_error = 0.01302350067761177 - - assert( abs(l2_error - expected_l2_error) < 1.e-7) - assert( abs(h1_error - expected_h1_error) < 1.e-7) - +# 2D Poisson's equation on quarter annulus #============================================================================== -def test_api_poisson_2d_dir_zero_neu_nonzero_identity_2(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') +def test_poisson_2d_quarter_annulus_dir0_1234(): - from sympy.abc import x,y + filename = os.path.join(mesh_dir, 'quarter_annulus.h5') + c = pi / (1. - 0.5**2) + r2 = 1. - x**2 - y**2 + solution = x*y*sin(c * r2) + f = 4.*c**2*x*y*(x**2 + y**2)*sin(c * r2) + 12.*c*x*y*cos(c * r2) - solution = sin(pi*x)*sin(pi*y) - f = 2*pi**2*solution + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries() - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, - [{'axis': 0, 'ext': 1}]) + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.00021786960671761908 - expected_h1_error = 0.01302350067761177 + expected_l2_error = 0.00010289930281268989 + expected_h1_error = 0.009473407914765117 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dir_zero_neu_nonzero_identity_3(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') +# 2D Poisson's equation on circle +#============================================================================== +def test_poisson_2d_circle_dir0(): - from sympy.abc import x,y + filename = os.path.join(mesh_dir, 'circle.h5') + solution = (1 - (x**2 + y**2)) * cos(2*pi*x) * cos(2*pi*y) + f = -laplace(solution) - solution = sin(pi*x)*sin(pi*y) - f = 2*pi**2*solution + dir_zero_boundary = get_boundaries(2) # only boundary is at r = r_max + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries() - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, - [{'axis': 1, 'ext': -1}]) + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.00021786960671761908 - expected_h1_error = 0.01302350067761177 + expected_l2_error = 0.0015245737751297718 + expected_h1_error = 0.06653900724243668 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_poisson_2d_dir_zero_neu_nonzero_identity_4(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') +# 2D "Laplace-like" equation +#============================================================================== +def test_laplace_2d_identity_neu0_1234(): - from sympy.abc import x,y + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = cos(pi*x)*cos(pi*y) + f = (2.*pi**2 + 1.)*solution - solution = sin(pi*x)*sin(pi*y) - f = 2*pi**2*solution + dir_zero_boundary = get_boundaries() + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1, 2, 3, 4) - l2_error, h1_error = run_poisson_2d_dirneu(filename, solution, f, - [{'axis': 1, 'ext': 1}]) + l2_error, h1_error = run_laplace_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.00021786960671761908 - expected_h1_error = 0.01302350067761177 + expected_l2_error = 0.00021728465388208586 + expected_h1_error = 0.012984852988123631 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) -#============================================================================== -def test_api_laplace_2d_neu_identity(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_laplace_2d_collela_neu0_1234(): + filename = os.path.join(mesh_dir, 'collela_2d.h5') solution = cos(pi*x)*cos(pi*y) f = (2.*pi**2 + 1.)*solution - l2_error, h1_error = run_laplace_2d_neu(filename, solution, f) + dir_zero_boundary = get_boundaries() + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1, 2, 3, 4) - expected_l2_error = 0.00021728465388208586 - expected_h1_error = 0.012984852988123631 + l2_error, h1_error = run_laplace_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.029603335241478155 + expected_h1_error = 0.4067746760978581 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) #============================================================================== -def test_api_biharmonic_2d_dir_identity(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') - - from sympy.abc import x,y +# 2D biharmonic equation +#============================================================================== +def test_biharmonic_2d_identity_dir0_1234(): + filename = os.path.join(mesh_dir, 'identity_2d.h5') solution = (sin(pi*x)*sin(pi*y))**2 f = laplace(laplace(solution)) - l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f) + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary) expected_l2_error = 0.015086415626060034 expected_h1_error = 0.08773346232941553 @@ -695,16 +728,18 @@ def test_api_biharmonic_2d_dir_identity(): assert( abs(h1_error - expected_h1_error) < 1.e-7) assert( abs(h2_error - expected_h2_error) < 1.e-7) -#============================================================================== -def test_api_biharmonic_2d_dir_collela(): - filename = os.path.join(mesh_dir, 'collela_2d.h5') - - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_biharmonic_2d_collela_dir0_1234(): + filename = os.path.join(mesh_dir, 'collela_2d.h5') solution = (cos(pi*x/2)*cos(pi*y/2))**2 f = laplace(laplace(solution)) - l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f) + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary) expected_l2_error = 0.10977627980052021 expected_h1_error = 0.32254511059711766 @@ -714,19 +749,21 @@ def test_api_biharmonic_2d_dir_collela(): assert( abs(h1_error - expected_h1_error) < 1.e-7) assert( abs(h2_error - expected_h2_error) < 1.e-7) -#============================================================================== -def test_api_biharmonic_2d_dir_quarter_annulus(): - filename = os.path.join(mesh_dir, 'quarter_annulus.h5') +#------------------------------------------------------------------------------ +def test_biharmonic_2d_quarter_annulus_dir0_1234(): - from sympy.abc import x,y - - r_in = 0.5 - r_out = 1 - kappa = 1 / 0.00643911127175763 + filename = os.path.join(mesh_dir, 'quarter_annulus.h5') + r_in = 0.5 + r_out = 1 + kappa = 1 / 0.00643911127175763 solution = kappa * (x * y * (x**2 + y**2 - r_in**2) * (x**2 + y**2 - r_out**2))**2 f = laplace(laplace(solution)) - l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f) + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + + l2_error, h1_error, h2_error = run_biharmonic_2d_dir(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary) expected_l2_error = 0.016730298635551484 expected_h1_error = 0.21243295522291714 @@ -736,40 +773,25 @@ def test_api_biharmonic_2d_dir_quarter_annulus(): assert( abs(h1_error - expected_h1_error) < 1.e-7) assert( abs(h2_error - expected_h2_error) < 1.e-7) -##============================================================================== -## TODO DEBUG, not working since merge with devel -#def test_api_laplace_2d_neu_collela(): -# filename = os.path.join(mesh_dir, 'collela_2d.h5') -# -# from sympy.abc import x,y -# -# solution = cos(pi*x)*cos(pi*y) -# f = (2.*pi**2 + 1.)*solution -# -# l2_error, h1_error = run_laplace_2d_neu(filename, solution, f) -# -# expected_l2_error = 0.029603335241478155 -# expected_h1_error = 0.4067746760978581 -# -# assert( abs(l2_error - expected_l2_error) < 1.e-7) -# assert( abs(h1_error - expected_h1_error) < 1.e-7) - ############################################################################### # PARALLEL TESTS ############################################################################### #============================================================================== @pytest.mark.parallel -def test_api_poisson_2d_dir_identity_parallel(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') - - from sympy.abc import x,y +def test_poisson_2d_identity_dir0_1234_parallel(): + filename = os.path.join(mesh_dir, 'identity_2d.h5') solution = sin(pi*x)*sin(pi*y) f = 2*pi**2*sin(pi*x)*sin(pi*y) - l2_error, h1_error = run_poisson_2d_dir(filename, solution, f, - comm=MPI.COMM_WORLD) + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary, + comm=MPI.COMM_WORLD) expected_l2_error = 0.00021808678604159413 expected_h1_error = 0.013023570720357957 @@ -777,7 +799,6 @@ def test_api_poisson_2d_dir_identity_parallel(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) - #============================================================================== # CLEAN UP SYMPY NAMESPACE #============================================================================== From fdbac701e361f948ff37f9c42886b34c44a4266a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 28 Nov 2019 18:02:15 +0100 Subject: [PATCH 17/23] New 2D tests w/ mapping & inhomog. Dirichlet/Neumann BCs --- .../api/tests/test_api_2d_scalar_mapping.py | 274 +++++++++++++++++- 1 file changed, 270 insertions(+), 4 deletions(-) diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py index 51b98b4fb..83c75e819 100644 --- a/psydac/api/tests/test_api_2d_scalar_mapping.py +++ b/psydac/api/tests/test_api_2d_scalar_mapping.py @@ -8,8 +8,11 @@ # larger square [-1, 1]^2, with deformations going as sin(pi x) * sin(pi y) # # - 'quarter_annulus.h5' is a NURBS transformation from the unit square [0, 1]^2 -# to the quarter annulus in the lower-left quadrant of the Cartesian place -# (hence both x and y are negative), with r_min = 0.5 and r_max = 0.5 +# to the quarter annulus in the lower-left quadrant of the Cartesian plane +# (hence both x and y are negative), with r_min = 0.5 and r_max = 1 +# +# Please note that the logical coordinates (x1, x2) correspond to the polar +# coordinates (r, theta), but with reversed order: hence x1=theta and x2=r from mpi4py import MPI from sympy import pi, cos, sin @@ -515,6 +518,46 @@ def test_poisson_2d_identity_dir0_123_neui_4(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_123_diri_4(): + + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(pi * x) * sin(0.5*pi * y) + f = 5/4*pi**2 * solution + + dir_zero_boundary = get_boundaries(1, 2, 3) + dir_nonzero_boundary = get_boundaries(4) + neumann_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.0001529221571156830 + expected_h1_error = 0.009293161646612863 + + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_13_diri_24(): + + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(3*pi/2 * x) * sin(3*pi/2 * y) + f = 9/2*pi**2 * solution + + dir_zero_boundary = get_boundaries(1, 3) + dir_nonzero_boundary = get_boundaries(2, 4) + neumann_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.0007786454571731944 + expected_h1_error = 0.0449669071240554 + + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + #============================================================================== # 2D Poisson's equation with "Collela" map #============================================================================== @@ -617,6 +660,146 @@ def test_poisson_2d_collela_dir0_123_neu0_4(): assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_234_neui_1(): + + filename = os.path.join(mesh_dir, 'collela_2d.h5') + solution = sin(pi/3 * (1 - x)) * cos(pi/2 * y) + f = (13/36)*pi**2 * solution + + dir_zero_boundary = get_boundaries(2, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.06606115297958418 + expected_h1_error = 0.1808088420439245 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_134_neui_2(): + + filename = os.path.join(mesh_dir, 'collela_2d.h5') + solution = sin(pi/3 * (1 + x)) * cos(pi/2 * y) + f = (13/36)*pi**2 * solution + + dir_zero_boundary = get_boundaries(1, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(2) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.06581556358643213 + expected_h1_error = 0.1810360332114653 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_124_neui_3(): + + filename = os.path.join(mesh_dir, 'collela_2d.h5') + solution = cos(pi/2 * x) * sin(pi/3 * (1 - y)) + f = (13/36)*pi**2 * solution + + dir_zero_boundary = get_boundaries(1, 2, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(3) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.06606115297960084 + expected_h1_error = 0.1808088420440391 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_123_neui_4(): + + filename = os.path.join(mesh_dir, 'collela_2d.h5') + solution = cos(pi/2 * x) * sin(pi/3 * (1 + y)) + f = (13/36)*pi**2 * solution + + dir_zero_boundary = get_boundaries(1, 2, 3) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(4) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.06581556358642991 + expected_h1_error = 0.1810360332114519 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_123_diri_4(): + + filename = os.path.join(mesh_dir, 'collela_2d.h5') + solution = cos(pi/2 * x) * sin(pi/3 * (1 + y)) + f = (13/36)*pi**2 * solution + + dir_zero_boundary = get_boundaries(1, 2, 3) + dir_nonzero_boundary = get_boundaries(4) + neumann_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.0025850223987204306 + expected_h1_error = 0.04401691486495642 + + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_13_diri_24(): + + filename = os.path.join(mesh_dir, 'collela_2d.h5') + solution = sin(pi/3 * (1 + x)) * sin(pi/3 * (1 + y)) + f = (2/9)*pi**2 * solution + + dir_zero_boundary = get_boundaries(1, 3) + dir_nonzero_boundary = get_boundaries(2, 4) + neumann_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.0012801077606328381 + expected_h1_error = 0.02314405549486328 + + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_diri_1234(): + + filename = os.path.join(mesh_dir, 'collela_2d.h5') + solution = cos(pi/3 * x) * cos(pi/3 * y) + f = (2/9)*pi**2 * solution + + dir_zero_boundary = get_boundaries() + dir_nonzero_boundary = get_boundaries(1, 2, 3, 4) + neumann_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.0014604231101091047 + expected_h1_error = 0.025023352115363873 + + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + #============================================================================== # 2D Poisson's equation on quarter annulus #============================================================================== @@ -638,8 +821,91 @@ def test_poisson_2d_quarter_annulus_dir0_1234(): expected_l2_error = 0.00010289930281268989 expected_h1_error = 0.009473407914765117 - assert( abs(l2_error - expected_l2_error) < 1.e-7) - assert( abs(h1_error - expected_h1_error) < 1.e-7) + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + +#------------------------------------------------------------------------------ +def test_poisson_2d_quarter_annulus_dir0_12_diri_34(): + + filename = os.path.join(mesh_dir, 'quarter_annulus.h5') + solution = sin(pi * x) * sin(pi * y) + f = 2*pi**2 * solution + + dir_zero_boundary = get_boundaries(1, 2) + dir_nonzero_boundary = get_boundaries(3, 4) + neumann_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.0005982761090480573 + expected_h1_error = 0.021271053089631443 + + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + +#------------------------------------------------------------------------------ +def test_poisson_2d_quarter_annulus_diri_1234(): + + filename = os.path.join(mesh_dir, 'quarter_annulus.h5') + solution = sin(pi*x + pi/4) * sin(pi*y + pi/4) + f = 2*pi**2 * solution + + dir_zero_boundary = get_boundaries() + dir_nonzero_boundary = get_boundaries(1, 2, 3, 4) + neumann_boundary = get_boundaries() + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.0006536882827670037 + expected_h1_error = 0.02592026831798558 + + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + +#------------------------------------------------------------------------------ +def test_poisson_2d_quarter_annulus_diri_34_neui_12(): + + filename = os.path.join(mesh_dir, 'quarter_annulus.h5') + solution = sin(pi*x + pi/4) * sin(pi*y + pi/4) + f = 2*pi**2 * solution + + dir_zero_boundary = get_boundaries() + dir_nonzero_boundary = get_boundaries(3, 4) + neumann_boundary = get_boundaries(1, 2) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + expected_l2_error = 0.005574858010884482 + expected_h1_error = 0.09514303219069051 + + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 + +#------------------------------------------------------------------------------ +@pytest.mark.xfail +def test_poisson_2d_quarter_annulus_diri_12_neui_34(): + + filename = os.path.join(mesh_dir, 'quarter_annulus.h5') + solution = sin(pi*x + pi/4) * sin(pi*y + pi/4) + f = 2*pi**2 * solution + + dir_zero_boundary = get_boundaries() + dir_nonzero_boundary = get_boundaries(1, 2) + neumann_boundary = get_boundaries(3, 4) + + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + print() + print() + print(l2_error) + print(h1_error) + print() + + assert False #============================================================================== # 2D Poisson's equation on circle From f821e44ac2557fb550168797fd8397802613830c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 28 Nov 2019 19:21:30 +0100 Subject: [PATCH 18/23] Bugfix: normalize vector normal to boundary --- psydac/api/ast/utilities.py | 34 +++++++++++-------- .../api/tests/test_api_2d_scalar_mapping.py | 31 ++++++++--------- 2 files changed, 34 insertions(+), 31 deletions(-) diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index 461f8ee7b..442f67a92 100644 --- a/psydac/api/ast/utilities.py +++ b/psydac/api/ast/utilities.py @@ -705,25 +705,31 @@ def compute_normal_vector(parent_namespace, vector, boundary, mapping=None): else: raise TypeError(boundary) - + # If there is no mapping, normal vector has only one non-zero component, + # which is +1 or -1 according to the orientation of the boundary. if mapping is None: - values = [(ext if i==axis else 0) for i, v in enumerate(vector)] - - else: - # Given the Jacobian matrix J, we need to extract the (i=axis) row - # of J^(-1). For efficiency we separately compute det(J) and the - # cofactors C[i=0...dim] of the (j=axis) column of J. - # NOTE: we also change the orientation according to 'ext' - J = SymbolicExpr(mapping.jacobian) - det = J.det() - cof = [J.cofactor(i, j=axis) for i in range(J.shape[0])] - - inv_jac = parent_namespace['inv_jac'] - values = [ext * inv_jac * cof[i] for i in range(J.shape[0])] + return [Assign(v, ext if i==axis else 0) for i, v in enumerate(vector)] + + # Given the Jacobian matrix J, we need to extract the (i=axis) row of + # J^(-1) and then normalize it. We recall that J^(-1)[i, j] is equal to + # the cofactor of J[i, j] divided by det(J). For efficiency we only + # compute the cofactors C[i=0:dim] of the (j=axis) column of J, and we + # do not divide them by det(J) because the normal vector will need to + # be normalized anyway. + # + # NOTE: we also change the vector orientation according to 'ext' + J = SymbolicExpr(mapping.jacobian) + values = [ext * J.cofactor(i, j=axis) for i in range(J.shape[0])] # Create statements for computing normal vector components stmts = [Assign(lhs, rhs) for lhs, rhs in zip(vector, values)] + # Normalize vector + inv_norm_variable = Symbol('inv_norm') + inv_norm_value = 1 / sympy_sqrt(sum(v**2 for v in values)) + stmts += [Assign(inv_norm_variable, inv_norm_value)] + stmts += [AugAssign(v, '*', inv_norm_variable) for v in vector] + return stmts #============================================================================== diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py index 83c75e819..d4fa4a753 100644 --- a/psydac/api/tests/test_api_2d_scalar_mapping.py +++ b/psydac/api/tests/test_api_2d_scalar_mapping.py @@ -674,8 +674,8 @@ def test_poisson_2d_collela_dir0_234_neui_1(): l2_error, h1_error = run_poisson_2d(filename, solution, f, dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.06606115297958418 - expected_h1_error = 0.1808088420439245 + expected_l2_error = 0.002701799327716594 + expected_h1_error = 0.043091889730461796 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) @@ -694,8 +694,8 @@ def test_poisson_2d_collela_dir0_134_neui_2(): l2_error, h1_error = run_poisson_2d(filename, solution, f, dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.06581556358643213 - expected_h1_error = 0.1810360332114653 + expected_l2_error = 0.0026061796931093253 + expected_h1_error = 0.04400143055955451 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) @@ -714,8 +714,8 @@ def test_poisson_2d_collela_dir0_124_neui_3(): l2_error, h1_error = run_poisson_2d(filename, solution, f, dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.06606115297960084 - expected_h1_error = 0.1808088420440391 + expected_l2_error = 0.002701799327724046 + expected_h1_error = 0.04309188973046055 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) @@ -734,8 +734,8 @@ def test_poisson_2d_collela_dir0_123_neui_4(): l2_error, h1_error = run_poisson_2d(filename, solution, f, dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.06581556358642991 - expected_h1_error = 0.1810360332114519 + expected_l2_error = 0.0026061796931066174 + expected_h1_error = 0.04400143055955377 assert( abs(l2_error - expected_l2_error) < 1.e-7) assert( abs(h1_error - expected_h1_error) < 1.e-7) @@ -878,14 +878,13 @@ def test_poisson_2d_quarter_annulus_diri_34_neui_12(): l2_error, h1_error = run_poisson_2d(filename, solution, f, dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - expected_l2_error = 0.005574858010884482 - expected_h1_error = 0.09514303219069051 + expected_l2_error = 0.0006527836834289991 + expected_h1_error = 0.025919435390680808 assert abs(l2_error - expected_l2_error) < 1.e-7 assert abs(h1_error - expected_h1_error) < 1.e-7 #------------------------------------------------------------------------------ -@pytest.mark.xfail def test_poisson_2d_quarter_annulus_diri_12_neui_34(): filename = os.path.join(mesh_dir, 'quarter_annulus.h5') @@ -899,13 +898,11 @@ def test_poisson_2d_quarter_annulus_diri_12_neui_34(): l2_error, h1_error = run_poisson_2d(filename, solution, f, dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - print() - print() - print(l2_error) - print(h1_error) - print() + expected_l2_error = 0.0006663772402662598 + expected_h1_error = 0.025906309642232804 - assert False + assert abs(l2_error - expected_l2_error) < 1.e-7 + assert abs(h1_error - expected_h1_error) < 1.e-7 #============================================================================== # 2D Poisson's equation on circle From fe0bd7d324b497223776b9adb8171c994537d79d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 3 Dec 2019 08:57:56 +0100 Subject: [PATCH 19/23] Improve handling of settings to linear solver --- psydac/api/discretization.py | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index ec6d81298..5d0f322f3 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -200,7 +200,10 @@ def assemble(self, **kwargs): self._linear_system = LinearSystem(M, rhs) def solve(self, **kwargs): - settings = kwargs.pop('settings', _default_solver.copy()) + + settings = _default_solver.copy() + settings.update(kwargs.pop('settings', {})) + rhs = kwargs.pop('rhs', None) if rhs: kwargs['assemble_rhs'] = False @@ -259,11 +262,15 @@ def solve(self, **kwargs): rhs_expr = sum(integral(i.boundary, product(i.rhs, factor(i))) for i in idbcs) equation = find(u, forall=v, lhs=lhs_expr, rhs=rhs_expr) - # Discretize weak form and find inhomogeneous solution + # Discretize weak form domain_h = self.domain Vh = self.trial_space equation_h = discretize(equation, domain_h, [Vh, Vh]) - X = equation_h.solve() + + # Find inhomogeneous solution (use CG as system is symmetric) + loc_settings = settings.copy() + loc_settings['solver'] = 'cg' + X = equation_h.solve(**loc_settings) # Use inhomogeneous solution as initial guess to solver settings['x0'] = X From f86fa079a6b08cd850855ff00f5511993c708018 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Thu, 19 Dec 2019 18:39:16 +0100 Subject: [PATCH 20/23] refactor: Add constructor 'Geometry.from_topological_domain' --- psydac/api/discretization.py | 42 +++++++++++++----------------------- psydac/cad/geometry.py | 29 +++++++------------------ 2 files changed, 23 insertions(+), 48 deletions(-) diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index 5d0f322f3..33b26f2f8 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -391,7 +391,7 @@ def discretize_derham(V, domain_h, *args, **kwargs): return DiscreteDerham(*spaces) #============================================================================== # TODO multi patch -# TODO bounds and knots +# TODO knots def discretize_space(V, domain_h, *args, **kwargs): degree = kwargs.pop('degree', None) normalize = kwargs.pop('normalize', True) @@ -424,13 +424,16 @@ def discretize_space(V, domain_h, *args, **kwargs): elif not( degree is None ): assert(hasattr(domain_h, 'ncells')) - ncells = domain_h.ncells + ncells = domain_h.ncells + min_coords = domain_h.domain.min_coords + max_coords = domain_h.domain.max_coords assert(isinstance( degree, (list, tuple) )) assert( len(degree) == ldim ) # Create uniform grid - grids = [np.linspace( 0., 1., num=ne+1 ) for ne in ncells] + grids = [np.linspace(xmin, xmax, num=ne + 1) + for xmin, xmax, ne in zip(min_coords, max_coords, ncells)] # Create 1D finite element spaces and precompute quadrature data spaces = [SplineSpace( p, grid=grid ) for p,grid in zip(degree, grids)] @@ -495,34 +498,19 @@ def discretize_space(V, domain_h, *args, **kwargs): return Vh #============================================================================== -def discretize_domain(domain, *args, **kwargs): - filename = kwargs.pop('filename', None) - ncells = kwargs.pop('ncells', None) - comm = kwargs.pop('comm', None) +def discretize_domain(domain, *, filename=None, ncells=None, comm=None): - if not( ncells is None ): - dtype = domain.dtype + if not (filename or ncells): + raise ValueError("Must provide either 'filename' or 'ncells'") - if dtype['type'].lower() == 'line' : - return Geometry.as_line(ncells, comm=comm) + elif filename and ncells: + raise ValueError("Cannot provide both 'filename' and 'ncells'") - elif dtype['type'].lower() == 'square' : - return Geometry.as_square(ncells, comm=comm) - - elif dtype['type'].lower() == 'cube' : - return Geometry.as_cube(ncells, comm=comm) - - else: - msg = 'no corresponding discrete geometry is available, given {}' - msg = msg.format(dtype['type']) - - raise ValueError(msg) - - elif not( filename is None ): - geometry = Geometry(filename=filename, comm=comm) - - return geometry + elif filename: + return Geometry(filename=filename, comm=comm) + elif ncells: + return Geometry.from_topological_domain(domain, ncells, comm) #============================================================================== def discretize(a, *args, **kwargs): diff --git a/psydac/cad/geometry.py b/psydac/cad/geometry.py index bbccfab39..d88604a80 100644 --- a/psydac/cad/geometry.py +++ b/psydac/cad/geometry.py @@ -89,34 +89,21 @@ def from_discrete_mapping( cls, mapping, comm=None ): #-------------------------------------------------------------------------- # Option [3]: discrete topological line/square/cube #-------------------------------------------------------------------------- - # TODO shall we create a discrete line/square/cube? @classmethod - def as_line( cls, ncells=None, comm=None ): - domain = Line(name='Omega') - mappings = {'Omega': None} + def from_topological_domain(cls, domain, ncells, comm=None): - geo = Geometry(domain=domain, mappings=mappings, comm=comm) - setattr(geo, 'ncells', ncells) - return geo - - @classmethod - def as_square( cls, ncells=None, comm=None ): - domain = Square(name='Omega') - mappings = {'Omega': None} + if not isinstance(domain, (Line, Square, Cube)): + msg = "Topological domain must be Line, Square or Cube;"\ + " got {} instead.".format(type(domain)) + raise TypeError(msg) + mappings = {domain.interior.name: None} geo = Geometry(domain=domain, mappings=mappings, comm=comm) - setattr(geo, 'ncells', ncells) - return geo + geo.ncells = ncells - @classmethod - def as_cube( cls, ncells=None, comm=None ): - domain = Cube(name='Omega') - mappings = {'Omega': None} - - geo = Geometry(domain=domain, mappings=mappings, comm=comm) - setattr(geo, 'ncells', ncells) return geo + #-------------------------------------------------------------------------- @property def ldim(self): return self._ldim From 6665983ac6bab319e65018226fc3e804f7dccea7 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Mon, 13 Jan 2020 17:01:10 +0100 Subject: [PATCH 21/23] Do not use Sympy version 1.5 --- setup.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/setup.py b/setup.py index 701b70881..2e47d3edb 100644 --- a/setup.py +++ b/setup.py @@ -48,7 +48,7 @@ # Third-party packages from PyPi 'numpy>=1.16', 'scipy>=0.18', - 'sympy>=1.2', + 'sympy>=1.2,<1.5', 'matplotlib', 'pytest', 'pyyaml', From 978cdf98912a7ce7cf6d42bb3f617c3e80cb06ba Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 14 Jan 2020 11:06:26 +0100 Subject: [PATCH 22/23] ci: add Bandit configuration file --- bandit.yml | 7 +++++++ 1 file changed, 7 insertions(+) create mode 100644 bandit.yml diff --git a/bandit.yml b/bandit.yml new file mode 100644 index 000000000..fec715d85 --- /dev/null +++ b/bandit.yml @@ -0,0 +1,7 @@ +# Do not check for use of 'assert' statements (which are standard in unit tests) +# See https://bandit.readthedocs.io/en/latest/plugins/b101_assert_used.html +# +# NOTE: ideally, we would like to only skip this check in our unit tests, but +# we do not know if this is possible. +skips: + - B101 # Ignore assert statements From d4bfd730019ed2590459d52028351bfc8575fe4f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Yaman=20G=C3=BC=C3=A7l=C3=BC?= Date: Tue, 14 Jan 2020 11:09:44 +0100 Subject: [PATCH 23/23] chore: remove unused variables and imports --- psydac/api/ast/utilities.py | 7 ++----- psydac/api/discretization.py | 2 +- psydac/api/grid.py | 1 - 3 files changed, 3 insertions(+), 7 deletions(-) diff --git a/psydac/api/ast/utilities.py b/psydac/api/ast/utilities.py index 442f67a92..32bc149d6 100644 --- a/psydac/api/ast/utilities.py +++ b/psydac/api/ast/utilities.py @@ -674,10 +674,7 @@ def fusion_loops(loops): def compute_boundary_jacobian(parent_namespace, boundary, mapping=None): # Sanity check on arguments - if isinstance(boundary, Boundary): - axis = boundary.axis - ext = boundary.ext - else: + if not isinstance(boundary, Boundary): raise TypeError(boundary) if mapping is None: @@ -686,7 +683,7 @@ def compute_boundary_jacobian(parent_namespace, boundary, mapping=None): else: # Compute metric determinant g on manifold J = SymbolicExpr(mapping.jacobian) - Jm = J[:, [i for i in range(J.shape[1]) if i != axis]] + Jm = J[:, [i for i in range(J.shape[1]) if i != boundary.axis]] g = (Jm.T * Jm).det() # Create statements for computing sqrt(g) diff --git a/psydac/api/discretization.py b/psydac/api/discretization.py index 33b26f2f8..44ee1d18a 100644 --- a/psydac/api/discretization.py +++ b/psydac/api/discretization.py @@ -236,7 +236,7 @@ def solve(self, **kwargs): from sympde.expr import integral from sympde.expr import find - from sympde.topology import element_of, ScalarTestFunction + from sympde.topology import element_of #, ScalarTestFunction # Extract trial functions from model equation u = self.expr.trial_functions diff --git a/psydac/api/grid.py b/psydac/api/grid.py index db41e8d5a..661849bc3 100644 --- a/psydac/api/grid.py +++ b/psydac/api/grid.py @@ -257,7 +257,6 @@ def __init__( self, V, grid, nderiv ): # Modify data for boundary grid if isinstance(grid, BoundaryQuadratureGrid): axis = grid.axis - ext = grid.ext if isinstance(V, ProductFemSpace): for i in range(len(V.spaces)): sp_space = V.spaces[i].spaces[axis]