diff --git a/fenicsxprecice/adapter_core.py b/fenicsxprecice/adapter_core.py index 17db94f..91f93e9 100644 --- a/fenicsxprecice/adapter_core.py +++ b/fenicsxprecice/adapter_core.py @@ -2,7 +2,8 @@ This module consists of helper functions used in the Adapter class. Names of the functions are self explanatory """ -from dolfinx.fem import FunctionSpace, Function +#from dolfinx.fem import FunctionSpace, Function +from dolfinx import fem, geometry import numpy as np from enum import Enum import logging @@ -69,12 +70,12 @@ def determine_function_type(input_obj): tag : bool 0 if input_function is SCALAR and 1 if input_function is VECTOR. """ - if isinstance(input_obj, FunctionSpace): # scalar-valued functions have rank 0 is FEniCSx + if isinstance(input_obj, fem.FunctionSpace): # scalar-valued functions have rank 0 is FEniCSx if input_obj.num_sub_spaces == 0: return FunctionType.SCALAR elif input_obj.num_sub_spaces == 2: return FunctionType.VECTOR - elif isinstance(input_obj, Function): + elif isinstance(input_obj, fem.Function): if len(input_obj.x.array.shape) == 1: return FunctionType.SCALAR elif input_obj.x.array.shape[1] > 1: @@ -85,16 +86,16 @@ def determine_function_type(input_obj): raise Exception("Error determining type of given dolfin FunctionSpace") -def convert_fenicsx_to_precice(fenicsx_function, local_ids): +def convert_fenicsx_to_precice_coordinateBased(fenicsx_function, local_coords): """ - Converts data of type dolfin.Function into Numpy array for all x and y coordinates on the boundary. + Converts data of type dolfinx.Function into Numpy array for all x and y coordinates on the boundary. Parameters ---------- fenicsx_function : FEniCSx function A FEniCSx function referring to a physical variable in the problem. - local_ids: numpy array - Array of local indices of vertices on the coupling interface and owned by this rank. + local_coords: numpy array + Array of local coordinates of vertices on the coupling interface and owned by this rank. Returns ------- @@ -102,32 +103,30 @@ def convert_fenicsx_to_precice(fenicsx_function, local_ids): Array of FEniCSx function values at each point on the boundary. """ - if not isinstance(fenicsx_function, Function): + if not isinstance(fenicsx_function, fem.Function): raise Exception("Cannot handle data type {}".format(type(fenicsx_function))) - precice_data = [] - # sampled_data = fenicsx_function.x.array # that works only for 1st order elements where dofs = grid points - # TODO begin dirty fix. See https://github.com/precice/fenicsx-adapter/issues/17 for details. - x_mesh = fenicsx_function.function_space.mesh.geometry.x - x_dofs = fenicsx_function.function_space.tabulate_dof_coordinates() - mask = [] # where dof coordinate == mesh coordinate - for i in range(x_dofs.shape[0]): - for j in range(x_mesh.shape[0]): - if np.allclose(x_dofs[i, :], x_mesh[j, :], 1e-15): - mask.append(i) - break - sampled_data = fenicsx_function.x.array[mask] - # end dirty fix - - if len(local_ids): - if fenicsx_function.function_space.num_sub_spaces > 0: # function space is VectorFunctionSpace - raise Exception("Functions from VectorFunctionSpaces are currently not supported.") - else: # function space is FunctionSpace (scalar) - for lid in local_ids: - precice_data.append(sampled_data[lid]) - else: - precice_data = np.array([]) - + mesh = fenicsx_function.function_space.mesh + + # this evaluation is a bit annoying, see: + # https://github.com/FEniCS/dolfinx/blob/main/python/test/unit/fem/test_function.py#L63 + + # for fast function evaluation + bb_tree = geometry.bb_tree(mesh, mesh.geometry.dim) # TODO: as long as the domain didn't change, we could store that tree somewhere + + cells = [] + points = [] + + # Find cells whose bounding-box collide with the the points + cell_candidates = geometry.compute_collisions_points(bb_tree, local_coords) + # Choose one of the cells that contains the point + colliding_cells = geometry.compute_colliding_cells(mesh, cell_candidates, local_coords) + for i, point in enumerate(local_coords): + if len(colliding_cells.links(i)) > 0: + points.append(point) + cells.append(colliding_cells.links(i)[0]) + + precice_data = fenicsx_function.eval(points, cells) return np.array(precice_data) @@ -156,13 +155,12 @@ def get_fenicsx_vertices(function_space, coupling_subdomain, dims): # Get mesh from FEniCSx function space mesh = function_space.mesh - + # Get coordinates and IDs of all vertices of the mesh which lie on the coupling boundary. try: - on_subdomain = coupling_subdomain(mesh.geometry.x.T) - ids, = np.where(on_subdomain) + ids = fem.locate_dofs_geometrical(function_space, coupling_subdomain) if dims == 2: - coords = mesh.geometry.x[ids][:, :2] + coords = function_space.tabulate_dof_coordinates()[ids] # we get 3d coordinates here else: coords = np.array([]) except Exception as e: # fall back to old method # TODO is that too general? Better use, e.g., IndexError here? @@ -177,4 +175,5 @@ def get_fenicsx_vertices(function_space, coupling_subdomain, dims): coords.append([v[0], v[1]]) ids = np.array(ids) coords = np.array(coords) + return ids, coords diff --git a/fenicsxprecice/fenicsxprecice.py b/fenicsxprecice/fenicsxprecice.py index cf05b38..0b91a28 100644 --- a/fenicsxprecice/fenicsxprecice.py +++ b/fenicsxprecice/fenicsxprecice.py @@ -7,7 +7,7 @@ from .config import Config import logging import precice -from .adapter_core import FunctionType, determine_function_type, convert_fenicsx_to_precice, get_fenicsx_vertices, CouplingMode, Vertices +from .adapter_core import FunctionType, determine_function_type, get_fenicsx_vertices, CouplingMode, Vertices, convert_fenicsx_to_precice_coordinateBased from .expression_core import SegregatedRBFInterpolationExpression from .solverstate import SolverState from dolfinx.fem import Function, FunctionSpace @@ -52,10 +52,14 @@ def __init__(self, mpi_comm, adapter_config_filename='precice-adapter-config.jso # Setup up MPI communicator self._comm = mpi_comm - - self._interface = precice.Interface(self._config.get_participant_name(), self._config.get_config_file_name(), - self._comm.Get_rank(), self._comm.Get_size()) - + + self._participant = precice.Participant( + self._config.get_participant_name(), + self._config.get_config_file_name(), + self._comm.Get_rank(), + self._comm.Get_size() + ) + # FEniCSx related quantities self._read_function_space = None # initialized later self._write_function_space = None # initialized later @@ -83,6 +87,9 @@ def __init__(self, mpi_comm, adapter_config_filename='precice-adapter-config.jso # Problem dimension in FEniCSx self._fenicsx_dims = None + + + self._empty_rank = True def create_coupling_expression(self): """ @@ -112,7 +119,7 @@ def update_coupling_expression(self, coupling_expression, data): coupling_expression : Object of class dolfinx.functions.expression.Expression Reference to object of class GeneralInterpolationExpression or ExactInterpolationExpression. data : dict_like - The coupling data. A dictionary containing nodal data with vertex coordinates as key and associated data as + The coupling data. A dictionary containing the values of the vertex coordinates as key and associated data as value. """ vertices = np.array(list(data.keys())) @@ -122,7 +129,7 @@ def update_coupling_expression(self, coupling_expression, data): def get_point_sources(self, data): raise Exception("PointSources are not implemented for the FEniCSx adapter.") - def read_data(self): + def read_data(self, dt): """ Read data from preCICE. Data is generated depending on the type of the read function (Scalar or Vector). For a scalar read function the data is a numpy array with shape (N) where N = number of coupling vertices @@ -141,18 +148,21 @@ def read_data(self): assert (self._coupling_type is CouplingMode.UNI_DIRECTIONAL_READ_COUPLING or CouplingMode.BI_DIRECTIONAL_COUPLING) - read_data_id = self._interface.get_data_id(self._config.get_read_data_name(), - self._interface.get_mesh_id(self._config.get_coupling_mesh_name())) - - read_data = None - - if self._read_function_type is FunctionType.SCALAR: - read_data = self._interface.read_block_scalar_data(read_data_id, self._precice_vertex_ids) - elif self._read_function_type is FunctionType.VECTOR: - read_data = self._interface.read_block_vector_data(read_data_id, self._precice_vertex_ids) - - read_data = {tuple(key): value for key, value in zip(self._fenicsx_vertices.get_coordinates(), read_data)} + read_data=None + + if not self._empty_rank: + read_data = self._participant.read_data( + self._config.get_coupling_mesh_name(), + self._config.get_read_data_name(), + self._precice_vertex_ids, + dt + ) + #TODO: MPI stuff + read_data = {tuple(key): value for key, value in zip(self._fenicsx_vertices.get_coordinates(), read_data)} + else: + pass + return copy.deepcopy(read_data) def write_data(self, write_function): @@ -169,28 +179,25 @@ def write_data(self, write_function): assert (self._coupling_type is CouplingMode.UNI_DIRECTIONAL_WRITE_COUPLING or CouplingMode.BI_DIRECTIONAL_COUPLING) - w_func = write_function + w_func = write_function.copy() # Check that the function provided lives on the same function space provided during initialization assert (self._write_function_type == determine_function_type(w_func)) - # TODO this raises AssertionError, not sure why. I just commented it out, still works... - # assert (write_function.function_space == self._write_function_space) + assert (write_function.function_space == self._write_function_space) - write_data_id = self._interface.get_data_id(self._config.get_write_data_name(), - self._interface.get_mesh_id(self._config.get_coupling_mesh_name())) + #write_data_id = self._participant.get_data_id(self._config.get_write_data_name(), + # self._participant.get_mesh_id(self._config.get_coupling_mesh_name())) + write_function_type = determine_function_type(write_function) assert (write_function_type in list(FunctionType)) - write_data = convert_fenicsx_to_precice(write_function, self._fenicsx_vertices.get_ids()) - if write_function_type is FunctionType.SCALAR: - assert (write_function.function_space.num_sub_spaces == 0) - write_data = np.squeeze(write_data) # TODO dirty solution - self._interface.write_block_scalar_data(write_data_id, self._precice_vertex_ids, write_data) - elif write_function_type is FunctionType.VECTOR: - assert (write_function.function_space.num_sub_spaces > 0) - self._interface.write_block_vector_data(write_data_id, self._precice_vertex_ids, write_data) - else: - raise Exception("write_function provided is neither VECTOR nor SCALAR type") + write_data = convert_fenicsx_to_precice_coordinateBased(write_function, self._fenicsx_vertices.get_coordinates()) + self._participant.write_data( + self._config.get_coupling_mesh_name(), + self._config.get_write_data_name(), + self._precice_vertex_ids, + write_data + ) def initialize(self, coupling_subdomain, read_function_space=None, write_object=None): """ @@ -268,8 +275,26 @@ def initialize(self, coupling_subdomain, read_function_space=None, write_object= # Ensure that function spaces of read and write functions are defined using the same mesh self._write_function_type = determine_function_type(write_function_space) self._write_function_space = write_function_space - + + + # Set vertices on the coupling subdomain for this rank self._fenicsx_dims = function_space.mesh.geometry.dim + ids, coords = get_fenicsx_vertices(function_space, coupling_subdomain, self._fenicsx_dims) # returns 3d coordinates (necessary later for writing the data!) + # this isnt a problem in update_coupling_expression, because in this function + # , the two first dimensions are extracted. Exactly what we want! + self._fenicsx_vertices.set_ids(ids) + self._fenicsx_vertices.set_coordinates(coords) + + # Set up mesh in preCICE + self._precice_vertex_ids = self._participant.set_mesh_vertices( + self._config.get_coupling_mesh_name(), self._fenicsx_vertices.get_coordinates()[:, :2]) # give preCICE only 2D coordinates + + + if self._fenicsx_vertices.get_ids().size > 0: + self._empty_rank = False + else: + print("Rank {} has no part of coupling boundary.".format(self._comm.Get_rank())) + # Ensure that function spaces of read and write functions use the same mesh if self._coupling_type is CouplingMode.BI_DIRECTIONAL_COUPLING: @@ -279,37 +304,23 @@ def initialize(self, coupling_subdomain, read_function_space=None, write_object= if self._fenicsx_dims != 2: raise Exception("Currently the fenicsx-adapter only supports 2D cases") - if self._fenicsx_dims != self._interface.get_dimensions(): + if self._fenicsx_dims != self._participant.get_mesh_dimensions(self._config.get_coupling_mesh_name()): raise Exception("Dimension of preCICE setup and FEniCSx do not match") - # Set vertices on the coupling subdomain for this rank - ids, coords = get_fenicsx_vertices(function_space, coupling_subdomain, self._fenicsx_dims) - self._fenicsx_vertices.set_ids(ids) - self._fenicsx_vertices.set_coordinates(coords) - - # Set up mesh in preCICE - self._precice_vertex_ids = self._interface.set_mesh_vertices(self._interface.get_mesh_id( - self._config.get_coupling_mesh_name()), self._fenicsx_vertices.get_coordinates()) - - precice_dt = self._interface.initialize() - - if self._interface.is_action_required(precice.action_write_initial_data()): + if self._participant.requires_initial_data(): if not write_function: - raise Exception("Non-standard initialization requires a write_function") + raise Exception("preCICE requires you to write initial data. Please provide a write_function to initialize(...)") self.write_data(write_function) - self._interface.mark_action_fulfilled(precice.action_write_initial_data()) - self._interface.initialize_data() + self._participant.initialize() - return precice_dt - - def store_checkpoint(self, user_u, t, n): + def store_checkpoint(self, payload, t, n): """ Defines an object of class SolverState which stores the current state of the variable and the time stamp. Parameters ---------- - user_u : FEniCSx Function + payload : FEniCSx Function Current state of the physical variable of interest for this participant. t : double Current simulation time. @@ -320,11 +331,10 @@ def store_checkpoint(self, user_u, t, n): assert (self.is_time_window_complete()) logger.debug("Store checkpoint") - my_u = user_u.copy() + my_u = payload.copy() # making sure that the FEniCSx function provided by user is not directly accessed by the Adapter - assert (my_u != user_u) + assert (my_u != payload) self._checkpoint = SolverState(my_u, t, n) - self._interface.mark_action_fulfilled(self.action_write_iteration_checkpoint()) def retrieve_checkpoint(self): """ @@ -341,7 +351,6 @@ def retrieve_checkpoint(self): """ assert (not self.is_time_window_complete()) logger.debug("Restore solver state") - self._interface.mark_action_fulfilled(self.action_read_iteration_checkpoint()) return self._checkpoint.get_state() def advance(self, dt): @@ -363,7 +372,7 @@ def advance(self, dt): Maximum length of timestep to be computed by solver. """ self._first_advance_done = True - max_dt = self._interface.advance(dt) + max_dt = self._participant.advance(dt) return max_dt def finalize(self): @@ -374,7 +383,7 @@ def finalize(self): ----- Refer finalize() in https://github.com/precice/python-bindings/blob/develop/precice.pyx """ - self._interface.finalize() + self._participant.finalize() def get_participant_name(self): """ @@ -398,7 +407,7 @@ def is_coupling_ongoing(self): tag : bool True if coupling is still going on and False if coupling has finished. """ - return self._interface.is_coupling_ongoing() + return self._participant.is_coupling_ongoing() def is_time_window_complete(self): """ @@ -413,55 +422,13 @@ def is_time_window_complete(self): tag : bool True if implicit coupling in the time window has converged and False if not converged yet. """ - return self._interface.is_time_window_complete() - - def is_action_required(self, action): - """ - Tag to check if a particular preCICE action is required. + return self._participant.is_time_window_complete() - Parameters - ---------- - action : string - Name of the preCICE action. - - Notes - ----- - Refer is_action_required(action) in https://github.com/precice/python-bindings/blob/develop/precice.pyx - - Returns - ------- - tag : bool - True if action is required and False if action is not required. - """ - return self._interface.is_action_required(action) - - def action_write_iteration_checkpoint(self): - """ - Get name of action to convey to preCICE that a checkpoint has been written. - - Notes - ----- - Refer action_write_iteration_checkpoint() in https://github.com/precice/python-bindings/blob/develop/precice.pyx - - Returns - ------- - action : string - Name of action related to writing a checkpoint. - """ - return precice.action_write_iteration_checkpoint() - - def action_read_iteration_checkpoint(self): - """ - Get name of action to convey to preCICE that a checkpoint has been read and the state of the system has been - restored to that checkpoint. - - Notes - ----- - Refer action_read_iteration_checkpoint() in https://github.com/precice/python-bindings/blob/develop/precice.pyx - - Returns - ------- - action : string - Name of action related to reading a checkpoint. - """ - return precice.action_read_iteration_checkpoint() + def get_max_time_step_size(self): + return self._participant.get_max_time_step_size() + + def requires_writing_checkpoint(self): + return self._participant.requires_writing_checkpoint() + + def requires_reading_checkpoint(self): + return self._participant.requires_reading_checkpoint() \ No newline at end of file diff --git a/tutorials/partitioned-heat-conduction/fenicsx/errorcomputation.py b/tutorials/partitioned-heat-conduction/fenicsx/errorcomputation.py index 069dd87..6e079d6 100644 --- a/tutorials/partitioned-heat-conduction/fenicsx/errorcomputation.py +++ b/tutorials/partitioned-heat-conduction/fenicsx/errorcomputation.py @@ -1,16 +1,24 @@ -from ufl import dx -from dolfinx.fem import assemble_scalar, form +from dolfinx import fem import numpy as np from mpi4py import MPI +import ufl def compute_errors(u_approx, u_ref, total_error_tol=10 ** -4): mesh = u_ref.function_space.mesh + # Compute L2 error and error at nodes + #V_ex = fem.functionspace(mesh, ("Lagrange", 2)) + #u_ex = fem.Function(V_ex) + #u_ex.interpolate(u_ref) + error_L2 = np.sqrt(mesh.comm.allreduce(fem.assemble_scalar(fem.form((u_approx - u_ref)**2 * ufl.dx)), op=MPI.SUM)) + if mesh.comm.rank == 0: + print(f"L2-error: {error_L2:.2e}") - # compute total L2 error between reference and calculated solution - error_pointwise = form(((u_approx - u_ref) / u_ref) ** 2 * dx) - error_total = np.sqrt(mesh.comm.allreduce(assemble_scalar(error_pointwise), MPI.SUM)) + # Compute values at mesh vertices + error_max = mesh.comm.allreduce(np.max(np.abs(u_approx.x.array - u_ref.x.array)), op=MPI.MAX) + if mesh.comm.rank == 0: + print(f"Error_max: {error_max:.2e}") + + assert (error_L2 < total_error_tol) - assert (error_total < total_error_tol) - - return error_total + return (error_L2, error_max) diff --git a/tutorials/partitioned-heat-conduction/fenicsx/heat.py b/tutorials/partitioned-heat-conduction/fenicsx/heat.py index 7ecef9f..2a0640f 100644 --- a/tutorials/partitioned-heat-conduction/fenicsx/heat.py +++ b/tutorials/partitioned-heat-conduction/fenicsx/heat.py @@ -1,42 +1,27 @@ """ -The basic example is taken from "Langtangen, Hans Petter, and Anders Logg. Solving PDEs in Python: The FEniCS -Tutorial I. Springer International Publishing, 2016." - -The example code has been extended with preCICE API calls and mixed boundary conditions to allow for a Dirichlet-Neumann -coupling of two separate heat equations. It also has been adapted to be compatible with FEniCSx. - -The original source code can be found on https://jsdokken.com/dolfinx-tutorial/chapter2/heat_code.html. - -Heat equation with Dirichlet conditions. (Dirichlet problem) - u'= Laplace(u) + f in the unit square [0,1] x [0,1] - u = u_C on the coupling boundary at x = 1 - u = u_D on the remaining boundary - u = u_0 at t = 0 - u = 1 + x^2 + alpha*y^2 + \beta*t - f = beta - 2 - 2*alpha - -Heat equation with mixed boundary conditions. (Neumann problem) - u'= Laplace(u) + f in the shifted unit square [1,2] x [0,1] - du/dn = f_N on the coupling boundary at x = 1 - u = u_D on the remaining boundary - u = u_0 at t = 0 - u = 1 + x^2 + alpha*y^2 + \beta*t - f = beta - 2 - 2*alpha +This code is mostly taken from: https://jsdokken.com/dolfinx-tutorial/chapter2/heat_equation.html """ -from __future__ import print_function, division +import basix.ufl +from petsc4py import PETSc +import ufl +from dolfinx import mesh, fem +from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc, LinearProblem +from ufl import TrialFunction, TestFunction, inner, dx, grad, ds +import basix + +import argparse +import numpy as np from mpi4py import MPI -from dolfinx.fem import Function, FunctionSpace, Expression, Constant, dirichletbc, locate_dofs_geometrical -from dolfinx.fem.petsc import LinearProblem -from dolfinx.io import XDMFFile -from ufl import TrialFunction, TestFunction, dx, ds, dot, grad, inner, lhs, rhs, FiniteElement, VectorElement +import sympy as sp + from fenicsxprecice import Adapter from errorcomputation import compute_errors + from my_enums import ProblemType, DomainPart -import argparse -import numpy as np from problem_setup import get_geometry +from dolfinx.io import VTXWriter def determine_gradient(V_g, u): """ @@ -48,206 +33,220 @@ def determine_gradient(V_g, u): w = TrialFunction(V_g) v = TestFunction(V_g) - a = inner(w, v) * dx - L = inner(grad(u), v) * dx + a = inner(w, v) * ufl.dx + L = inner(grad(u), v) * ufl.dx problem = LinearProblem(a, L) return problem.solve() +# Parse arguments parser = argparse.ArgumentParser(description="Solving heat equation for simple or complex interface case") -command_group = parser.add_mutually_exclusive_group(required=True) -command_group.add_argument("-d", "--dirichlet", help="create a dirichlet problem", dest="dirichlet", - action="store_true") -command_group.add_argument("-n", "--neumann", help="create a neumann problem", dest="neumann", action="store_true") -parser.add_argument("-e", "--error-tol", help="set error tolerance", type=float, default=10**-6,) - +parser.add_argument("participantName", help="Name of the solver.", type=str, choices=[p.value for p in ProblemType]) +parser.add_argument("-e", "--error-tol", help="set error tolerance", type=float, default=10**-8,) args = parser.parse_args() - -fenics_dt = .1 # time step size -# Error is bounded by coupling accuracy. In theory we would obtain the analytical solution. +# Init variables with arguments +participant_name = args.participantName error_tol = args.error_tol -alpha = 3 # parameter alpha -beta = 1.2 # parameter beta +t = 0 +fenics_dt = 0.1 +alpha = 3 +beta = 1.2 + -if args.dirichlet and not args.neumann: +# define the domain + +if participant_name == ProblemType.DIRICHLET.value: problem = ProblemType.DIRICHLET domain_part = DomainPart.LEFT -elif args.neumann and not args.dirichlet: +elif participant_name == ProblemType.NEUMANN.value: problem = ProblemType.NEUMANN domain_part = DomainPart.RIGHT -mesh, coupling_boundary, remaining_boundary = get_geometry(MPI.COMM_WORLD, domain_part) -# Define function space using mesh -scalar_element = FiniteElement("P", mesh.ufl_cell(), 2) -vector_element = VectorElement("P", mesh.ufl_cell(), 1) -V = FunctionSpace(mesh, scalar_element) -V_g = FunctionSpace(mesh, vector_element) -W = V_g.sub(0).collapse()[0] +# create domain and function space +domain, coupling_boundary, remaining_boundary = get_geometry(domain_part) +V = fem.functionspace(domain, ("Lagrange", 2)) +element = basix.ufl.element("Lagrange", domain.topology.cell_name(), 1, shape=(domain.geometry.dim,)) +V_g = fem.functionspace(domain, element) +W, map_to_W = V_g.sub(0).collapse() + +# Define the exact solution + -# Define boundary conditions +class exact_solution(): + def __init__(self, alpha, beta, t): + self.alpha = alpha + self.beta = beta + self.t = t + def __call__(self, x): + return 1 + x[0]**2 + self.alpha * x[1]**2 + self.beta * self.t -class Expression_u_D: - def __init__(self): - self.t = 0.0 - def eval(self, x): - return np.full(x.shape[1], 1 + x[0] * x[0] + alpha * x[1] * x[1] + beta * self.t) +u_exact = exact_solution(alpha, beta, t) +# Define the boundary condition +bcs = [] +u_D = fem.Function(V) +u_D.interpolate(u_exact) +tdim = domain.topology.dim +fdim = tdim - 1 +domain.topology.create_connectivity(fdim, tdim) +# dofs for the coupling boundary +dofs_coupling = fem.locate_dofs_geometrical(V, coupling_boundary) +# dofs for the remaining boundary. Can be directly set to u_D +dofs_remaining = fem.locate_dofs_geometrical(V, remaining_boundary) +bc_D = fem.dirichletbc(u_D, dofs_remaining) +bcs.append(bc_D) -u_D = Expression_u_D() -u_D_function = Function(V) -u_D_function.interpolate(u_D.eval) if problem is ProblemType.DIRICHLET: # Define flux in x direction - class Expression_f_N: - def __init__(self): - pass + f_N = fem.Function(W) + f_N.interpolate(lambda x: 2 * x[0]) - def eval(self, x): - return np.full(x.shape[1], 2 * x[0]) - f_N = Expression_f_N() - f_N_function = Function(V) - f_N_function.interpolate(f_N.eval) - -# Define initial value -u_n = Function(V, name="Temperature") -u_n.interpolate(u_D.eval) +u_n = fem.Function(V) # IV and solution u for the n-th time step +u_n.interpolate(u_exact) +# initialise precice precice, precice_dt, initial_data = None, 0.0, None - -# Initialize the adapter according to the specific participant if problem is ProblemType.DIRICHLET: - precice = Adapter(MPI.COMM_WORLD, adapter_config_filename="precice-adapter-config-D.json") - precice_dt = precice.initialize(coupling_boundary, read_function_space=V, write_object=f_N_function) -elif problem is ProblemType.NEUMANN: - precice = Adapter(MPI.COMM_WORLD, adapter_config_filename="precice-adapter-config-N.json") - precice_dt = precice.initialize(coupling_boundary, read_function_space=W, write_object=u_D_function) + precice = Adapter(adapter_config_filename="precice-adapter-config-D.json", mpi_comm=MPI.COMM_WORLD) +else: + precice = Adapter(adapter_config_filename="precice-adapter-config-N.json", mpi_comm=MPI.COMM_WORLD) -dt = Constant(mesh, 0.0) -dt.value = np.min([fenics_dt, precice_dt]) -# Define variational problem -u = TrialFunction(V) -v = TestFunction(V) +if problem is ProblemType.DIRICHLET: + precice.initialize(coupling_boundary, read_function_space=V, write_object=f_N) +elif problem is ProblemType.NEUMANN: + precice.initialize(coupling_boundary, read_function_space=W, write_object=u_D) +# get precice's dt +precice_dt = precice.get_max_time_step_size() +dt = np.min([fenics_dt, precice_dt]) -class Expression_f: - def __init__(self): - self.t = 0.0 - def eval(self, x): - return np.full(x.shape[1], beta - 2 - 2 * alpha) +# Define the variational formualation +# As $f$ is a constant independent of $t$, we can define it as a constant. +f = fem.Constant(domain, beta - 2 - 2 * alpha) -f = Expression_f() -f_function = Function(V) -F = u * v / dt * dx + dot(grad(u), grad(v)) * dx - (u_n / dt + f_function) * v * dx +# We can now create our variational formulation, with the bilinear form `a` and linear form `L`. -dofs_remaining = locate_dofs_geometrical(V, remaining_boundary) -bcs = [dirichletbc(u_D_function, dofs_remaining)] +u, v = ufl.TrialFunction(V), ufl.TestFunction(V) +F = u * v * ufl.dx + dt * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx - (u_n + dt * f) * v * ufl.dx +# create a coupling expression for the coupling_boundary and modify variational problem accordingly coupling_expression = precice.create_coupling_expression() -read_data = precice.read_data() -precice.update_coupling_expression(coupling_expression, read_data) if problem is ProblemType.DIRICHLET: # modify Dirichlet boundary condition on coupling interface - dofs_coupling = locate_dofs_geometrical(V, coupling_boundary) - bcs.append(dirichletbc(coupling_expression, dofs_coupling)) + bc_coup = fem.dirichletbc(coupling_expression, dofs_coupling) + bcs.append(bc_coup) if problem is ProblemType.NEUMANN: # modify Neumann boundary condition on coupling interface, modify weak # form correspondingly - F += v * coupling_expression * ds - -a, L = lhs(F), rhs(F) - -# Time-stepping -u_np1 = Function(V, name="Temperature") -t = 0 + F += dt*coupling_expression * v * ufl.ds -# reference solution at t=0 -u_ref = Function(V, name="reference") -u_ref.interpolate(u_D_function) +a = fem.form(ufl.lhs(F)) +L = fem.form(ufl.rhs(F)) -# Generating output files -temperature_out = XDMFFile(MPI.COMM_WORLD, f"./output/{precice.get_participant_name()}.xdmf", "w") -temperature_out.write_mesh(mesh) -ref_out = XDMFFile(MPI.COMM_WORLD, f"./output/ref{precice.get_participant_name()}.xdmf", "w") -ref_out.write_mesh(mesh) +# ## Create the matrix and vector for the linear problem +# To ensure that we are solving the variational problem efficiently, we +# will create several structures which can reuse data, such as matrix +# sparisty patterns. Especially note as the bilinear form `a` is +# independent of time, we only need to assemble the matrix once. -# output solution at t=0, n=0 -n = 0 +A = assemble_matrix(a, bcs=bcs) +A.assemble() +b = create_vector(L) +uh = fem.Function(V) -temperature_out.write_function(u_n, t) -ref_out.write_function(u_ref, t) +# ## Define a linear variational solver +# We will use [PETSc](https://www.mcs.anl.gov/petsc/) to solve the +# resulting linear algebra problem. We use the Python-API `petsc4py` to +# define the solver. We will use a linear solver. -u_D.t = t + dt.value -u_D_function.interpolate(u_D.eval) -f.t = t + dt.value -f_function.interpolate(f.eval) +solver = PETSc.KSP().create(domain.comm) +solver.setOperators(A) +solver.setType(PETSc.KSP.Type.PREONLY) +solver.getPC().setType(PETSc.PC.Type.LU) if problem is ProblemType.DIRICHLET: - flux = Function(V_g, name="Heat-Flux") + flux = fem.Function(V_g) + +# boundaries point as always to the end of the timestep +u_exact.t += dt +u_D.interpolate(u_exact) +# create writer for output files +vtxwriter = VTXWriter(MPI.COMM_WORLD, f"output_{problem.name}.bp", [u_n]) +vtxwriter.write(t) + while precice.is_coupling_ongoing(): - # write checkpoint - if precice.is_action_required(precice.action_write_iteration_checkpoint()): - precice.store_checkpoint(u_n, t, n) + if precice.requires_writing_checkpoint(): + precice.store_checkpoint(u_n, t, 0) - read_data = precice.read_data() + precice_dt = precice.get_max_time_step_size() + dt = np.min([fenics_dt, precice_dt]) - # Update the coupling expression with the new read data + read_data = precice.read_data(dt) + + # Update the right hand side reusing the initial vector + with b.localForm() as loc_b: + loc_b.set(0) + assemble_vector(b, L) + precice.update_coupling_expression(coupling_expression, read_data) - dt.value = np.min([fenics_dt, precice_dt]) + # Apply Dirichlet boundary condition to the vector (according to the tutorial, the lifting operation is used to preserve the symmetry of the matrix) + # As far as I understood, the boundary condition bc is updated by + # u_D.interpolate above, since this function is wrapped into the bc object + apply_lifting(b, [a], [bcs]) + set_bc(b, bcs) - linear_problem = LinearProblem(a, L, bcs=bcs) - u_np1 = linear_problem.solve() + # Solve linear problem + solver.solve(b, uh.x.petsc_vec) # Write data to preCICE according to which problem is being solved if problem is ProblemType.DIRICHLET: # Dirichlet problem reads temperature and writes flux on boundary to Neumann problem - flux = determine_gradient(V_g, u_np1) - flux_x = Function(W) + flux = determine_gradient(V_g, uh) + flux_x = fem.Function(W) flux_x.interpolate(flux.sub(0)) precice.write_data(flux_x) + #precice.write_data(f_N) elif problem is ProblemType.NEUMANN: # Neumann problem reads flux and writes temperature on boundary to Dirichlet problem - precice.write_data(u_np1) + precice.write_data(uh) - precice_dt = precice.advance(dt.value) + precice.advance(dt) + precice_dt = precice.get_max_time_step_size() # roll back to checkpoint - if precice.is_action_required(precice.action_read_iteration_checkpoint()): - u_cp, t_cp, n_cp = precice.retrieve_checkpoint() - u_n.interpolate(u_cp) + if precice.requires_reading_checkpoint(): + u_cp, t_cp, _ = precice.retrieve_checkpoint() + u_n.x.array[:] = u_cp.x.array t = t_cp - n = n_cp else: # update solution - u_n.interpolate(u_np1) - t += dt.value - n += 1 + # Update solution at previous time step (u_n) + u_n.x.array[:] = uh.x.array + t += float(dt) + vtxwriter.write(t) if precice.is_time_window_complete(): - u_ref.interpolate(u_D_function) - error = compute_errors(u_n, u_ref, total_error_tol=error_tol) - print('n = %d, t = %.2f: L2 error on domain = %.3g' % (n, t, error)) - print('output u^%d and u_ref^%d' % (n, n)) - - temperature_out.write_function(u_n, t) - ref_out.write_function(u_ref, t) + u_ref = fem.Function(V) + u_ref.interpolate(u_D) + error, error_pointwise = compute_errors(u_n, u_ref, total_error_tol=1) + print("t = %.2f: L2 error on domain = %.3g" % (t, error)) + + # Update Dirichlet BC + u_exact.t += dt + u_D.interpolate(u_exact) + # TODO: update time dependent f (as soon as it is time dependent)! - # Update Dirichlet BC - u_D.t = t + dt.value - u_D_function.interpolate(u_D.eval) - f.t = t + dt.value - f_function.interpolate(f.eval) - - -# Hold plot precice.finalize() + +vtxwriter.close() diff --git a/tutorials/partitioned-heat-conduction/fenicsx/my_enums.py b/tutorials/partitioned-heat-conduction/fenicsx/my_enums.py index 5708bad..d5e74db 100644 --- a/tutorials/partitioned-heat-conduction/fenicsx/my_enums.py +++ b/tutorials/partitioned-heat-conduction/fenicsx/my_enums.py @@ -5,8 +5,8 @@ class ProblemType(Enum): """ Enum defines problem type. Details see above. """ - DIRICHLET = 1 # Dirichlet problem - NEUMANN = 2 # Neumann problem + DIRICHLET = "Dirichlet" # Dirichlet problem + NEUMANN = "Neumann" # Neumann problem class DomainPart(Enum): diff --git a/tutorials/partitioned-heat-conduction/fenicsx/problem_setup.py b/tutorials/partitioned-heat-conduction/fenicsx/problem_setup.py index 9dbbf74..37127a2 100644 --- a/tutorials/partitioned-heat-conduction/fenicsx/problem_setup.py +++ b/tutorials/partitioned-heat-conduction/fenicsx/problem_setup.py @@ -2,8 +2,10 @@ Problem setup for partitioned-heat-conduction/fenicsx tutorial """ from dolfinx.mesh import DiagonalType, create_rectangle +import dolfinx.mesh from my_enums import DomainPart import numpy as np +from mpi4py import MPI y_bottom, y_top = 0, 1 @@ -24,11 +26,8 @@ def straight_boundary(x): return np.isclose(x[0], x_coupling, tol) -def get_geometry(mpi_comm, domain_part): +def get_geometry(domain_part): nx = ny = 9 - low_resolution = 5 - high_resolution = 5 - n_vertices = 20 if domain_part is DomainPart.LEFT: p0 = (x_left, y_bottom) @@ -38,8 +37,7 @@ def get_geometry(mpi_comm, domain_part): p1 = (x_right, y_top) else: raise Exception("invalid domain_part: {}".format(domain_part)) - - mesh = create_rectangle(mpi_comm, [p0, p1], [nx, ny], diagonal=DiagonalType.left) + mesh = create_rectangle(MPI.COMM_WORLD, [np.asarray(p0), np.asarray(p1)], [nx, ny], dolfinx.mesh.CellType.triangle) coupling_boundary = straight_boundary remaining_boundary = exclude_straight_boundary diff --git a/tutorials/partitioned-heat-conduction/fenicsx/run.sh b/tutorials/partitioned-heat-conduction/fenicsx/run.sh index e31f07a..59efc2c 100755 --- a/tutorials/partitioned-heat-conduction/fenicsx/run.sh +++ b/tutorials/partitioned-heat-conduction/fenicsx/run.sh @@ -4,10 +4,10 @@ set -e -u while getopts ":dn" opt; do case ${opt} in d) - python3 heat.py -d --error-tol 10e-3 + python3 heat.py Dirichlet --error-tol 10e-3 ;; n) - python3 heat.py -n --error-tol 10e-3 + python3 heat.py Neumann --error-tol 10e-3 ;; \?) echo "Usage: cmd [-d] [-n]"