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 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/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..6752cb3c5 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 @@ -64,15 +65,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 +90,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 +127,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): @@ -205,9 +209,17 @@ 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, - 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 +256,14 @@ 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): + 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 +279,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 +525,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 +564,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 +610,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 @@ -631,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] @@ -706,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') + # ... # ... @@ -745,6 +752,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 @@ -754,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, @@ -769,174 +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.discrete_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.discrete_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: - inv_jac = Symbol('inv_jac') - det_jac = Symbol('det_jac') - - if not isinstance(self.target, Boundary): - - # ... inv jacobian - jac = mapping.det_jacobian - jac = SymbolicExpr(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)] - - init_map = OrderedDict(sorted(init_map.items())) - for stmt in list(init_map.values()): - body += [stmt.subs(1/jac, inv_jac)] + 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)] - else: - body += [Assign(coordinates[i],positions[i][indices_quad[i]]) - for i in range(dim)] - # ... - - # ... - weighted_vol = filter_product(indices_quad, weighted_vols, self.discrete_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])) - # ... + # Convert logical derivatives to physical derivatives + body += [stmt.subs(1/jac, inv_jac) for stmt in self._map_stmts_fields.values()] 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)] - # ... - + 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) + + # Multiply by the correct metric determinant # TODO use positive mapping all the time? Abs? if mapping: - weighted_vol = weighted_vol * Abs(det_jac) - - body.append(Assign(wvol,weighted_vol)) + if isinstance(self.target, Boundary): + weighted_vol *= Abs(det_jac_bnd) + else: + 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 - init_stmts += init_loop_quadrature( indices_quad, ranges_quad, - self.discrete_boundary ) + # If we are assemblying on a boundary, some indices are set to zero + init_stmts += init_loop_quadrature( indices_quad, ranges_quad, self.boundary ) - body = select_loops( indices_quad, ranges_quad, body, - self.discrete_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)]) @@ -947,32 +934,32 @@ 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.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) # ... @@ -997,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 @@ -1016,8 +1001,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 @@ -1029,18 +1013,42 @@ 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 - math_elements = math_atoms_as_str(self.kernel_expr, 'numpy') + #------------------------------------------------------------------ + # 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.update(new) + elif isinstance(i, (Assign, AugAssign)): + new = math_atoms_as_str(i.rhs, lib) + math_elements.update(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] + self._imports.update(math_imports) + + #------------------------------------------------------------------ + # Identify function arguments + #------------------------------------------------------------------ - imports += math_imports - self._imports = 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 - func_args = self.build_arguments(func_args) + #------------------------------------------------------------------ + # Add header and decorators, if needed + #------------------------------------------------------------------ + decorators = {} header = None if self.backend['name'] == 'pyccel': @@ -1051,13 +1059,22 @@ def _initialize(self): 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): @@ -1078,11 +1095,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 +1180,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 +1427,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..32bc149d6 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,85 +671,68 @@ def fusion_loops(loops): return loops_cp #============================================================================== -def compute_normal_vector(vector, discrete_boundary, mapping): - dim = len(vector) - pdim = dim - len(discrete_boundary) - if len(discrete_boundary) > 1: raise NotImplementedError('TODO') +def compute_boundary_jacobian(parent_namespace, boundary, mapping=None): - face = discrete_boundary[0] - axis = face[0] ; ext = face[1] + # Sanity check on arguments + if not isinstance(boundary, Boundary): + raise TypeError(boundary) - map_stmts = [] - body = [] - - if not mapping: - - values = np.zeros(dim) - values[axis] = ext + if mapping is None: + stmts = [] else: - M = mapping - inv_jac = Symbol('inv_jac') - det_jac = Symbol('det_jac') - - # ... 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) + # Compute metric determinant g on manifold + J = SymbolicExpr(mapping.jacobian) + Jm = J[:, [i for i in range(J.shape[1]) if i != boundary.axis]] + g = (Jm.T * Jm).det() - 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*J[1], -inv_jac*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*v for v in values] + # Create statements for computing sqrt(g) + det_jac_bnd = parent_namespace['det_jac_bnd'] + stmts = [Assign(det_jac_bnd, sympy_sqrt(g))] + return stmts - # change the orientation - values = [ext*i for i in values] +#============================================================================== +def compute_normal_vector(parent_namespace, vector, boundary, mapping=None): - map_stmts += [Assign(det_jac, j)] - map_stmts += [Assign(inv_jac, 1./j)] + # Sanity check on arguments + if isinstance(boundary, Boundary): + axis = boundary.axis + ext = boundary.ext + 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: + 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])] - for i in range(0, dim): - body += [Assign(vector[i], values[i])] + # Create statements for computing normal vector components + stmts = [Assign(lhs, rhs) for lhs, rhs in zip(vector, values)] -# print(map_stmts) -# print(body) -# import sys; sys.exit(0) + # 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 map_stmts, body + return stmts #============================================================================== -def compute_tangent_vector(vector, discrete_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])') 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/discretization.py b/psydac/api/discretization.py index fc110c5a8..44ee1d18a 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) # ... @@ -138,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): @@ -154,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 @@ -194,7 +200,10 @@ def assemble(self, **kwargs): self._linear_system = LinearSystem(M, rhs) def solve(self, **kwargs): - settings = kwargs.pop('settings', _default_solver) + + settings = _default_solver.copy() + settings.update(kwargs.pop('settings', {})) + rhs = kwargs.pop('rhs', None) if rhs: kwargs['assemble_rhs'] = False @@ -206,6 +215,68 @@ 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)) + + # 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 + 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 + domain_h = self.domain + Vh = self.trial_space + equation_h = discretize(equation, domain_h, [Vh, Vh]) + + # 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 + + #---------------------------------------------------------------------- + return driver_solve(self.linear_system, **settings) #============================================================================== @@ -320,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) @@ -353,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)] @@ -424,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) - - if not( ncells is None ): - dtype = domain.dtype - - if dtype['type'].lower() == 'line' : - return Geometry.as_line(ncells, comm=comm) - - 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']) +def discretize_domain(domain, *, filename=None, ncells=None, comm=None): - raise ValueError(msg) + if not (filename or ncells): + raise ValueError("Must provide either 'filename' or 'ncells'") - elif not( filename is None ): - geometry = Geometry(filename=filename, comm=comm) + elif filename and ncells: + raise ValueError("Cannot provide both 'filename' and 'ncells'") - 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/api/fem.py b/psydac/api/fem.py index e51ab217c..0fb384c0e 100644 --- a/psydac/api/fem.py +++ b/psydac/api/fem.py @@ -73,11 +73,17 @@ 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) # ... - test_space = self.spaces[0] - trial_space = self.spaces[1] + trial_space = self.spaces[0] + test_space = self.spaces[1] # ... # ... @@ -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..661849bc3 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,22 @@ 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 + 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 diff --git a/psydac/api/tests/test_api_2d_scalar.py b/psydac/api/tests/test_api_2d_scalar.py index 01ea0e0bb..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 @@ -20,365 +21,268 @@ 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() - # ... + if not args: + return () + else: + assert all(1 <= a <= 4 for a in args) + assert len(set(args)) == len(args) - # ... - phi = FemField( Vh, x ) - # ... + boundaries = {1: {'axis': 0, 'ext': -1}, + 2: {'axis': 0, 'ext': 1}, + 3: {'axis': 1, 'ext': -1}, + 4: {'axis': 1, 'ext': 1}} - # ... 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(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 + #+++++++++++++++++++++++++++++++ + # 1. 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: 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, 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 #============================================================================== -def run_laplace_2d_neu(solution, f, ncells, degree, comm=None): - - # ... abstract model - domain = Square() - - V = ScalarFunctionSpace('V', domain) +def run_laplace_2d(solution, f, dir_zero_boundary, dir_nonzero_boundary, + ncells, degree, comm=None): - B_neumann = domain.boundary + assert isinstance(dir_zero_boundary , (list, tuple)) + assert isinstance(dir_nonzero_boundary, (list, tuple)) - x,y = domain.coordinates - - 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) + # Bilinear form a: V x V --> R + a = BilinearForm((u, v), integral(domain, dot(grad(u), grad(v)) + u * v)) - expr = dot(grad(v), grad(u)) + v*u - a = BilinearForm((v,u), int_0(expr)) - - 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 #============================================================================== -def run_biharmonic_2d_dir(solution, f, ncells, degree, comm=None): - - # ... 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 = laplace(v) * laplace(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') - - 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) - # ... +def run_biharmonic_2d_dir(solution, f, dir_zero_boundary, ncells, degree, comm=None): - # ... 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() - # ... + assert isinstance(dir_zero_boundary, (list, tuple)) - # ... - 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 - - -#============================================================================== -def run_poisson_user_function_2d_dir(f, solution, ncells, degree, comm=None): - - # ... abstract model + #+++++++++++++++++++++++++++++++ + # 1. Abstract model + #+++++++++++++++++++++++++++++++ domain = Square() - x,y = domain.coordinates - f = implemented_function('f', f) + B_dirichlet_0 = Union(*[domain.get_boundary(**kw) for kw in dir_zero_boundary]) + B_dirichlet_i = domain.boundary.complement(B_dirichlet_0) - V = ScalarFunctionSpace('V', domain) - - 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 = dot(grad(v), grad(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(x,y)*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') - bc = EssentialBC(u, 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() - # ... - - # ... - phi = FemField( Vh, x ) - # ... + #+++++++++++++++++++++++++++++++ + # 3. Solution + #+++++++++++++++++++++++++++++++ - # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) - # ... + # Solve linear system + x = equation_h.solve() + uh = FemField( Vh, x ) - return l2_error, h1_error + # 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 ############################################################################### # SERIAL TESTS ############################################################################### #============================================================================== -def test_api_poisson_2d_dir_1(): - - 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) - 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 @@ -386,16 +290,17 @@ def test_api_poisson_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_2d_dirneu_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 - 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 @@ -403,16 +308,17 @@ def test_api_poisson_2d_dirneu_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_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 - 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 @@ -420,16 +326,17 @@ def test_api_poisson_2d_dirneu_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_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 - 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 @@ -437,16 +344,17 @@ def test_api_poisson_2d_dirneu_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_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 - 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 @@ -454,18 +362,17 @@ def test_api_poisson_2d_dirneu_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_dirneu_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 - 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 @@ -473,18 +380,17 @@ def test_api_poisson_2d_dirneu_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_dirneu_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 - 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 @@ -492,19 +398,17 @@ def test_api_poisson_2d_dirneu_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_dirneu_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 - 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 @@ -512,65 +416,135 @@ 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_poisson_2d_dir0_234_neui_1(): -#============================================================================== -def test_api_laplace_2d_neu(): + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution - from sympy.abc import x,y + dir_zero_boundary = get_boundaries(2, 3, 4) + dir_nonzero_boundary = get_boundaries() - solution = cos(pi*x)*cos(pi*y) - f = (2.*pi**2 + 1.)*solution + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) - l2_error, h1_error = run_laplace_2d_neu(solution, f, ncells=[2**3,2**3], degree=[2,2]) + expected_l2_error = 0.00021786960672322118 + expected_h1_error = 0.01302350067761091 - expected_l2_error = 0.0002172846538950129 - expected_h1_error = 0.012984852988125026 + 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_134_neui_2(): + + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution + + 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 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_1(): +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_124_neui_3(): + + solution = sin(pi*x)*sin(pi*y) + f = 2*pi**2*solution - from sympy.abc import x,y - from sympde.expr import TerminalExpr + dir_zero_boundary = get_boundaries(1, 2, 4) + dir_nonzero_boundary = get_boundaries() - solution = (sin(pi*x)*sin(pi*y))**2 + l2_error, h1_error = run_poisson_2d(solution, f, dir_zero_boundary, + dir_nonzero_boundary, ncells=[2**3, 2**3], degree=[2, 2]) - # compute the analytical solution - f = laplace(laplace(solution)) - f = TerminalExpr(f, dim=2) + expected_l2_error = 0.00021786960672322118 + expected_h1_error = 0.01302350067761091 - l2_error, h1_error = run_biharmonic_2d_dir(solution, f, - ncells=[2**3,2**3], degree=[2,2]) + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) - expected_l2_error = 0.015086415626061608 - expected_h1_error = 0.08773346232942228 +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_123_neui_4(): + + 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() + + 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 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_123_diri_4(): + solution = sin(pi * x) * sin(0.5*pi * y) + f = 5/4*pi**2 * solution -#============================================================================== -def test_api_poisson_user_function_2d_dir_1(): + 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_poisson_2d_dir0_13_diri_24(): + + 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 - from sympy.abc import x,y +#------------------------------------------------------------------------------ +def test_poisson_2d_dir0_1234_user_function(): solution = sin(pi*x)*sin(pi*y) # ... - def f(x,y): - from numpy import pi - from numpy import cos - from numpy import sin + # 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) - value = 2*pi**2*sin(pi*x)*sin(pi*y) - return value + # Python function is converted to Sympy's "implemented function" and then + # called with symbolic arguments (x, y): + f = implemented_function('f', f)(x, y) # ... - l2_error, h1_error = run_poisson_user_function_2d_dir(f, solution, - 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 @@ -578,22 +552,104 @@ def f(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 + + 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 + + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) + +#============================================================================== +# 2D biharmonic equation +#============================================================================== +def test_biharmonic_2d_dir0_1234(): + + solution = sin(pi * x)**2 * sin(pi * y)**2 + f = laplace(laplace(solution)) + + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + + 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.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_biharmonic_2d_dir0_123_diri_4(): + + 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_biharmonic_2d_dir0_13_diri_24(): + + 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 + ############################################################################### # PARALLEL TESTS ############################################################################### #============================================================================== @pytest.mark.parallel -def test_api_poisson_2d_dir_1_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) - 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 diff --git a/psydac/api/tests/test_api_2d_scalar_mapping.py b/psydac/api/tests/test_api_2d_scalar_mapping.py index 2706f7959..d4fa4a753 100644 --- a/psydac/api/tests/test_api_2d_scalar_mapping.py +++ b/psydac/api/tests/test_api_2d_scalar_mapping.py @@ -1,7 +1,22 @@ # -*- 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 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 +from sympy.abc import x, y import pytest import os @@ -30,512 +45,574 @@ # ... #============================================================================== -def run_poisson_2d_dir(filename, solution, f, dir_boundary=None, comm=None): +def get_boundaries(*args): - # ... abstract model - domain = Domain.from_file(filename) + if not args: + return () + else: + assert all(1 <= a <= 4 for a in args) + assert len(set(args)) == len(args) - V = ScalarFunctionSpace('V', domain) + boundaries = {1: {'axis': 0, 'ext': -1}, + 2: {'axis': 0, 'ext': 1}, + 3: {'axis': 1, 'ext': -1}, + 4: {'axis': 1, 'ext': 1}} - if dir_boundary is None: - B_dirichlet = domain.boundary - elif len(dir_boundary) == 1: - B_dirichlet = domain.get_boundary(**dir_boundary[0]) - else: - B_dirichlet = Union(*[domain.get_boundary(**kw) for kw in dir_boundary]) + return tuple(boundaries[i] for i in args) - x,y = domain.coordinates +#============================================================================== +def run_poisson_2d(filename, solution, f, dir_zero_boundary, + dir_nonzero_boundary, neumann_boundary, comm=None): - F = element_of(V, name='F') + assert isinstance( dir_zero_boundary, (list, tuple)) + assert isinstance(dir_nonzero_boundary, (list, tuple)) + assert isinstance( neumann_boundary, (list, tuple)) - v = element_of(V, name='v') - u = element_of(V, name='u') + #+++++++++++++++++++++++++++++++ + # 1. Abstract model + #+++++++++++++++++++++++++++++++ + domain = Domain.from_file(filename) - int_0 = lambda expr: integral(domain , expr) + 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]) - expr = dot(grad(v), grad(u)) - a = BilinearForm((v,u), int_0(expr)) + V = ScalarFunctionSpace('V', domain) + u = element_of(V, name='u') + v = element_of(V, name='v') + nn = NormalVector('nn') - expr = f*v - l = LinearForm(v, int_0(expr)) + # Bilinear form a: V x V --> R + a = BilinearForm((u, v), integral(domain, dot(grad(u), grad(v)))) - error = F - solution + # 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 + + # 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)] + + # 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') - bc = EssentialBC(u, 0, B_dirichlet) - 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) - # ... - # ... 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_poisson_2d_dirneu(filename, solution, f, boundary, comm=None): +def run_laplace_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)) + 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') - 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) +def run_biharmonic_2d_dir(filename, solution, f, dir_zero_boundary, + dir_nonzero_boundary, comm=None): - V = ScalarFunctionSpace('V', domain) + assert isinstance( dir_zero_boundary, (list, tuple)) + assert isinstance(dir_nonzero_boundary, (list, tuple)) - B_neumann = domain.boundary - - x,y = domain.coordinates - - 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]) + 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) + # Bilinear form a: V x V --> R + a = BilinearForm((u, v), integral(domain, laplace(u) * laplace(v))) - expr = dot(grad(v), grad(u)) + v*u - a = BilinearForm((v,u), int_0(expr)) + # Linear form l: V --> R + l = LinearForm(v, integral(domain, f * v)) - expr = f*v - l0 = 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)] - expr = v*dot(grad(solution), nn) - l_B_neumann = LinearForm(v, int_1(expr)) + # Variational model + equation = find(u, forall=v, lhs=a(u, v), rhs=l(v), bc=bc) - expr = l0(v) + l_B_neumann(v) - l = LinearForm(v, expr) - - 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') - 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) - # ... + 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 +############################################################################### +# SERIAL TESTS +############################################################################### + +#============================================================================== +# 2D Poisson's equation with identity map #============================================================================== -def run_biharmonic_2d_dir(filename, solution, f, comm=None): +def test_poisson_2d_identity_dir0_1234(): - # ... abstract model - domain = Domain.from_file(filename) + 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) - V = ScalarFunctionSpace('V', domain) + dir_zero_boundary = get_boundaries(1, 2, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries() - F = element_of(V, name='F') + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - v = element_of(V, name='v') - u = element_of(V, name='u') + expected_l2_error = 0.00021808678604159413 + expected_h1_error = 0.013023570720357957 - int_0 = lambda expr: integral(domain , expr) + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) - expr = laplace(v) * laplace(u) - a = BilinearForm((v,u),int_0(expr)) +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_234_neu0_1(): - expr = f*v - l = LinearForm(v, int_0(expr)) + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = cos(0.5*pi*x)*sin(pi*y) + f = (5./4.)*pi**2*solution - error = F - solution - l2norm = Norm(error, domain, kind='l2') - h1norm = Norm(error, domain, kind='h1') - h2norm = Norm(error, domain, kind='h2') + dir_zero_boundary = get_boundaries(2, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1) - 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) - # ... + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - # ... create the computational domain from a topological domain - domain_h = discretize(domain, filename=filename, comm=comm) - # ... + expected_l2_error = 0.00015546057795986509 + expected_h1_error = 0.009269302784527035 - # ... discrete spaces - Vh = discretize(V, domain_h) - # ... + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) - # ... dsicretize the equation using Dirichlet bc - equation_h = discretize(equation, domain_h, [Vh, Vh]) - # ... +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_134_neu0_2(): - # ... discretize norms - l2norm_h = discretize(l2norm, domain_h, Vh) - h1norm_h = discretize(h1norm, domain_h, Vh) - h2norm_h = discretize(h2norm, domain_h, Vh) - # ... + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(0.5*pi*x)*sin(pi*y) + f = (5./4.)*pi**2*solution - # ... solve the discrete equation - x = equation_h.solve() - # ... + dir_zero_boundary = get_boundaries(1, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(2) - # ... - phi = FemField( Vh, x ) - # ... + l2_error, h1_error = run_poisson_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) - # ... compute norms - l2_error = l2norm_h.assemble(F=phi) - h1_error = h1norm_h.assemble(F=phi) - h2_error = h2norm_h.assemble(F=phi) - # ... + expected_l2_error = 0.00015546057795095866 + expected_h1_error = 0.009269302784528054 - return l2_error, h1_error, h2_error + assert( abs(l2_error - expected_l2_error) < 1.e-7) + assert( abs(h1_error - expected_h1_error) < 1.e-7) -############################################################################### -# SERIAL TESTS -############################################################################### +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_124_neu0_3(): -#============================================================================== -def test_api_poisson_2d_dir_identity(): filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = sin(pi*x)*cos(0.5*pi*y) + f = (5./4.)*pi**2*solution - from sympy.abc import x,y - - solution = sin(pi*x)*sin(pi*y) - f = 2*pi**2*sin(pi*x)*sin(pi*y) + 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_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.00021808678604159413 - expected_h1_error = 0.013023570720357957 + 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_dir_collela(): - filename = os.path.join(mesh_dir, 'collela_2d.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_123_neu0_4(): - from sympy.abc import x,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 - 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) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(4) - 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.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_dir_quart_circle(): - filename = os.path.join(mesh_dir, 'quart_circle.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 - 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(2, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1, 3) - 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 = 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_dir_circle(): - filename = os.path.join(mesh_dir, 'circle.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_4_neu0_123(): - from sympy.abc import x,y + filename = os.path.join(mesh_dir, 'identity_2d.h5') + solution = cos(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(4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1, 2, 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.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_1(): - 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 = cos(0.5*pi*x)*sin(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': 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.00015546057795986509 - expected_h1_error = 0.009269302784527035 + 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_2(): +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_134_neui_2(): + 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, 3, 4) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(2) - solution = sin(0.5*pi*x)*sin(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}]) - - expected_l2_error = 0.00015546057795095866 - expected_h1_error = 0.009269302784528054 + 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_3(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_124_neui_3(): - 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)*cos(0.5*pi*y) - f = (5./4.)*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_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.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_4(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') +#------------------------------------------------------------------------------ +def test_poisson_2d_identity_dir0_123_neui_4(): - 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(1, 2, 3) + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(4) - 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(): +#------------------------------------------------------------------------------ +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 - from sympy.abc import x,y + dir_zero_boundary = get_boundaries(1, 2, 3) + dir_nonzero_boundary = get_boundaries(4) + neumann_boundary = get_boundaries() - solution = cos(0.5*pi*x)*cos(0.5*pi*y) - f = (1./2.)*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': 1, 'ext': -1}]) + expected_l2_error = 0.0001529221571156830 + expected_h1_error = 0.009293161646612863 - 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 - 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(): -#============================================================================== -def test_api_poisson_2d_dirneu_identity_123(): 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 - from sympy.abc import x,y + dir_zero_boundary = get_boundaries(1, 3) + dir_nonzero_boundary = get_boundaries(2, 4) + neumann_boundary = get_boundaries() - 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.0007786454571731944 + expected_h1_error = 0.0449669071240554 - 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 + +#============================================================================== +# 2D Poisson's equation with "Collela" map +#============================================================================== +def test_poisson_2d_collela_dir0_1234(): + + 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) +#------------------------------------------------------------------------------ +def test_poisson_2d_collela_dir0_234_neu0_1(): -##============================================================================== -# 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_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 + + 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 + + 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_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 @@ -543,34 +620,39 @@ 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) -##============================================================================== -## 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_poisson_2d_collela_dir0_124_neu0_3(): -#============================================================================== -def test_api_poisson_2d_dirneu_collela_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 + + 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 + + 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_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 @@ -578,16 +660,287 @@ 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_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.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) + +#------------------------------------------------------------------------------ +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.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) + +#------------------------------------------------------------------------------ +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.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) + +#------------------------------------------------------------------------------ +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.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) + +#------------------------------------------------------------------------------ +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 + #============================================================================== -def test_api_laplace_2d_neu_identity(): - filename = os.path.join(mesh_dir, 'identity_2d.h5') +# 2D Poisson's equation on quarter annulus +#============================================================================== +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) + + 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.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_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.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 + +#------------------------------------------------------------------------------ +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) + + expected_l2_error = 0.0006663772402662598 + expected_h1_error = 0.025906309642232804 + + 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 +#============================================================================== +def test_poisson_2d_circle_dir0(): + + 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) + + 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(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) + + 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) + +#============================================================================== +# 2D "Laplace-like" equation +#============================================================================== +def test_laplace_2d_identity_neu0_1234(): + + filename = os.path.join(mesh_dir, 'identity_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) + + l2_error, h1_error = run_laplace_2d(filename, solution, f, + dir_zero_boundary, dir_nonzero_boundary, neumann_boundary) expected_l2_error = 0.00021728465388208586 expected_h1_error = 0.012984852988123631 @@ -595,44 +948,93 @@ def test_api_laplace_2d_neu_identity(): 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, 'quart_circle.h5') +#------------------------------------------------------------------------------ +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 + + dir_zero_boundary = get_boundaries() + dir_nonzero_boundary = get_boundaries() + neumann_boundary = get_boundaries(1, 2, 3, 4) + + 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) - 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() - expected_l2_error = 0.3092792236008727 - expected_h1_error = 1.3320441589030887 - expected_h2_error = 6.826223014197719 + 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 + 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_biharmonic_2d_collela_dir0_1234(): -##============================================================================== -## 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) + filename = os.path.join(mesh_dir, 'collela_2d.h5') + solution = (cos(pi*x/2)*cos(pi*y/2))**2 + f = laplace(laplace(solution)) + + 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 + 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_biharmonic_2d_quarter_annulus_dir0_1234(): + + 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)) + + 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 + 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) ############################################################################### # PARALLEL TESTS @@ -640,16 +1042,19 @@ def test_api_biharmonic_2d_dir_identity(): #============================================================================== @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 @@ -657,7 +1062,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 #============================================================================== 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() 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 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',