Skip to content

Commit

Permalink
Port the dunefem solver for flow-over-heated-plate simulation for pre…
Browse files Browse the repository at this point in the history
…CICE v3 (#470)

* port the dunefem solver for flow-over-heated-plate for preCICE v3

* Add a comment for the mesh creation

* Formatting

* Formatting

* Use venv for dune-fem case

---------

Co-authored-by: Jun Chen <jun.chen@ipvs.uni-stuttgart.de>
Co-authored-by: Benjamin Uekermann <benjamin.uekermann@ipvs.uni-stuttgart.de>
  • Loading branch information
3 people authored Mar 21, 2024
1 parent 3f139d3 commit 383e824
Show file tree
Hide file tree
Showing 4 changed files with 35 additions and 26 deletions.
2 changes: 1 addition & 1 deletion flow-over-heated-plate/README.md
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ Solid participant:

* Nutils. For more information, have a look at the [Nutils adapter documentation](https://precice.org/adapter-nutils.html).

* Dune-Fem. For more information, have a look at the [official documentation of Dune-Fem](https://www.dune-project.org/sphinx/dune-fem/). The solver can be installed through [PyPI](https://pypi.org/project/dune-fem/). Make sure that you are in a Python virtual environment first, which you can create inside the `solid-dune` directory and load again before running (you may need to install some tools again in this environment). Please note that Dune-Fem uses just-in-time compilation: The first time you run the solver script, it will take some time.
* Dune-Fem. For more information, have a look at the [official documentation of Dune-Fem](https://www.dune-project.org/sphinx/dune-fem/). The `run.sh` script installs the solver from [PyPI](https://pypi.org/project/dune-fem/) into a Python virtual environment. Please note that Dune-Fem uses just-in-time compilation: The first time you run the solver script, it will take some time.

It is also possible to use CalculiX as solid solver. In that case, two coupling meshes are needed: CalculiX read/writes temperatures on nodes, but read/writes heat-fluxes on face centers. This requires some adaptation of the `precice-config.xml` file, and [a separate tutorial](tutorials-flow-over-heated-plate-two-meshes.html) has been designed for it.

Expand Down
2 changes: 2 additions & 0 deletions flow-over-heated-plate/solid-dunefem/requirements.txt
Original file line number Diff line number Diff line change
@@ -0,0 +1,2 @@
dune-fem>=2.8
pyprecice==3
4 changes: 3 additions & 1 deletion flow-over-heated-plate/solid-dunefem/run.sh
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
#!/bin/bash
set -e -u

python3 -m venv .venv
. .venv/bin/activate
pip install -r requirements.txt
python3 solid.py

53 changes: 29 additions & 24 deletions flow-over-heated-plate/solid-dunefem/solid.py
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
from dune.fem.operator import galerkin
from dune.fem.utility import Sampler
from dune.grid import cartesianDomain
from dune.alugrid import aluSimplexGrid
from dune.alugrid import aluGrid
from dune.fem.function import uflFunction, gridFunction
from dune.ufl import expression2GF

Expand All @@ -32,19 +32,18 @@
x_right = x_left + 1

# preCICE setup
interface = precice.Interface("Solid", "../precice-config.xml", 0, 1)
participant = precice.Participant("Solid", "../precice-config.xml", 0, 1)

# define coupling mesh
mesh_name = "Solid-Mesh"
mesh_id = interface.get_mesh_id(mesh_name)
interface_x_coordinates = np.linspace(0, 1, nx + 1) # meshpoints for interface values
interface_x_coordinates = np.linspace(
0, 1, nx + 1) # meshpoints for interface values
vertices = [[x0, 0] for x0 in interface_x_coordinates]
vertex_ids = interface.set_mesh_vertices(mesh_id, vertices)
temperature_id = interface.get_data_id("Temperature", mesh_id)
flux_id = interface.get_data_id("Heat-Flux", mesh_id)

vertex_ids = participant.set_mesh_vertices(mesh_name, vertices)

domain = cartesianDomain([x_left, y_bottom], [x_right, y_top], [nx, ny])
mesh = aluSimplexGrid(domain, serial=True)
mesh = aluGrid(domain, 2, 2, elementType="simplex") # create a simple mesh with dimGrid=2 and dimWorld=2
space = solutionSpace(mesh, order=1)
u = ufl.TrialFunction(space)
v = ufl.TestFunction(space)
Expand Down Expand Up @@ -80,32 +79,36 @@ def u_gamma(xg):
scheme = solutionScheme([A == b, *bcs], solver='cg')

# Weak form of the flux
flux_expr = -(u - uold) * v / dt * ufl.dx - k * ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
flux_expr = -(u - uold) * v / dt * ufl.dx - k * \
ufl.dot(ufl.grad(u), ufl.grad(v)) * ufl.dx
flux_expr_operator = galerkin(flux_expr)
flux_sol = space.interpolate(u0, name='flux_sol')

# Sample the flux on the top edge
flux_sol_expr = expression2GF(flux_sol.space.grid, flux_sol, flux_sol.space.order)
flux_sol_expr = expression2GF(
flux_sol.space.grid, flux_sol, flux_sol.space.order)
sampler_weak_flux = Sampler(flux_sol_expr)
def flux_f_weak(): return sampler_weak_flux.lineSample([0., 0.], [1., 0.], nx + 1)[1] * nx

def flux_f_weak(): return sampler_weak_flux.lineSample(
[0., 0.], [1., 0.], nx + 1)[1] * nx

if not os.path.exists("output"):
os.makedirs("output")
vtk = mesh.sequencedVTK("output/heat", pointdata=[unew])

precice_dt = interface.initialize()
participant.initialize()
precice_dt = participant.get_max_time_step_size()
dt.assign(min(float(dt), precice_dt))

t = float(dt)

while interface.is_coupling_ongoing():
while participant.is_coupling_ongoing():

if interface.is_action_required(precice.action_write_iteration_checkpoint()): # write checkpoint
# write checkpoint
if participant.requires_writing_checkpoint():
t_check = t
ucheckpoint.assign(uold)
interface.mark_action_fulfilled(precice.action_write_iteration_checkpoint())

ug = interface.read_block_scalar_data(temperature_id, vertex_ids)
ug = participant.read_data(mesh_name, "Temperature", vertex_ids, dt)
ug_interp.y = ug

scheme.model.dt = dt
Expand All @@ -115,25 +118,27 @@ def flux_f_weak(): return sampler_weak_flux.lineSample([0., 0.], [1., 0.], nx +
flux_expr_operator(unew, flux_sol)
flux_values = flux_f_weak() # sample the flux function

interface.write_block_scalar_data(flux_id, vertex_ids, flux_values)
participant.write_data(mesh_name, "Heat-Flux", vertex_ids, flux_values)

precice_dt = interface.advance(dt)
participant.advance(dt)
precice_dt = participant.get_max_time_step_size()
dt.assign(min(float(dt), precice_dt))

if interface.is_action_required(precice.action_read_iteration_checkpoint()): # roll back to checkpoint
# roll back to checkpoint
if participant.requires_reading_checkpoint():
t = t_check
uold.assign(ucheckpoint)
interface.mark_action_fulfilled(precice.action_read_iteration_checkpoint())

else: # update solution

uold.assign(unew)
t += float(dt)

if interface.is_time_window_complete():
tol = 10e-5 # we need some tolerance, since otherwise output might be skipped.
if participant.is_time_window_complete():
# we need some tolerance, since otherwise output might be skipped.
tol = 10e-5
if abs((t + tol) % dt_out) < 2 * tol: # output if t is a multiple of dt_out
print("output vtk for time = {}".format(float(t)))
vtk()

interface.finalize()
participant.finalize()

0 comments on commit 383e824

Please sign in to comment.