FSI with FenicsX and preCICE 3

Dear all, I dont know if here is the rigth place to post this, but as I need quick help I think that here is a good better place than in github isses.

I need to implement FSi coupling preCICE+OpenFOAM+FenicsX, but as I understand the fenicsx adapter have all the limitantions that I need :grin:

1- don’t support FSI at the moment
2- is only for 2D
3- is not ported yet for preCICEv3

So I decideto doit by myself based on the current fenicsx adapter but with a lot of simplifications. Here is my code at the moment:

import precice
import os
import json
import sys
import numpy as np
import dolfinx as dfx
from petsc4py import PETSc
import logging

logger = logging.getLogger("precice")


class Config:
    """
    Handles reading of config. parameters of the fenicsxadapter based on JSON
    configuration file. Initializer calls read_json() method. Instance attributes
    can be accessed by provided getter functions.
    """

    def __init__(self, config_path: str):
        self._config_file = config_path
        self._precice_config = None
        self._participant = None
        self._mesh = None
        self._read_data = None
        self._write_data = None
        self._patch_tags = []

        self.read_json(config_path)

    def read_json(self, config_path: str):
        """
        Reads JSON adapter configuration file and saves the data to the respective instance attributes.

        Parameters
        ----------
        config_path : string
            Name of the JSON configuration file
        """
        read_file = open(config_path, "r")
        data = json.load(read_file)
        folder = os.path.dirname(
            os.path.join(os.getcwd(), os.path.dirname(sys.argv[0]), config_path)
        )
        self._precice_config = os.path.join(folder, data["precice_config"])
        self._participant = data["participant"]
        self._mesh = data["interface"]["mesh"]

        try:
            self._write_data = data["interface"]["write_data"]
        except KeyError:
            # not required for one-way coupling, if this participant reads data
            self._write_data = None

        try:
            self._read_data = data["interface"]["read_data"]
        except KeyError:
            # not required for one-way coupling, if this participant writes data
            self._read_data = None

        read_file.close()

    @property
    def precice_config(self):
        return self._precice_config

    @property
    def participant(self):
        return self._participant

    @property
    def mesh(self):
        return self._mesh

    @property
    def patch_tags(self):
        return self._patch_tags

    @property
    def read_data(self):
        return self._read_data

    @property
    def write_data(self):
        return self._write_data


class FunctionData:
    def __init__(self, vector_space, interface_dof):
        self.V = vector_space
        self._dim = self.V.mesh.geometry.dim
        self._v_space_coords = self.V.tabulate_dof_coordinates()[:, : self._dim]
        self._interface_dof = interface_dof
        self._vector_size = len(self._v_space_coords.flatten())

        self._vector = PETSc.Vec().createMPI(self._vector_size, comm=PETSc.COMM_WORLD)
        self._vector.setFromOptions()
        self._vector.setUp()
        self.reset()

    def reset(self):
        self._vector.set(0.0)
        self._vector.assemble()

    def setValues(self, values):
        self.reset()
        arr = self.vector.getArray()
        arr = arr.reshape(-1, self._dim)
        arr[self._interface_dof] = values

        start, end = self._vector.getOwnershipRange()
        self._vector[start:end] = arr[start:end]
        self._vector.assemble()

    def __str__(self):
        if self._dim == 3:
            Fx, Fy, Fz = self.array.sum(axis=0)
            return f"Sum Fx: {Fx:.2f}, Sum Fy: {Fy:.2f}, Sum Fz:{Fz:.2f}"
        if self._dim == 2:
            Fx, Fy = self.array.sum(axis=0)
            return f"Sum Fx: {Fx:.2f}, Sum Fy: {Fy:.2f}"

    @property
    def array(self):
        return self._vector.getArray().reshape(-1, self._dim)

    @property
    def coords(self):
        fs_coords = self.V.tabulate_dof_coordinates()
        return fs_coords[self._interface_dof]

    @property
    def vector(self):
        return self._vector


class SolverState:
    def __init__(self, u, v, a, t, n):
        self.u = u
        self.v = v
        self.a = a
        self.t = t
        self.n = n

    def get_state(self):
        return self.u, self.v, self.a, self.t, self.n

    def __str__(self):
        u, v, a, t, n = self.get_state()
        return f"u={u}, v={u}, a={u}, t={t}, n={n}"


