From 383e824fa924a44aed8e5c95d00fa79a94999848 Mon Sep 17 00:00:00 2001 From: June <94080048+Fujikawas@users.noreply.github.com> Date: Thu, 21 Mar 2024 09:48:32 +0100 Subject: [PATCH] Port the dunefem solver for flow-over-heated-plate simulation for preCICE 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 Co-authored-by: Benjamin Uekermann --- flow-over-heated-plate/README.md | 2 +- .../solid-dunefem/requirements.txt | 2 + flow-over-heated-plate/solid-dunefem/run.sh | 4 +- flow-over-heated-plate/solid-dunefem/solid.py | 53 ++++++++++--------- 4 files changed, 35 insertions(+), 26 deletions(-) create mode 100644 flow-over-heated-plate/solid-dunefem/requirements.txt diff --git a/flow-over-heated-plate/README.md b/flow-over-heated-plate/README.md index 8d651fbec..a83b6abc9 100644 --- a/flow-over-heated-plate/README.md +++ b/flow-over-heated-plate/README.md @@ -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. diff --git a/flow-over-heated-plate/solid-dunefem/requirements.txt b/flow-over-heated-plate/solid-dunefem/requirements.txt new file mode 100644 index 000000000..b0b7bdf42 --- /dev/null +++ b/flow-over-heated-plate/solid-dunefem/requirements.txt @@ -0,0 +1,2 @@ +dune-fem>=2.8 +pyprecice==3 diff --git a/flow-over-heated-plate/solid-dunefem/run.sh b/flow-over-heated-plate/solid-dunefem/run.sh index 21bd5c3b1..8678a56d1 100755 --- a/flow-over-heated-plate/solid-dunefem/run.sh +++ b/flow-over-heated-plate/solid-dunefem/run.sh @@ -1,5 +1,7 @@ #!/bin/bash set -e -u +python3 -m venv .venv +. .venv/bin/activate +pip install -r requirements.txt python3 solid.py - diff --git a/flow-over-heated-plate/solid-dunefem/solid.py b/flow-over-heated-plate/solid-dunefem/solid.py index 24ee12c16..fc56c4e5b 100644 --- a/flow-over-heated-plate/solid-dunefem/solid.py +++ b/flow-over-heated-plate/solid-dunefem/solid.py @@ -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 @@ -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) @@ -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 @@ -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()