diff --git a/partitioned-heat-conduction-3d/README.md b/partitioned-heat-conduction-3d/README.md new file mode 100644 index 000000000..bc193ea26 --- /dev/null +++ b/partitioned-heat-conduction-3d/README.md @@ -0,0 +1,57 @@ +--- +title: 3D partitioned heat conduction +permalink: tutorials-partitioned-heat-conduction-3d.html +keywords: FEniCSx, Heat conduction +summary: We solve a simple heat equation on a 3D domain. The domain is partitioned and the coupling is established in a Dirichlet-Neumann fashion. +--- + +{% note %} +Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/develop/partitioned-heat-conduction-3d), as continuously rendered here, or see the [latest released version](https://github.com/precice/tutorials/tree/master/partitioned-heat-conduction-3d) (if there is already one). Read how in the [tutorials introduction](https://precice.org/tutorials.html). +{% endnote %} + +## Setup + +We solve a partitioned heat equation. For information on the two dimensional non-partitioned case, please refer to [1, p.37ff]. In this tutorial the computational domain is partitioned and coupled via preCICE. The coupling roughly follows the approach described in [2]. + +![Case setup of partitioned-heat-conduction case](images/tutorials-partitioned-heat-conduction-setup.png) + +The computational domain can be seen in the picture above. To get a better overview of the domain, the domains of the participants are translated. +The domain of the Neumann participant is the small box, and the Dirichlet participant's domain is the rest. In the simulation, the small box is fully contained in the Dirichlet domain. This means, the coupling interface consists of five sides of the box. + +## Configuration + +preCICE configuration (image generated using the [precice-config-visualizer](https://precice.org/tooling-config-visualization.html)): + +![preCICE configuration visualization](images/tutorials-partitioned-heat-conduction-precice-config.png) + +## Available solvers and dependencies + +You can either couple a solver with itself or different solvers with each other. In any case you will need to have preCICE and the python bindings installed on your system. + +* FEniCSx. Install [FEniCS](https://fenicsproject.org/download/) and the [FEniCSx-adapter](https://github.com/precice/fenicsx-adapter). The code is adapted from the existing [fenics-tutorial](https://github.com/hplgit/fenics-tutorial/blob/master/pub/python/vol1/ft03_heat.py) from [1]. + +## Running the simulation + +You can find the corresponding `run.sh` script for running the case in the folders corresponding to the participant you want to use: + +```bash +cd dirichlet-fenicsx +./run.sh +``` + +and + +```bash +cd neumann-fenicsx +./run.sh +``` + +## Visualization + +Output is written into the folders of the fenicsx solvers (`neumann-fenicsx/output-neumann.bp` and `dirichlet-fenicsx/output-dirichlet.bp`). + +It is sufficient to import the folders to paraview to get the visualization of the simulation. + +## References + +[1] Hans Petter Langtangen and Anders Logg. "Solving PDEs in Minutes-The FEniCS Tutorial Volume I." (2016). [pdf](https://fenicsproject.org/pub/tutorial/pdf/fenics-tutorial-vol1.pdf) diff --git a/partitioned-heat-conduction-3d/clean-tutorial.sh b/partitioned-heat-conduction-3d/clean-tutorial.sh new file mode 100755 index 000000000..56e40d689 --- /dev/null +++ b/partitioned-heat-conduction-3d/clean-tutorial.sh @@ -0,0 +1,10 @@ +#!/usr/bin/env sh +set -e -u + +# shellcheck disable=SC1091 +. ../tools/cleaning-tools.sh + +clean_tutorial . +clean_precice_logs . +rm -fv ./*.log + diff --git a/partitioned-heat-conduction-3d/dirichlet-fenicsx/clean.sh b/partitioned-heat-conduction-3d/dirichlet-fenicsx/clean.sh new file mode 100755 index 000000000..9628d477a --- /dev/null +++ b/partitioned-heat-conduction-3d/dirichlet-fenicsx/clean.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env sh +set -e -u + +rm -rfv ./precice-profiling +rm -rfv ./output-dirichlet.bp +rm -fv ./*.log diff --git a/partitioned-heat-conduction-3d/dirichlet-fenicsx/precice-adapter-config-D.json b/partitioned-heat-conduction-3d/dirichlet-fenicsx/precice-adapter-config-D.json new file mode 100644 index 000000000..2dcb3a9c1 --- /dev/null +++ b/partitioned-heat-conduction-3d/dirichlet-fenicsx/precice-adapter-config-D.json @@ -0,0 +1,10 @@ +{ + "participant_name": "Dirichlet", + "config_file_name": "../precice-config.xml", + "interfaces": { + "Dirichlet-Mesh": { + "write_data_names": ["Heat-Flux"], + "read_data_names": ["Temperature"] + } + } +} diff --git a/partitioned-heat-conduction-3d/dirichlet-fenicsx/run.sh b/partitioned-heat-conduction-3d/dirichlet-fenicsx/run.sh new file mode 100755 index 000000000..74b46579e --- /dev/null +++ b/partitioned-heat-conduction-3d/dirichlet-fenicsx/run.sh @@ -0,0 +1,8 @@ +#!/bin/sh +set -e -u + +python3 -m venv --system-site-packages .venv +. .venv/bin/activate +pip install -r ../solver-fenicsx/requirements.txt + +python3 ../solver-fenicsx/heat.py Dirichlet --error-tol 10e-3 diff --git a/partitioned-heat-conduction-3d/images/tutorials-partitioned-heat-conduction-3d-setup.png b/partitioned-heat-conduction-3d/images/tutorials-partitioned-heat-conduction-3d-setup.png new file mode 100644 index 000000000..70a1d12df Binary files /dev/null and b/partitioned-heat-conduction-3d/images/tutorials-partitioned-heat-conduction-3d-setup.png differ diff --git a/partitioned-heat-conduction-3d/neumann-fenicsx/clean.sh b/partitioned-heat-conduction-3d/neumann-fenicsx/clean.sh new file mode 100755 index 000000000..1982adbce --- /dev/null +++ b/partitioned-heat-conduction-3d/neumann-fenicsx/clean.sh @@ -0,0 +1,6 @@ +#!/usr/bin/env sh +set -e -u + +rm -rfv ./precice-profiling +rm -rfv ./output-neumann.bp +rm -fv ./*.log diff --git a/partitioned-heat-conduction-3d/neumann-fenicsx/precice-adapter-config-N.json b/partitioned-heat-conduction-3d/neumann-fenicsx/precice-adapter-config-N.json new file mode 100644 index 000000000..f3987fe72 --- /dev/null +++ b/partitioned-heat-conduction-3d/neumann-fenicsx/precice-adapter-config-N.json @@ -0,0 +1,10 @@ +{ + "participant_name": "Neumann", + "config_file_name": "../precice-config.xml", + "interfaces": { + "Neumann-Mesh": { + "write_data_names": ["Temperature"], + "read_data_names": ["Heat-Flux"] + } + } +} diff --git a/partitioned-heat-conduction-3d/neumann-fenicsx/run.sh b/partitioned-heat-conduction-3d/neumann-fenicsx/run.sh new file mode 100755 index 000000000..008e81ee6 --- /dev/null +++ b/partitioned-heat-conduction-3d/neumann-fenicsx/run.sh @@ -0,0 +1,7 @@ +#!/bin/sh +set -e -u + +python3 -m venv --system-site-packages .venv +. .venv/bin/activate +pip install -r ../solver-fenicsx/requirements.txt +python3 ../solver-fenicsx/heat.py Neumann --error-tol 10e-3 diff --git a/partitioned-heat-conduction-3d/precice-config.xml b/partitioned-heat-conduction-3d/precice-config.xml new file mode 100644 index 000000000..00961dea6 --- /dev/null +++ b/partitioned-heat-conduction-3d/precice-config.xml @@ -0,0 +1,80 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/partitioned-heat-conduction-3d/solver-fenicsx/clean.sh b/partitioned-heat-conduction-3d/solver-fenicsx/clean.sh new file mode 100755 index 000000000..3a8b4619d --- /dev/null +++ b/partitioned-heat-conduction-3d/solver-fenicsx/clean.sh @@ -0,0 +1,6 @@ +#!/bin/sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_fenics . diff --git a/partitioned-heat-conduction-3d/solver-fenicsx/errorcomputation.py b/partitioned-heat-conduction-3d/solver-fenicsx/errorcomputation.py new file mode 100644 index 000000000..c3d239216 --- /dev/null +++ b/partitioned-heat-conduction-3d/solver-fenicsx/errorcomputation.py @@ -0,0 +1,21 @@ +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 ** -7): + mesh = u_ref.function_space.mesh + # Compute L2 error and error at nodes + 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 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) + + return (error_L2, error_max) diff --git a/partitioned-heat-conduction-3d/solver-fenicsx/heat.py b/partitioned-heat-conduction-3d/solver-fenicsx/heat.py new file mode 100644 index 000000000..310e164ae --- /dev/null +++ b/partitioned-heat-conduction-3d/solver-fenicsx/heat.py @@ -0,0 +1,311 @@ +""" +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, mixed boundary conditions to allow for a Dirichlet-Neumann +coupling of two separate heat equations and the domain has been changed to three dimensions. +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_equation.html. + +Heat equation with Dirichlet conditions. (Dirichlet problem) + u'= Laplace(u) + f in the outer domain + u = u_C on the coupling boundary + 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 inner domain + du/dn = f_N on the coupling boundary + 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 +""" +import argparse +import numpy as np +from mpi4py import MPI +import basix.ufl +from petsc4py import PETSc +import ufl +from dolfinx import fem, io +from dolfinx.fem.petsc import assemble_matrix, assemble_vector, apply_lifting, create_vector, set_bc +import basix +from fenicsxprecice import Adapter, CouplingMesh +from errorcomputation import compute_errors +from my_enums import ProblemType, DomainPart +from problem_setup import get_geometry + + +class GradientSolver: + """ + compute flux following http://hplgit.github.io/INF5620/doc/pub/fenics_tutorial1.1/tu2.html#tut-poisson-gradu + The solver has been changed since the original version from the link above introduces larger errors + + :param V_g: Vector function space + :param u: solution where gradient is to be determined + """ + + def __init__(self, domain, V_g): + self.domain = domain, + self.V_g = V_g + + w = ufl.TrialFunction(V_g) + self.v = ufl.TestFunction(V_g) + a = fem.form(ufl.inner(w, self.v) * ufl.dx) + self.A = assemble_matrix(a) + self.A.assemble() + + self.solver = PETSc.KSP().create(domain.comm) + self.solver.setOperators(self.A) + self.solver.setType(PETSc.KSP.Type.PREONLY) + self.solver.getPC().setType(PETSc.PC.Type.LU) + + self.returnValue = fem.Function(V_g) + + def compute(self, u): + L = fem.form(ufl.inner(ufl.grad(u), self.v) * ufl.dx) + b = create_vector(fem.extract_function_spaces(L)) + assemble_vector(b, L) + # In the assembly above, ranks did not communicate and therefore, the cells have incorrect values + # To resolve this, the values in the ghost region must be added to the values from the owner + b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + self.solver.solve(b, self.returnValue.x.petsc_vec) + return self.returnValue + + +# Parse arguments +parser = argparse.ArgumentParser(description="Solving heat equation for simple or complex interface case") +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() +# Init variables with arguments +participant_name = args.participantName +error_tol = args.error_tol + +comm = MPI.COMM_WORLD +rank = comm.Get_rank() + +t = 0 +fenics_dt = 0.1 +alpha = 3 +beta = .25 +gamma = 1.2 + +# define the domain +if participant_name == ProblemType.DIRICHLET.value: + problem = ProblemType.DIRICHLET + domain_part = DomainPart.OUTER +elif participant_name == ProblemType.NEUMANN.value: + problem = ProblemType.NEUMANN + domain_part = DomainPart.INNER + +# create domain and function space +domain, coupling_boundary, remaining_boundary = get_geometry(domain_part, comm) +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) +V_coup = None +if problem is ProblemType.DIRICHLET: + V_coup = V +else: + V_coup = V_g + + +class exact_solution(): + # Define the exact solution + def __init__(self, alpha, beta, gamma, t): + self.alpha = alpha + self.beta = beta + self.gamma = gamma + self.t = t + + def __call__(self, x): + return 1 + x[0]**2 + self.alpha * x[1]**2 + self.beta * x[2]**2 + self.gamma * self.t + + +u_exact = exact_solution(alpha, beta, gamma, t) + +gradient_solver = GradientSolver(domain, V_g) + +# 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 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) + +# dofs (and the corresponding coordinates) for the coupling boundary (coupling function space) +# this is later used to read data from preCICE and update boundary condition +dofs_coupling = fem.locate_dofs_geometrical(V_coup, coupling_boundary) +dofs_coupling_coordinates = V_coup.tabulate_dof_coordinates()[dofs_coupling] + +if problem is ProblemType.DIRICHLET: + f_N = fem.Function(V_g) + f_N.interpolate(gradient_solver.compute(u_D)) + +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 +if problem is ProblemType.DIRICHLET: + precice = Adapter(adapter_config_filename="precice-adapter-config-D.json", mpi_comm=comm) +else: + precice = Adapter(adapter_config_filename="precice-adapter-config-N.json", mpi_comm=comm, digit_cutoff=12) + +coupling_mesh = None +if problem is ProblemType.DIRICHLET: + coupling_mesh = CouplingMesh("Dirichlet-Mesh", coupling_boundary, {"Temperature": V_coup}, {"Heat-Flux": f_N}) + precice.initialize([coupling_mesh]) +elif problem is ProblemType.NEUMANN: + coupling_mesh = CouplingMesh("Neumann-Mesh", coupling_boundary, {"Heat-Flux": V_coup}, {"Temperature": u_D}) + precice.initialize([coupling_mesh]) + +# get precice's dt +precice_dt = precice.get_max_time_step_size() +dt = np.min([fenics_dt, precice_dt]) + + +# Define the variational formualation +# As $f$ is a constant independent of $t$, we can define it as a constant. +f = fem.Constant(domain, gamma - 2 - 2 * alpha - 2 * beta) +# We can now create our variational formulation, with the bilinear form `a` and linear form `L`. +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 + +name_read = None +if problem is ProblemType.DIRICHLET: + name_read = "Temperature" + read_mesh = "Dirichlet-Mesh" +else: + name_read = "Heat-Flux" + read_mesh = "Neumann-Mesh" + +coupling_function = fem.Function(V_coup) + +if problem is ProblemType.DIRICHLET: + # modify Dirichlet boundary condition on coupling interface + bc_coup = fem.dirichletbc(coupling_function, dofs_coupling) + bcs.append(bc_coup) +if problem is ProblemType.NEUMANN: + # modify Neumann boundary condition on coupling interface, modify weak + # form correspondingly + n = ufl.FacetNormal(domain) + F -= dt * ufl.inner(coupling_function, n) * v * ufl.ds + coupling_function.interpolate(gradient_solver.compute(u_D)) +a = fem.form(ufl.lhs(F)) +L = fem.form(ufl.rhs(F)) + +# ## 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. +A = assemble_matrix(a, bcs=bcs) +A.assemble() +b = create_vector(fem.extract_function_spaces(L)) +uh = fem.Function(V) + +# ## 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. +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 = fem.Function(V_g) + +# boundaries point as always to the end of the timestep +u_exact.t += dt +u_D.interpolate(u_exact) + +f_err = fem.Function(V) +# create writer for output files +vtxwriter = io.VTXWriter(MPI.COMM_WORLD, f"output-{problem.name}.bp".lower(), [f_err]) +vtxwriter.write(t) + +if problem is ProblemType.NEUMANN: + mappingToOriginalSpace = None + dofs_coupling_coordinates = None + Tmp, mappingToOriginalSpace = V_coup.sub(0).collapse() + dofs_subspace = fem.locate_dofs_geometrical(Tmp, coupling_boundary) + # to get the vector dof's (true) coordinates, it is sufficient to use the dof's coordinates of an arbitrary subspace + dofs_coupling_coordinates = V_coup.tabulate_dof_coordinates()[dofs_subspace] + +while precice.is_coupling_ongoing(): + + if precice.requires_writing_checkpoint(): + precice.store_checkpoint(u_n, t, 0) + + precice_dt = precice.get_max_time_step_size() + dt = np.min([fenics_dt, precice_dt]) + + # Update the coupling expression with the new read data + precice.read_data(read_mesh, name_read, dt, coupling_function) + + # Update the right hand side reusing the initial vector + with b.localForm() as loc_b: + loc_b.set(0) + assemble_vector(b, L) + + # Apply Dirichlet boundary condition to the vector (according to the tutorial, the lifting operation is used to preserve the symmetry of the matrix) + # Boundary condition bc should be updated by u_D.interpolate above, since + # this function is wrapped into the bc object + apply_lifting(b, [a], [bcs]) + b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE) + set_bc(b, bcs) + + # Solve linear problem + solver.solve(b, uh.x.petsc_vec) + uh.x.scatter_forward() + + # 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 = gradient_solver.compute(uh) + flux.x.scatter_forward() + precice.write_data(coupling_mesh.get_name(), "Heat-Flux", flux) + elif problem is ProblemType.NEUMANN: + # Neumann problem reads flux and writes temperature on boundary to Dirichlet problem + precice.write_data(coupling_mesh.get_name(), "Temperature", uh) + + precice.advance(dt) + precice_dt = precice.get_max_time_step_size() + + # roll back to checkpoint + if precice.requires_reading_checkpoint(): + u_cp, t_cp, _ = precice.retrieve_checkpoint() + u_n.x.array[:] = u_cp.x.array + t = t_cp + else: # update solution + # Update solution at previous time step (u_n) + u_n.x.array[:] = uh.x.array + f_err.x.array[:] = np.abs(u_n.x.array - u_D.x.array) + t += float(dt) + vtxwriter.write(t) + + if precice.is_time_window_complete(): + u_ref = fem.Function(V) + u_ref.interpolate(u_D) + error, error_pointwise = compute_errors(u_n, u_ref, 1e-4) + 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)! + +precice.finalize() + +vtxwriter.close() diff --git a/partitioned-heat-conduction-3d/solver-fenicsx/my_enums.py b/partitioned-heat-conduction-3d/solver-fenicsx/my_enums.py new file mode 100644 index 000000000..666ae348e --- /dev/null +++ b/partitioned-heat-conduction-3d/solver-fenicsx/my_enums.py @@ -0,0 +1,17 @@ +from enum import Enum + + +class ProblemType(Enum): + """ + Enum defines problem type. Details see above. + """ + DIRICHLET = "Dirichlet" # Dirichlet problem + NEUMANN = "Neumann" # Neumann problem + + +class DomainPart(Enum): + """ + Enum defines which part of the domain [x_left, x_right] x [y_bottom, y_top] we compute. + """ + OUTER = 1 # left part of domain in simple interface case + INNER = 2 # right part of domain in simple interface case diff --git a/partitioned-heat-conduction-3d/solver-fenicsx/problem_setup.py b/partitioned-heat-conduction-3d/solver-fenicsx/problem_setup.py new file mode 100644 index 000000000..35719f452 --- /dev/null +++ b/partitioned-heat-conduction-3d/solver-fenicsx/problem_setup.py @@ -0,0 +1,131 @@ +""" +Problem setup for partitioned-heat-conduction/fenicsx tutorial +""" +from dolfinx.io import gmsh as gmshio +from dolfinx.mesh import DiagonalType, create_box +import dolfinx.mesh +from my_enums import DomainPart +import numpy as np +from mpi4py import MPI +import gmsh + + +def get_geometry(domain_part, communicator): + gmsh.initialize() + gmsh.model.add("HeatMesh") + + if domain_part is DomainPart.OUTER: + # In this section the outer box of the domain is defined. + # It is a cuboid from (0,0,0) to (2,1,2) and has a cuboid cutout from (1, .25, .5) to (2, .75, 1.5) + # The coupling interface is the surface from the cutout. + + # Definition of the coupling boundary + def coupling_bc(x): + tol = 1E-14 + top = np.logical_and.reduce(( + np.isclose(x[1], 0.75, atol=tol), + np.logical_or(x[0] >= 1, np.isclose(x[0], 1, tol)), + np.logical_or(x[2] >= 0.5, np.isclose(x[2], 0.5, tol)), + np.logical_or(x[2] <= 1.5, np.isclose(x[2], 1.5, tol)) + )) + + bot = np.logical_and.reduce(( + np.isclose(x[1], 0.25, tol), + np.logical_or(x[0] >= 1, np.isclose(x[0], 1, tol)), + np.logical_or(x[2] >= 0.5, np.isclose(x[2], 0.5, tol)), + np.logical_or(x[2] <= 1.5, np.isclose(x[2], 1.5, tol)) + )) + + sideCenter = np.logical_and.reduce(( + np.isclose(x[0], 1, tol), + x[1] >= 0.25, + x[1] <= 0.75, + x[2] >= 0.5, + x[2] <= 1.5 + )) + + sideLeft = np.logical_and.reduce(( + x[0] >= 1, + x[1] >= 0.25, + x[1] <= 0.75, + np.isclose(x[2], 0.5, tol) + )) + + sideRight = np.logical_and.reduce(( + x[0] >= 1, + x[1] >= 0.25, + x[1] <= 0.75, + np.isclose(x[2], 1.5, tol) + )) + + return np.logical_or.reduce((top, bot, sideCenter, sideLeft, sideRight)) + + # definition of the boundary excluding the coupling boundary + + def boundary_bc(x): + tol = 1E-14 + or_part = np.logical_or.reduce(( + np.isclose(x[0], 0, tol), + np.isclose(x[0], 2, tol), + np.isclose(x[1], 0, tol), + np.isclose(x[1], 1, tol), + np.isclose(x[2], 0, tol), + np.isclose(x[2], 2, tol) + )) + and_part = np.logical_and.reduce(( + np.isclose(x[0], 2, tol), + x[1] >= 0.25, + x[1] <= 0.75, + x[2] >= 0.5, + x[2] <= 1.5, + + )) + return np.logical_and(or_part, ~and_part) + + # for creating the mesh, gmsh is used + # create full cuboid + gmsh.model.occ.addBox(0, 0, 0, 2, 1, 2, 1) + gmsh.model.occ.synchronize() + # create cuboid for cutout + gmsh.model.occ.addBox(1, 0.25, .5, 1, .5, 1, 2) + gmsh.model.occ.synchronize() + + # remove smaller cuboid from large cuboid + _, _ = gmsh.model.occ.cut([(3, 1)], [(3, 2)], 3, True, True) + gmsh.model.occ.synchronize() + gmsh.model.addPhysicalGroup(3, [3], 4, "foo") + # generate mesh with max. edge length of 0.15 + gmsh.model.mesh.setOrder(2) + gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 0.15) + gmsh.model.mesh.generate(3) + # convert from gmsh to a dolfinx representation + outer_mesh = gmshio.model_to_mesh(gmsh.model, communicator, 0, 3).mesh + gmsh.finalize() + return outer_mesh, coupling_bc, boundary_bc + else: + # inner part of the domain + # it is defined by a cuboid from (1, .25, .5) to (2, .75, 1.5) + + # 5 sides of the cuboid define the coupling interface + def coupling_bc(x): + tol = 1E-14 + return np.logical_or( + np.logical_or(np.isclose(x[1], 0.75, tol), np.isclose(x[1], 0.25, tol)), + np.logical_or(np.isclose(x[0], 1, tol), + np.logical_or(np.isclose(x[2], 0.5, tol), np.isclose(x[2], 1.5, tol)))) + + def boundary_bc(x): + tol = 1E-14 + return np.isclose(x[0], 2, tol) + + # create cuboid + gmsh.model.occ.addBox(1, 0.25, .5, 1, .5, 1, 1) + gmsh.model.occ.synchronize() + gmsh.model.addPhysicalGroup(3, [1], 2, "foo") + # generate mesh + gmsh.model.mesh.setOrder(2) + gmsh.model.mesh.generate(3) + # convert from gmsh to dolfinx representation + inner_mesh = gmshio.model_to_mesh(gmsh.model, communicator, 0, 3).mesh + gmsh.finalize() + return inner_mesh, coupling_bc, boundary_bc diff --git a/partitioned-heat-conduction-3d/solver-fenicsx/requirements.txt b/partitioned-heat-conduction-3d/solver-fenicsx/requirements.txt new file mode 100644 index 000000000..f8c5b4b8f --- /dev/null +++ b/partitioned-heat-conduction-3d/solver-fenicsx/requirements.txt @@ -0,0 +1,4 @@ +numpy >1, <2 +fenicsxprecice +scipy +sympy diff --git a/partitioned-heat-conduction-3d/solver-fenicsx/run.sh b/partitioned-heat-conduction-3d/solver-fenicsx/run.sh new file mode 100755 index 000000000..297b728ed --- /dev/null +++ b/partitioned-heat-conduction-3d/solver-fenicsx/run.sh @@ -0,0 +1,21 @@ +#!/bin/sh +set -e -u + +python3 -m venv --system-site-packages .venv +. .venv/bin/activate +pip install ../../.. +pip install -r requirements.txt + +while getopts ":dn" opt; do + case ${opt} in + d) + python3 heat.py Dirichlet --error-tol 10e-3 + ;; + n) + python3 heat.py Neumann --error-tol 10e-3 + ;; + \?) + echo "Usage: cmd [-d] [-n]" + ;; + esac +done