diff --git a/README.md b/README.md index 8a1f5b9..f9bbd5e 100644 --- a/README.md +++ b/README.md @@ -20,7 +20,8 @@ Notes: * This adapter is a fork from the original [FEniCS-preCICE adapter](https://github.com/precice/fenics-adapter). Based on [v1.2.0](https://github.com/precice/fenics-adapter/releases/tag/v1.2.0). * This adapter is currently under development and experimental. -* Target version: dolfinx v0.5.2. Other version that have been tested and might also work: v0.4.1 +* Target version: dolfinx v0.8.0. +* Target version: preCICE v3. ## Installing the package @@ -44,7 +45,7 @@ Make sure to install the following dependencies: * python3 (this adapter **only supports python3**) * [the python language bindings for preCICE](https://github.com/precice/python-bindings) * :construction: [FEniCSx](https://fenicsproject.org/) (with python interface, installed by default) (under construction refer to notes on FEniCSx below) :construction: -* and scipy (`pip3 install scipy`) + #### Build and install the adapter @@ -90,4 +91,6 @@ If you are using FEniCSx, please also consider the information on [the official 2021: For development of FEniCSx support, `precice/fenics-adapter@v1.2.0` was forked as `precice/fenicsx-adapter`. The required modifications were carried out by [Benjamin Rodenberg](https://www.cs.cit.tum.de/sccs/personen/benjamin-rodenberg/) and [Ishaan Desai](https://www.ipvs.uni-stuttgart.de/institute/team/Desai/). -2023: [Philip Hildebrand](https://github.com/PhilipHildebrand) updated the adapter to a [first minimal working version](https://github.com/precice/fenicsx-adapter/pull/15) and contributed a [first tutorial](https://github.com/precice/tutorials/pull/317) in the scope of his Bachelor's thesis ["Extending the FEniCSx Adapter for the Coupling Library preCICE"](https://mediatum.ub.tum.de/node?id=1706280) under supervision of [Benjamin Rodenberg](https://www.cs.cit.tum.de/sccs/personen/benjamin-rodenberg/) and [Ishaan Desai](https://www.ipvs.uni-stuttgart.de/institute/team/Desai/). +2023: [Philip Hildebrand](https://github.com/PhilipHildebrand) updated the adapter to a first minimal working version (https://github.com/precice/fenicsx-adapter/pull/15) and contributed a first tutorial (https://github.com/precice/tutorials/pull/317) in the scope of his Bachelor's thesis ["Extending the FEniCSx Adapter for the Coupling Library preCICE"](https://mediatum.ub.tum.de/node?id=1706280) under supervision of [Benjamin Rodenberg](https://www.cs.cit.tum.de/sccs/personen/benjamin-rodenberg/) and [Ishaan Desai](https://www.ipvs.uni-stuttgart.de/institute/team/Desai/). + +2024: [Eduardo Firvida](https://github.com/efirvida) updated the adapter to a minimal working version (https://github.com/precice/fenicsx-adapter/pull/26) with FSI and 3D support for preCICEv3. In the scope of his PhD Thesis. diff --git a/fenicsxprecice/adapter_core.py b/fenicsxprecice/adapter_core.py deleted file mode 100644 index 17db94f..0000000 --- a/fenicsxprecice/adapter_core.py +++ /dev/null @@ -1,180 +0,0 @@ -""" -This module consists of helper functions used in the Adapter class. Names of the functions are self explanatory -""" - -from dolfinx.fem import FunctionSpace, Function -import numpy as np -from enum import Enum -import logging -import copy - -logger = logging.getLogger(__name__) -logger.setLevel(level=logging.INFO) - - -class Vertices: - """ - Vertices class provides a generic skeleton for vertices. A set of vertices has a set of IDs and - coordinates as defined in FEniCSx. - """ - - def __init__(self): - self._ids = None - self._coordinates = None - - def set_ids(self, ids): - self._ids = ids - - def set_coordinates(self, coords): - self._coordinates = coords - - def get_ids(self): - return copy.deepcopy(self._ids) - - def get_coordinates(self): - return copy.deepcopy(self._coordinates) - - -class FunctionType(Enum): - """ - Defines scalar- and vector-valued function. - Used in assertions to check if a FEniCSx function is scalar or vector. - """ - SCALAR = 0 # scalar valued function - VECTOR = 1 # vector valued function - - -class CouplingMode(Enum): - """ - Defines the type of coupling being used. - Options are: Bi-directional coupling, Uni-directional Write Coupling, Uni-directional Read Coupling - Used in assertions to check which type of coupling is done - """ - BI_DIRECTIONAL_COUPLING = 4 - UNI_DIRECTIONAL_WRITE_COUPLING = 5 - UNI_DIRECTIONAL_READ_COUPLING = 6 - - -def determine_function_type(input_obj): - """ - Determines if the function is scalar- or vector-valued based on rank evaluation. - - Parameters - ---------- - input_obj : - A FEniCSx function. - - Returns - ------- - 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 input_obj.num_sub_spaces == 0: - return FunctionType.SCALAR - elif input_obj.num_sub_spaces == 2: - return FunctionType.VECTOR - elif isinstance(input_obj, Function): - if len(input_obj.x.array.shape) == 1: - return FunctionType.SCALAR - elif input_obj.x.array.shape[1] > 1: - return FunctionType.VECTOR - else: - raise Exception("Error determining type of given dolfin Function") - else: - raise Exception("Error determining type of given dolfin FunctionSpace") - - -def convert_fenicsx_to_precice(fenicsx_function, local_ids): - """ - Converts data of type dolfin.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. - - Returns - ------- - precice_data : array_like - Array of FEniCSx function values at each point on the boundary. - """ - - if not isinstance(fenicsx_function, 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([]) - - return np.array(precice_data) - - -def get_fenicsx_vertices(function_space, coupling_subdomain, dims): - """ - Extracts vertices which FEniCSx accesses on this rank and which lie on the given coupling domain, from a given - function space. - - Parameters - ---------- - function_space : FEniCSx function space - Function space on which the finite element problem definition lives. - coupling_subdomain : FEniCSx Domain - Subdomain consists of only the coupling interface region. - dims : int - Dimension of problem. - - Returns - ------- - ids : numpy array - Array of ids of fenicsx vertices. - coords : numpy array - The coordinates of fenicsx vertices in a numpy array [N x D] where - N = number of vertices and D = dimensions of geometry. - """ - - # 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) - if dims == 2: - coords = mesh.geometry.x[ids][:, :2] - else: - coords = np.array([]) - except Exception as e: # fall back to old method # TODO is that too general? Better use, e.g., IndexError here? - print("Caught the following exception in the detection of the coupling subdomain:\n{e}") - print("Falling back to old, point-wise method.") - ids, coords = [], [] - for idx in range(mesh.geometry.x.shape[0]): - v = mesh.geometry.x[idx] - if coupling_subdomain(v): - ids.append(idx) - if dims == 2: - coords.append([v[0], v[1]]) - ids = np.array(ids) - coords = np.array(coords) - return ids, coords diff --git a/fenicsxprecice/config.py b/fenicsxprecice/config.py index ba7c938..894f78f 100644 --- a/fenicsxprecice/config.py +++ b/fenicsxprecice/config.py @@ -35,7 +35,7 @@ def read_json(self, adapter_config_filename): """ folder = os.path.dirname(os.path.join(os.getcwd(), os.path.dirname(sys.argv[0]), adapter_config_filename)) path = os.path.join(folder, os.path.basename(adapter_config_filename)) - read_file = open(path, "r") + read_file = open(path) data = json.load(read_file) self._config_file_name = os.path.join(folder, data["config_file_name"]) self._participant_name = data["participant_name"] diff --git a/fenicsxprecice/expression_core.py b/fenicsxprecice/expression_core.py deleted file mode 100644 index 35e6ea4..0000000 --- a/fenicsxprecice/expression_core.py +++ /dev/null @@ -1,147 +0,0 @@ -""" -This module provides a mechanism to imterpolate point data acquired from preCICE into FEniCSx Expressions. -""" - -from .adapter_core import FunctionType -from scipy.interpolate import Rbf -from scipy.linalg import lstsq -from dolfinx.fem import Function -import numpy as np -from mpi4py import MPI - -import logging - -logger = logging.getLogger(__name__) -logger.setLevel(level=logging.INFO) - - -class CouplingExpression(Function): - """ - Creates functional representation (for FEniCSx) of nodal data provided by preCICE. - """ - - def __init__(self, function_space, function_type): - self._dimension = function_space.mesh.geometry.dim - self._function_type = function_type - super().__init__(function_space) - - def update_boundary_data(self, vals, coords_x, coords_y=None, coords_z=None): - """ - Update object of this class of type FEniCSx UserExpression with given point data. - - Parameters - ---------- - vals : double - Point data to be used to update the Expression. - coords_x : double - X coordinate of points of which point data is provided. - coords_y : double - Y coordinate of points of which point data is provided. - coords_z : double - Z coordinate of points of which point data is provided. - """ - self._coords_x = coords_x - if coords_y is None: - coords_y = np.zeros(self._coords_x.shape) - self._coords_y = coords_y - if coords_z is None: - coords_z = np.zeros(self._coords_x.shape) - - self._coords_y = coords_y - self._coords_z = coords_z - self._vals = vals - - interpolant = self._create_interpolant() - - if self._is_scalar_valued(): - assert (self._vals.shape[0] == self._coords_x.shape[0]) - elif self._is_vector_valued(): - assert (self._vals.shape[0] == self._coords_x.shape[0]) - - self.interpolate(interpolant) - - def _create_interpolant(self, x): - # TODO: the correct way to deal with this would be using an abstract class. Since this is technically more - # complex and the current implementation is a workaround anyway, we do not - # use the proper solution, but this hack. - """ - Creates interpolant from boundary data that has been provided before. - - Parameters - ---------- - x : double - Point. - - Returns - ------- - list : python list - Interpolant as a list. If scalar function is interpolated this list has a single - element. If a vector function is interpolated the list has self._dimensions elements. - """ - raise Exception("Please use one of the classes derived from this class, that implements an actual strategy for" - "interpolation.") - - def _is_scalar_valued(self): - """ - Determines if function being interpolated is scalar-valued based on dimension of provided vector self._vals. - - Returns - ------- - tag : bool - True if function being interpolated is scalar-valued, False otherwise. - """ - return self._function_type is FunctionType.SCALAR - - def _is_vector_valued(self): - """ - Determines if function being interpolated is vector-valued based on dimension of provided vector self._vals. - - Returns - ------- - tag : bool - True if function being interpolated is vector-valued, False otherwise. - """ - return self._function_type is FunctionType.VECTOR - - -class SegregatedRBFInterpolationExpression(CouplingExpression): - """ - Uses polynomial quadratic fit + RBF interpolation for implementation of CustomExpression.interpolate. Allows for - arbitrary coupling interfaces. - - See Lindner, F., Mehl, M., & Uekermann, B. (2017). Radial basis function interpolation for black-box multi-physics - simulations. - """ - - def _segregated_interpolant_2d(self, coords_x, coords_y, data): - assert (coords_x.shape == coords_y.shape) - # create least squares system to approximate a * x ** 2 + b * x + c ~= y - - def lstsq_interp(x, y, w): return w[0] * x ** 2 + w[1] * y ** 2 + w[2] * x * y + w[3] * x + w[4] * y + w[5] - - A = np.empty((coords_x.shape[0], 0)) - n_unknowns = 6 - for i in range(n_unknowns): - w = np.zeros([n_unknowns]) - w[i] = 1 - column = lstsq_interp(coords_x, coords_y, w).reshape((coords_x.shape[0], 1)) - A = np.hstack([A, column]) - - # solve system - w, _, _, _ = lstsq(A, data) - # create fit - - # compute remaining error - res = data - lstsq_interp(coords_x, coords_y, w) - # add RBF for error - rbf_interp = Rbf(coords_x, coords_y, res) - - return lambda x, y: rbf_interp(x, y) + lstsq_interp(x, y, w) - - def _create_interpolant(self): - """ - See base class description. - """ - assert (self._is_scalar_valued()) # this implementation only supports scalar valued functions - assert (self._dimension == 2) # this implementation only supports two dimensions - return lambda x: self._segregated_interpolant_2d(self._coords_x, self._coords_y, self._vals)(x[0], x[1]) diff --git a/fenicsxprecice/fenicsxprecice.py b/fenicsxprecice/fenicsxprecice.py index cf05b38..d84b5ba 100644 --- a/fenicsxprecice/fenicsxprecice.py +++ b/fenicsxprecice/fenicsxprecice.py @@ -1,44 +1,36 @@ -""" -FEniCSx - preCICE Adapter. API to help users couple FEniCS with other solvers using the preCICE library. -:raise ImportError: if PRECICE_ROOT is not defined -""" -from os import write -import numpy as np -from .config import Config +import copy import logging +from typing import List, NamedTuple + +import dolfinx as dfx +import numpy as np import precice -from .adapter_core import FunctionType, determine_function_type, convert_fenicsx_to_precice, get_fenicsx_vertices, CouplingMode, Vertices -from .expression_core import SegregatedRBFInterpolationExpression + +from .config import Config from .solverstate import SolverState -from dolfinx.fem import Function, FunctionSpace -import copy logger = logging.getLogger(__name__) logger.setLevel(level=logging.INFO) -class Adapter: - """ - This adapter class provides an interface to the preCICE coupling library for setting up a coupling case which has - FEniCSx as a participant for 2D problems. - The user can create and manage a dolfinx.UserExpression and/or dolfinx.PointSource at the coupling boundary. - Reading data from preCICE and writing data to preCICE is also managed via functions of this class. - If the user wants to perform implicit coupling then a steering mechanism for checkpointing is also provided. +class Dimensions(NamedTuple): + interface: int + mesh: int - For more information on setting up a coupling case using dolfinx.UserExpression at the coupling boundary please have - a look at this tutorial: - TODO - For more information on setting up a coupling case using dolfinx.PointSource at the coupling boundary please have a - look at this tutorial: - TODO +class Adapter: + """This adapter class provides an interface to the preCICE v3 coupling library for setting up a + coupling case with FEniCSx as a participant in 2D and 3D problems. + + The coupling is achieved using the FunctionSpace degrees of freedom (DOFs) on the interface. + Data interchange between participants (read/write) is performed by accessing the values on + the interface DOFs. This approach allows us to leverage the full range of FEniCSx functionalities + while ensuring seamless communication with other participants in the coupling process. - NOTE: dolfinx.PointSource use only works in serial """ - def __init__(self, mpi_comm, adapter_config_filename='precice-adapter-config.json'): - """ - Constructor of Adapter class. + def __init__(self, mpi_comm, adapter_config_filename="precice-adapter-config.json"): + """Constructor of Adapter class. Parameters ---------- @@ -47,30 +39,22 @@ def __init__(self, mpi_comm, adapter_config_filename='precice-adapter-config.jso adapter_config_filename : string Name of the JSON adapter configuration file (to be provided by the user) """ - - self._config = Config(adapter_config_filename) - - # Setup up MPI communicator self._comm = mpi_comm + self._config = Config(adapter_config_filename) - self._interface = precice.Interface(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 - self._dofmap = None # initialized later using function space provided by user + self._interface = precice.Participant( + self._config.get_participant_name(), + self._config.get_config_file_name(), + self._comm.Get_rank(), + self._comm.Get_size(), + ) # coupling mesh related quantities - self._fenicsx_vertices = Vertices() - self._precice_vertex_ids = None # initialized later + self._fenicsx_vertices = None + self._precice_vertex_ids = None - # read data related quantities (read data is read from preCICE and applied in FEniCSx) - self._read_function_type = None # stores whether read function is scalar or vector valued - self._write_function_type = None # stores whether write function is scalar or vector valued - - # Interpolation strategy - self._my_expression = SegregatedRBFInterpolationExpression + # problem Function Space + self._function_space = None # Solver state used by the Adapter internally to handle checkpointing self._checkpoint = None @@ -84,127 +68,78 @@ def __init__(self, mpi_comm, adapter_config_filename='precice-adapter-config.jso # Problem dimension in FEniCSx self._fenicsx_dims = None - def create_coupling_expression(self): - """ - Creates a FEniCSx Expression in the form of an object of class GeneralInterpolationExpression or - ExactInterpolationExpression. The adapter will hold this object till the coupling is on going. + def interpolation_points_in_vector_space(self): + """Determine interpolation points in the vector space. Returns ------- - coupling_expression : Object of class dolfinx.functions.expression.Expression - Reference to object of class GeneralInterpolationExpression or ExactInterpolationExpression. + tuple + Unrolled degrees of freedom and their coordinates. """ + V = self._function_space + bs = V.dofmap.bs - if not (self._read_function_type is FunctionType.SCALAR or self._read_function_type is FunctionType.VECTOR): - raise Exception("No valid read_function is provided in initialization. Cannot create coupling expression") + fs_coords = V.tabulate_dof_coordinates() + fs_coords = fs_coords[:, : self.dimensions.mesh] + boundary_dofs = dfx.fem.locate_dofs_topological(V, self.dimensions.mesh - 1, self._tags) + unrolled_dofs = np.empty(len(boundary_dofs) * bs, dtype=np.int32) - coupling_expression = self._my_expression(self._read_function_space, self._read_function_type) + for i, dof in enumerate(boundary_dofs): + for b in range(bs): + unrolled_dofs[i * bs + b] = dof * bs + b - return coupling_expression - - def update_coupling_expression(self, coupling_expression, data): - """ - Updates the given FEniCSx Expression using provided data. The boundary data is updated. - User needs to explicitly call this function in each time step. - - Parameters - ---------- - 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 - value. - """ - vertices = np.array(list(data.keys())) - nodal_data = np.array(list(data.values())) - coupling_expression.update_boundary_data(nodal_data, vertices[:, 0], vertices[:, 1]) - - def get_point_sources(self, data): - raise Exception("PointSources are not implemented for the FEniCSx adapter.") + return unrolled_dofs.reshape(-1, bs), fs_coords[boundary_dofs] def read_data(self): - """ - 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 - For a vector read function the data is a numpy array with shape (N, D) where - N = number of coupling vertices and D = dimensions of FEniCSx setup + """Read data from preCICE. - Note: For quasi 2D-3D coupled simulation (FEniCSx participant is 2D) the Z-component of the data and vertices - is deleted. + Incoming data is a ndarray where the shape of the array depends on the dimensions of the problem. + For scalar problems, this will be a 1D array (vector), while for vector problems, + it will be an Mx2 array (in 2D) or an Mx3 array (in 3D), where M is the number of interface nodes. Returns ------- - data : dict_like - The coupling data. A dictionary containing nodal data with vertex coordinates as key and associated data as - value. - """ - 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)} - + np.ndarray + The incoming data containing nodal data ordered according to _fenicsx_vertices + """ + mesh_name = self._config.get_coupling_mesh_name() + data_name = self._config.get_read_data_name() + # For more information about readDara see and and time start here: + # https://precice.org/couple-your-code-porting-v2-3.html#add-relativereadtime-for-all-read-data-calls + read_data = self._interface.read_data( + mesh_name, data_name, self._precice_vertex_ids, self.dt + ) return copy.deepcopy(read_data) def write_data(self, write_function): - """ - Writes data to preCICE. Depending on the dimensions of the simulation (2D-3D Coupling, 2D-2D coupling or - Scalar/Vector write function) write_data is first converted into a format needed for preCICE. + """Writes data to preCICE. Depending on the dimensions of the simulation. + For scalar problems, this will be a 1D array (vector), while for vector problems, + it will be an Mx2 array (in 2D) or an Mx3 array (in 3D), where M is the number of interface nodes. + Parameters ---------- - write_function : Object of class dolfinx.functions.function.Function - A FEniCSx function consisting of the data which this participant will write to preCICE in every time step. + write_function : dolfinx.fem.Function + A FEniCSx function consisting of the data which this participant will write to preCICE + in every time step. """ - - assert (self._coupling_type is CouplingMode.UNI_DIRECTIONAL_WRITE_COUPLING or - CouplingMode.BI_DIRECTIONAL_COUPLING) - - w_func = write_function - - # 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) - - 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_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") + mesh_name = self._config.get_coupling_mesh_name() + write_data_name = self._config.get_write_data_name() + write_data = write_function.x.array[self.interface_dof] + self._interface.write_data(mesh_name, write_data_name, self._precice_vertex_ids, write_data) def initialize(self, coupling_subdomain, read_function_space=None, write_object=None): - """ - Initializes the coupling and sets up the mesh where coupling happens in preCICE. + """Initializes the coupling and sets up the mesh where coupling happens in preCICE. Parameters ---------- - coupling_subdomain : Object of class dolfinx.cpp.mesh.SubDomain - SubDomain of mesh which is the physical coupling boundary. - read_function_space : Object of class dolfinx.functions.functionspace.FunctionSpace + coupling_subdomain : List + Indices of entities representing the coupling interface normally face sets tags. + read_function_space : dolfinx.fem.FunctionSpace Function space on which the read function lives. If not provided then the adapter assumes that this participant is a write-only participant. - write_object : Object of class dolfinx.functions.functionspace.FunctionSpace / dolfinx.functions.function.Function - Function space on which the write function lives or FEniCSx function related to the quantity to be written + write_object : dolfinx.fem.Function + FEniCSx function related to the quantity to be written by FEniCSx during each coupling iteration. If not provided then the adapter assumes that this participant is a read-only participant. @@ -213,118 +148,35 @@ def initialize(self, coupling_subdomain, read_function_space=None, write_object= dt : double Recommended time step value from preCICE. """ - - write_function_space, write_function = None, None - if isinstance(write_object, Function): # precice.initialize_data() will be called using this Function - write_function_space = write_object.function_space - write_function = write_object - elif isinstance(write_object, FunctionSpace): # preCICE will use default zero values for initialization. - write_function_space = write_object - write_function = None - elif write_object is None: - pass - else: - raise Exception("Given write object is neither of type dolfinx.functions.function.Function or " - "dolfinx.functions.functionspace.FunctionSpace") - - if isinstance(read_function_space, FunctionSpace): - pass - elif read_function_space is None: - pass - else: - raise Exception("Given read_function_space is not of type dolfinx.functions.functionspace.FunctionSpace") - - if read_function_space is None and write_function_space: - self._coupling_type = CouplingMode.UNI_DIRECTIONAL_WRITE_COUPLING - assert (self._config.get_write_data_name()) - print("Participant {} is write-only participant".format(self._config.get_participant_name())) - function_space = write_function_space - elif read_function_space and write_function_space is None: - self._coupling_type = CouplingMode.UNI_DIRECTIONAL_READ_COUPLING - assert (self._config.get_read_data_name()) - print("Participant {} is read-only participant".format(self._config.get_participant_name())) - function_space = read_function_space - elif read_function_space and write_function_space: - self._coupling_type = CouplingMode.BI_DIRECTIONAL_COUPLING - assert (self._config.get_read_data_name() and self._config.get_write_data_name()) - function_space = read_function_space - elif read_function_space is None and write_function_space is None: - raise Exception("Neither read_function_space nor write_function_space is provided. Please provide a " - "write_object if this participant is used in one-way coupling and only writes data. " - "Please provide a read_function_space if this participant is used in one-way coupling and " - "only reads data. If two-way coupling is implemented then both read_function_space" - " and write_object need to be provided.") - else: - raise Exception("Incorrect read and write function space combination provided. Please check input " - "parameters in initialization") - - if self._coupling_type is CouplingMode.UNI_DIRECTIONAL_READ_COUPLING or \ - self._coupling_type is CouplingMode.BI_DIRECTIONAL_COUPLING: - self._read_function_type = determine_function_type(read_function_space) - self._read_function_space = read_function_space - - if self._coupling_type is CouplingMode.UNI_DIRECTIONAL_WRITE_COUPLING or \ - self._coupling_type is CouplingMode.BI_DIRECTIONAL_COUPLING: - # 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 - - self._fenicsx_dims = function_space.mesh.geometry.dim - - # Ensure that function spaces of read and write functions use the same mesh - if self._coupling_type is CouplingMode.BI_DIRECTIONAL_COUPLING: - assert (self._read_function_space.mesh is write_function_space.mesh - ), "read_function_space and write_object need to be defined using the same mesh" - - if self._fenicsx_dims != 2: - raise Exception("Currently the fenicsx-adapter only supports 2D cases") - - if self._fenicsx_dims != self._interface.get_dimensions(): - 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 not write_function: - raise Exception("Non-standard initialization requires a write_function") - self.write_data(write_function) - self._interface.mark_action_fulfilled(precice.action_write_initial_data()) - - self._interface.initialize_data() - - return precice_dt - - def store_checkpoint(self, user_u, 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 - Current state of the physical variable of interest for this participant. - t : double - Current simulation time. - n : int - Current time window (iteration) number. - """ + self._domain = read_function_space.mesh + self._function_space = read_function_space + self._tags = coupling_subdomain + + self._interface_dof, self._interface_dof_coords = ( + self.interpolation_points_in_vector_space() + ) + self._fenicsx_vertices = list(zip(self._interface_dof, self._interface_dof_coords)) + + self._precice_vertex_ids = self._interface.set_mesh_vertices( + self._config.get_coupling_mesh_name(), self._interface_dof_coords[:, : self.dimensions.interface] + ) + + if self._interface.requires_initial_data(): + self._interface.write_data( + self._config.get_coupling_mesh_name(), + self._config.get_write_data_name(), + self._precice_vertex_ids, + np.zeros(self._interface_dof_coords.shape), + ) + self._interface.initialize() + return self.dt + + def store_checkpoint(self, states: List): + """Defines an object of class SolverState which stores the current states of the variable and the time stamp.""" if self._first_advance_done: - assert (self.is_time_window_complete()) - + assert self.is_time_window_complete() logger.debug("Store checkpoint") - my_u = user_u.copy() - # making sure that the FEniCSx function provided by user is not directly accessed by the Adapter - assert (my_u != user_u) - self._checkpoint = SolverState(my_u, t, n) - self._interface.mark_action_fulfilled(self.action_write_iteration_checkpoint()) + self._checkpoint = SolverState(states) def retrieve_checkpoint(self): """ @@ -332,21 +184,15 @@ def retrieve_checkpoint(self): Returns ------- - u : FEniCSx Function - Current state of the physical variable of interest for this participant. - t : double - Current simulation time. - n : int - Current time window (iteration) number. + tuple + The stored checkpoint state (u, v, a, t). """ - assert (not self.is_time_window_complete()) + 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): - """ - Advances coupling in preCICE. + def advance(self, dt: float): + """Advances coupling in preCICE. Parameters ---------- @@ -401,67 +247,62 @@ def is_coupling_ongoing(self): return self._interface.is_coupling_ongoing() def is_time_window_complete(self): - """ - Tag to check if implicit iteration has converged. + """Tag to check if implicit iteration has converged. - Notes + Notes: ----- Refer is_time_window_complete() in https://github.com/precice/python-bindings/blob/develop/precice.pyx - Returns + Returns: ------- 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. - - 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 + def requires_reading_checkpoint(self): + """Checks if reading a checkpoint is required. - Returns + 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): + bool + True if reading a checkpoint is required, False otherwise. """ - Get name of action to convey to preCICE that a checkpoint has been written. + return self._interface.requires_reading_checkpoint() - Notes - ----- - Refer action_write_iteration_checkpoint() in https://github.com/precice/python-bindings/blob/develop/precice.pyx + def requires_writing_checkpoint(self): + """Checks if writing a checkpoint is required. - Returns + 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() + bool + True if writing a checkpoint is required, False otherwise. + """ + return self._interface.requires_writing_checkpoint() + + @property + def interface_dof(self): + """Returns the interface degrees of freedom.""" + return self._interface_dof + + @property + def interface_coordinates(self): + """Returns the interface coordinates.""" + return self._interface_dof_coords + + @property + def precice(self): + """Returns the preCICE interface object.""" + return self._interface + + @property + def dimensions(self) -> Dimensions: + """Returns the dimensions of the mesh and the interface.""" + return Dimensions( + self._interface.get_mesh_dimensions(self._config.get_coupling_mesh_name()), self._domain.geometry.dim + ) + + @property + def dt(self): + """Returns the maximum time step size allowed by preCICE.""" + return self._interface.get_max_time_step_size() diff --git a/fenicsxprecice/solverstate.py b/fenicsxprecice/solverstate.py index 26fe44e..1cf6082 100644 --- a/fenicsxprecice/solverstate.py +++ b/fenicsxprecice/solverstate.py @@ -1,36 +1,22 @@ -class SolverState: - def __init__(self, u, t, n): - """ - Solver state consists of a value u, associated time t and the timestep n +import copy +from collections.abc import Sequence - Parameters - ---------- - u : Object of class dolfinx.functions.function.Function - FEniCSx function related to the field during each coupling iteration. - t : double - Time stamp. - n : int - Iteration number. - """ - self.u = u - self.t = t - self.n = n +import dolfinx as dfx - def get_state(self): - """ - Returns the state variables value u, associated time t and timestep n - Returns - ------- - u : Object of class dolfinx.functions.function.Function - A copy of FEniCSx function related to the field during each coupling iteration. - t : double - Time stamp. - n : int - Iteration number. +class SolverState: + def __init__(self, states: Sequence): + """Store a list of function as states for those + iteration that not converge """ - return self.u.copy(), self.t, self.n + states_cp = [] + for state in states: + if isinstance(state, dfx.fem.Function): + states_cp.append(state.copy()) + else: + states_cp.append(copy.deepcopy(state)) + self.__state = states_cp - def print_state(self): - u, t, n = self.get_state() - return print("u={u}, t={t}, n={n}".format(u=u, t=t, n=n)) + def get_state(self): + """Returns the state of the solver.""" + return self.__state diff --git a/setup.py b/setup.py index 0a24b23..73f710b 100644 --- a/setup.py +++ b/setup.py @@ -39,7 +39,7 @@ author_email='info@precice.org', license='LGPL-3.0', packages=['fenicsxprecice'], - install_requires=['pyprecice~=2.0.0', 'scipy', 'numpy>=1.13.3', 'mpi4py'], + install_requires=['pyprecice~=3.0.0', 'numpy>=1.13.3'], tests_require=['sympy'], test_suite='tests', zip_safe=False) diff --git a/tutorials/perpendicular-flap/solid-fenicsx/solid.py b/tutorials/perpendicular-flap/solid-fenicsx/solid.py new file mode 100644 index 0000000..fc791d7 --- /dev/null +++ b/tutorials/perpendicular-flap/solid-fenicsx/solid.py @@ -0,0 +1,260 @@ +import os + +import dolfinx as dfx +import numpy as np +import ufl +from dolfinx.fem.petsc import ( + LinearProblem, + apply_lifting, + assemble_matrix_mat, + assemble_vector, + set_bc, +) +from dolfinx.mesh import CellType, create_rectangle +from fenicsxprecice import Adapter +from mpi4py import MPI +from petsc4py import PETSc + +# ------ # +# SOLVER # +# ------ # + +# EXTEND the linear petsc.LinearProblem to add values into the bilinear form matrix on +# desired DOFs. Wee need this due to the incoming Force are discrete vector, so wee need +# to avoid domain integration + + +class DiscreteLinearProblem(LinearProblem): + def __init__( + self, + a, + L, + bcs, + u, + point_dofs, + petsc_options=None, + form_compiler_options=None, + jit_options=None, + ): + if point_dofs is not None: + if not isinstance(type(point_dofs), np.ndarray): + point_dofs = np.array([point_dofs], dtype="int32") + self._point_dofs = point_dofs + + super().__init__(a, L, bcs, u, petsc_options, form_compiler_options, jit_options) + + def solve(self, values): + """Solve the problem.""" + # Assemble lhs + self._A.zeroEntries() + assemble_matrix_mat(self._A, self._a, bcs=self.bcs) + self._A.assemble() + + # Assemble rhs + with self._b.localForm() as b_loc: + b_loc.set(0) + + if self._point_dofs is not None and values is not None: + # apply load in dofs + self._b[self._point_dofs] = values + self._b.assemble() + + assemble_vector(self._b, self._L) + + # Apply boundary conditions to the rhs + apply_lifting(self._b, [self._a], bcs=[self.bcs]) + self._b.ghostUpdate(addv=PETSc.InsertMode.ADD, mode=PETSc.ScatterMode.REVERSE) + set_bc(self._b, self.bcs) + + # Solve linear system and update ghost values in the solution + self._solver.solve(self._b, self._x) + self.u.x.scatter_forward() + + return self.u + + +# --------- # +# CONSTANTS # +# --------- # + +MPI_COMM = MPI.COMM_WORLD +CURRENT_FOLDER = os.path.dirname(__file__) +PARTICIPANT_CONFIG = os.path.join(CURRENT_FOLDER, "precice-config.json") +RESULTS_DIR = os.path.join(CURRENT_FOLDER, "results") +os.makedirs(RESULTS_DIR, exist_ok=True) +os.chdir(CURRENT_FOLDER) + +WRITER = dfx.io.VTKFile(MPI_COMM, f"{RESULTS_DIR}/result.pvd", "w") + +WIDTH, HEIGHT = 0.1, 1 +NX, NY = 2, 15 * 8 + +E = 4000000.0 +NU = 0.3 +RHO = 3000.0 + +BETA_ = 0.25 +GAMMA_ = 0.5 + +# ------- # +# MESHING # +# ------- # + +domain = create_rectangle( + MPI_COMM, + [np.array([-WIDTH / 2, 0]), np.array([WIDTH / 2, HEIGHT])], + [NX, NY], + cell_type=CellType.quadrilateral, +) +dim = domain.topology.dim + +# -------------- # +# Function Space # +# -------------- # +degree = 2 +shape = (dim,) +V = dfx.fem.functionspace(domain, ("P", degree, shape)) +u = dfx.fem.Function(V, name="Displacement") +f = dfx.fem.Function(V, name="Force") + +# ------------------- # +# Boundary conditions # +# ------------------- # +tol = 1e-14 + + +def clamped_boundary(x): + return abs(x[1]) < tol + + +def neumann_boundary(x): + """Determines whether a node is on the coupling boundary.""" + return np.logical_or((np.abs(x[1] - HEIGHT) < tol), np.abs(np.abs(x[0]) - WIDTH / 2) < tol) + + +fixed_boundary = dfx.fem.locate_dofs_geometrical(V, clamped_boundary) +coupling_boundary = dfx.mesh.locate_entities_boundary(domain, dim - 1, neumann_boundary) + +bcs = [dfx.fem.dirichletbc(np.zeros((dim,)), fixed_boundary, V)] + +# ------------ # +# PRECICE INIT # +# ------------ # +participant = Adapter(MPI_COMM, PARTICIPANT_CONFIG) +participant.initialize(coupling_boundary, V, u) +dt = participant.dt + +# ------------------------ # +# linear elastic equations # +# ------------------------ # +E = dfx.fem.Constant(domain, E) +nu = dfx.fem.Constant(domain, NU) +rho = dfx.fem.Constant(domain, RHO) + +lmbda = E * nu / (1 + nu) / (1 - 2 * nu) +mu = E / 2 / (1 + nu) + + +def epsilon(v): + return ufl.sym(ufl.grad(v)) + + +def sigma(v): + return lmbda * ufl.tr(epsilon(v)) * ufl.Identity(dim) + 2 * mu * epsilon(v) + + +# ------------------- # +# Time discretization # +# ------------------- # +# prev time step +u_old = dfx.fem.Function(V) +v_old = dfx.fem.Function(V) +a_old = dfx.fem.Function(V) + +# current time step +a_new = dfx.fem.Function(V) +v_new = dfx.fem.Function(V) + +beta = dfx.fem.Constant(domain, BETA_) +gamma = dfx.fem.Constant(domain, GAMMA_) + +dx = ufl.Measure("dx", domain=domain) + +a = (1 / (beta * dt**2)) * (u - u_old - dt * v_old) - ((1 - 2 * beta) / (2 * beta)) * a_old +a_expr = dfx.fem.Expression(a, V.element.interpolation_points()) + +v = v_old + dt * ((1 - gamma) * a_old + gamma * a) +v_expr = dfx.fem.Expression(v, V.element.interpolation_points()) + +# ------------------ # +# mass, a stiffness # +# ------------------ # +u_ = ufl.TestFunction(V) +du = ufl.TrialFunction(V) + + +def mass(u, u_): + return rho * ufl.dot(u, u_) * dx + + +def stiffness(u, u_): + return ufl.inner(sigma(u), epsilon(u_)) * dx + + +Residual = mass(a, u_) + stiffness(u, u_) +Residual_du = ufl.replace(Residual, {u: du}) +a_form = ufl.lhs(Residual_du) +L_form = ufl.rhs(Residual_du) + + +problem = DiscreteLinearProblem( + a=a_form, L=L_form, u=u, bcs=bcs, point_dofs=participant.interface_dof +) + + +# parameters for Time-Stepping +t = 0.0 + +while participant.is_coupling_ongoing(): + if participant.requires_writing_checkpoint(): # write checkpoint + participant.store_checkpoint((u_old, v_old, a_old, t)) + + read_data = participant.read_data() + + u = problem.solve(read_data) + + # Write new displacements to preCICE + participant.write_data(u) + + # Call to advance coupling, also returns the optimum time step value + participant.advance(dt) + + # Either revert to old step if timestep has not converged or move to next timestep + if participant.requires_reading_checkpoint(): + ( + u_cp, + v_cp, + a_cp, + t_cp, + ) = participant.retrieve_checkpoint() + u_old.vector.copy(u_cp.vector) + v_old.vector.copy(v_cp.vector) + a_old.vector.copy(a_cp.vector) + t = t_cp + else: + v_new.interpolate(v_expr) + a_new.interpolate(a_expr) + u.vector.copy(u_old.vector) + v_new.vector.copy(v_old.vector) + a_new.vector.copy(a_old.vector) + t += dt + + if participant.is_time_window_complete(): + f.vector[participant.interface_dof] = read_data + f.vector.assemble() + WRITER.write_function(u, t) + WRITER.write_function(f, t) + WRITER.close() + +WRITER.close() +participant.finalize()