class Adapter:
    def __init__(self, mpi_comm, config_path: str, mesh) -> None:
        self._config = Config(config_path)

        self._comm = mpi_comm
        self._domain = mesh
        self._topology = self._domain.topology
        self._checkpoint = None

        self._interface_nodes_coords = None
        self._function_space = None

        self._interface = precice.Participant(
            self._config.participant,
            self._config.precice_config,
            self._comm.Get_rank(),
            self._comm.Get_size(),
        )
        self._dim = self._interface.get_mesh_dimensions(self._config.mesh)
        self._fdim = self._domain.geometry.dim - 1
        assert (
            self._dim == self._domain.geometry.dim
        ), "Mesh in the config file and the mesh provided have different dimensions"

    def interpolation_points_in_vector_space(self):
        V = self._function_space
        fs_coords = V.tabulate_dof_coordinates()
        fs_coords = fs_coords[:, : self._dim]
        boundary_dofs = dfx.fem.locate_dofs_topological(V, self._fdim, self._facesets)
        boundary_coords = fs_coords[boundary_dofs]
        return boundary_dofs, boundary_coords

    def initialize(self, function_space, facesets):
        self._function_space = function_space
        self._facesets = facesets
        self._facesets_tags = dfx.mesh.meshtags(self._domain, self._fdim, np.sort(facesets), 1)

        self._interface_dof, self._interface_dof_coords = (
            self.interpolation_points_in_vector_space()
        )

        self._precice_vertex_ids = self._interface.set_mesh_vertices(
            self._config.mesh, self._interface_dof_coords
        )

        self._function_data = FunctionData(self._function_space, self._interface_dof)

        if self._interface.requires_initial_data():
            self._interface.write_data(
                self._config.mesh,
                self._config.write_data,
                self._precice_vertex_ids,
                np.zeros(self._interface_dof_coords.shape),
            )
        self._interface.initialize()
        return self._function_data

    @property
    def interface_ids(self):
        return self._interface_dof

    @property
    def interface_coordinates(self):
        return self._interface_dof_coords

    @property
    def precice(self):
        return self._interface

    @property
    def dt(self):
        return self._interface.get_max_time_step_size()

     def is_coupling_ongoing(self):
        return self._interface.is_coupling_ongoing()

    def is_time_window_complete(self):
        return self._interface.is_time_window_complete()

    def requires_reading_checkpoint(self):
        return self._interface.requires_reading_checkpoint()

    def requires_writing_checkpoint(self):
        return self._interface.requires_writing_checkpoint()

    def advance(self, dt: float):
        return self._interface.advance(dt)

    def finalize(self):
        return self._interface.finalize()

    def read_data(self, dt):
        mesh_name = self._config.mesh
        data_name = self._config.read_data
        read_data = self._interface.read_data(mesh_name, data_name, self._precice_vertex_ids, dt)
        return read_data

    def write_data(self, write_function):
        mesh_name = self._config.mesh
        write_data_name = self._config.write_data
        write_data = write_function.x.array.reshape(-1, self._dim)[self.interface_ids]
        self._interface.write_data(mesh_name, write_data_name, self._precice_vertex_ids, write_data)

    def store_checkpoint(self, u, v, a, t, n):
        self._checkpoint = SolverState(u.copy(), v.copy(), a.copy(), t, n)

    def retrieve_checkpoint(self):
        assert not self.is_time_window_complete()
        return self._checkpoint.get_state()

And flap tutorial code will looks like this:

import numpy as np
import ufl
import os
from mpi4py import MPI

import dolfinx as dfx
import dolfinx.fem.petsc
from dolfinx.mesh import create_rectangle, CellType
from petsc4py import PETSc
from dolfinx.fem.petsc import LinearProblem

# --------- #
# 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 = 4, 26


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 
WRITER.write_mesh(domain)

# -------------- #
# Function Space #
# -------------- #
degree = 2
shape = (dim,)  # this means we want a vector field of size `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)
coupling_boundary_tags = dfx.mesh.meshtags(domain, dim-1, np.sort(coupling_boundary), 1)

bcs = [dfx.fem.dirichletbc(np.zeros((dim,)), fixed_boundary, V)]

# ------------ #
# PRECICE INIT #
# ------------ #
participant = Adapter(MPI_COMM, PARTICIPANT_CONFIG, domain)
_f = participant.initialize(V, coupling_boundary)
with open(f"{RESULTS_DIR}/interpolation_points.csv", "w") as p_file:
    p_file.write("X,Y,Z\n")
    np.savetxt(p_file, _f.coords, delimiter=",")
dt = dfx.fem.Constant(domain, PETSc.ScalarType(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)
ds = ufl.Measure("ds", subdomain_data=coupling_boundary_tags)

a = 1 / beta / dt**2 * (u - u_old - dt * v_old) + a_old * (1 - 1 / 2 / beta)
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

L = ufl.dot(f, u_) * ds(1)  # include wind loads over surface

Residual = mass(a, u_) + stiffness(u, u_) - L

Residual_du = ufl.replace(Residual, {u: du})
a_form = ufl.lhs(Residual_du)
L_form = ufl.rhs(Residual_du)

opts = PETSc.Options()  # type: ignore
opts["ksp_type"] = "cg"
opts["ksp_rtol"] = 1.0e-8
opts["pc_type"] = "gamg"

# Use Chebyshev smoothing for multigrid
opts["mg_levels_ksp_type"] = "chebyshev"
opts["mg_levels_pc_type"] = "jacobi"

# Improve estimate of eigenvalues for Chebyshev smoothing
opts["mg_levels_ksp_chebyshev_esteig_steps"] = 10

problem = LinearProblem(a_form, L_form, u=u, bcs=bcs, petsc_options=opts.getAll())


# parameters for Time-Stepping
t = 0.0
n = 0

while participant.is_coupling_ongoing():

    if participant.requires_writing_checkpoint():  # write checkpoint
        participant.store_checkpoint(u_old, v_old, a_old, t, n)

    dt.value = participant.dt

    read_data = participant.read_data(dt)

    start1, end1 = _f.vector.getOwnershipRange()
    start2, end2 = f.vector.getOwnershipRange()
    _f.setValues(read_data)
    f.vector[start2:end2] = _f.vector[start1: end1]
    f.vector.assemble()
    
    problem.solve()
    f.x.scatter_forward()
    u.x.scatter_forward()

    # Write new displacements to preCICE
    participant.write_data(u)

    # Call to advance coupling, also returns the optimum time step value
    participant.advance(float(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, n_cp = participant.retrieve_checkpoint()
        u_old.interpolate(u_cp)
        v_old.interpolate(v_cp)
        a_old.interpolate(a_cp)
        t = t_cp
        n = n_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
        n += 1
        print("===================================================================")
        u_values =u.x.array.reshape((-1, dim)).sum(axis=0)
        print(_f)
        if dim == 3:
            print(f"Sum Ux: {u_values[0]:.2e}, Sum Uy: {u_values[1]:.2e},Sum Uz: {u_values[2]:.2e}")
        if dim == 2:
            print(f"Sum Ux: {u_values[0]:.2e}, Sum Uy: {u_values[1]:.2e}")
        print("===================================================================")

    if participant.is_time_window_complete():
        if n % 10 == 0:
            WRITER.write_function(u, t)
            WRITER.write_function(f, t)
            WRITER.close()

WRITER.close()
participant.finalize()

And it runs until the end but the results are very inaccurate, and I can’t figure put where is the mistake.

can someone help me with this?

Update
seems that I have an scaling problem here, in this next chart I multiply the Displacement in fenix *100, and seems to be OK. I test with other points and seems to be reproducible.

I hav e a question about wathc point. the data of the wathc point are the received data frome the other participant or are they interpolated on the mesh receptor mesh point. Im, asking this, because If I plot the Forces for calculix and for FenicsX in te same point they are not to far. So my gess is that there is something wrong on my solver or in the way that I map the force to the boundary

1 Like

Hi @efirvida!

Nice to see that you are using FEniCSx with preCICE. Did you already have a look at the thesis of Philip Hildebrand? He developed the FEniCSx adapter. Unfortunately, the adapter is not in very good shape at the moment because I am mainly focusing on FEniCS until I finish my thesis.

Some ideas how we could help you better:

  • Feel free to open a pull request with your code from above. This helps us to compare it to the version of the FEniCSx adapter we have in the repository.
  • Similarly for the perpendicular flap that you used for testing: Feel free to open a pull request in the tutorials repository. After having a brief look at your code: I think this would be very useful and something we could integrate into the tutorials after some polishing and testing.

Some general comments:

  • Regarding your scaling problem: Did you consider using PointSource? I remember this existed in FEniCS and was not available in FEniCSx at the time when Philip wrote his thesis. The forces that the FEniCS participant reads from OpenFOAM are discrete forces on the mesh nodes and not a load distribution. Therefore, applying the forces as a Neumann boundary condtion (L = ufl.dot(f, u_) * ds(1)) is probably not the right approach here. We also wrote a paragraph about this topic in our FEniCS adapter paper.
  • You can also find the partitioned heat conduction case from Philip’s thesis on this branch.
  • If you want to come to the preCICE workshop: This would be a great opportunity to get things started and work together

Greetings!
Benjamin

Hi @BenjaminRodenberg

Thanks for your response, I will continue this threat on a pull request.

Thanks for the advice of applying the force as Newman you have a right, I found that by myself, but it’s good to have a second opinion.

Sadly the FenicsX doesn’t have a PointSource, so I use another approach of modifying the Function vector directly which works for a single point, but I still have problems mapping all loads over the interface. Further discussion in the pull request :).