diff --git a/.gitignore b/.gitignore index 927c110d4..3f91e9fa5 100644 --- a/.gitignore +++ b/.gitignore @@ -13,6 +13,7 @@ out *.vtk *.vtu *.pvd +*.series # preCICE precice-profiling/ diff --git a/changelog-entries/.keep b/changelog-entries/.keep new file mode 100644 index 000000000..8b1378917 --- /dev/null +++ b/changelog-entries/.keep @@ -0,0 +1 @@ + diff --git a/changelog-entries/599.md b/changelog-entries/599.md new file mode 100644 index 000000000..a4458a609 --- /dev/null +++ b/changelog-entries/599.md @@ -0,0 +1,2 @@ +- Updated the channel-transport tutorial to use a parallel-explicit coupling scheme with remeshing support. +- Added AMR support to the transport solver of the channel-transport case. diff --git a/channel-transport/README.md b/channel-transport/README.md index 267a70550..3f597ac24 100644 --- a/channel-transport/README.md +++ b/channel-transport/README.md @@ -1,7 +1,7 @@ --- title: Channel transport permalink: tutorials-channel-transport.html -keywords: volume coupling, chemistry, OpenFOAM, Nutils, species, transport +keywords: volume coupling, chemistry, OpenFOAM, Nutils, species, transport, remeshing summary: A CFD problem is coupled to a transport (of, e.g., a chemistry species) in a uni-directional way. --- @@ -22,6 +22,10 @@ The behavior of the blob over the full 200 timesteps looks as follows: Animation of blob over 200 timesteps. +The transport solver also supports the use of adaptive mesh refinement. + +![Setup with AMR](images/tutorials-channel-transport-amr.png) + ## Configuration preCICE configuration (image generated using the [precice-config-visualizer](https://precice.org/tooling-config-visualization.html)): @@ -39,10 +43,12 @@ Fluid participant: Transport participant: -* Nutils. For more information, have a look at the [Nutils adapter documentation](https://precice.org/adapter-nutils.html). This Nutils solver requires at least Nutils v7.0. +* Nutils with support for [adaptive mesh refinement](https://precice.org/couple-your-code-moving-or-changing-meshes.html#pseudo-reference-domain). For more information, have a look at the [Nutils adapter documentation](https://precice.org/adapter-nutils.html). This Nutils solver requires at least Nutils v7.0 and a preCICE release with [remeshing support](couple-your-code-moving-or-changing-meshes.html#remeshing-using-precice). ## Running the simulation +For the fluid solver, use Nutils for ease of installation and OpenFOAM for speed. + Open two separate terminals and start one fluid and one transport participant by calling the respective run scripts `run.sh` located in each of the participants' directory. For example: ```bash @@ -50,13 +56,27 @@ cd fluid-nutils ./run.sh ``` -and +and either the non-adaptive mesh transport solver ```bash cd transport-nutils ./run.sh ``` +or the adaptive mesh transport solver + +```bash +cd transport-nutils +./run.sh remesh=True +``` + +The mesh refinement frequency can also be adjusted: + +```bash +cd transport-nutils +./run.sh remesh=True frequency=4 # refine every 4 time-windows +``` + ## Post-processing All solvers generate `vtk` files which can be visualized using, e.g., ParaView. diff --git a/channel-transport/fluid-nutils/requirements.txt b/channel-transport/fluid-nutils/requirements.txt index 78eb2a14a..56c3d4a92 100644 --- a/channel-transport/fluid-nutils/requirements.txt +++ b/channel-transport/fluid-nutils/requirements.txt @@ -2,3 +2,4 @@ setuptools # required by nutils nutils==7 numpy >1, <2 pyprecice~=3.0 +setuptools diff --git a/channel-transport/images/tutorials-channel-transport-amr.png b/channel-transport/images/tutorials-channel-transport-amr.png new file mode 100644 index 000000000..21c16619a Binary files /dev/null and b/channel-transport/images/tutorials-channel-transport-amr.png differ diff --git a/channel-transport/images/tutorials-channel-transport-precice-config.png b/channel-transport/images/tutorials-channel-transport-precice-config.png index 3fd022358..f8c69ec51 100644 Binary files a/channel-transport/images/tutorials-channel-transport-precice-config.png and b/channel-transport/images/tutorials-channel-transport-precice-config.png differ diff --git a/channel-transport/precice-config.xml b/channel-transport/precice-config.xml index 9abf595ca..902bde439 100644 --- a/channel-transport/precice-config.xml +++ b/channel-transport/precice-config.xml @@ -1,5 +1,5 @@ - + - + - + diff --git a/channel-transport/transport-nutils/requirements.txt b/channel-transport/transport-nutils/requirements.txt index 78eb2a14a..65ebe15cd 100644 --- a/channel-transport/transport-nutils/requirements.txt +++ b/channel-transport/transport-nutils/requirements.txt @@ -1,4 +1,6 @@ setuptools # required by nutils nutils==7 numpy >1, <2 -pyprecice~=3.0 +# remeshing support will be released with preCICE version 3.2 +pyprecice @ git+https://github.com/precice/python-bindings.git@develop +setuptools diff --git a/channel-transport/transport-nutils/run.sh b/channel-transport/transport-nutils/run.sh index ab44dadce..a596de3c2 100755 --- a/channel-transport/transport-nutils/run.sh +++ b/channel-transport/transport-nutils/run.sh @@ -7,6 +7,6 @@ exec > >(tee --append "$LOGFILE") 2>&1 python3 -m venv .venv . .venv/bin/activate pip install -r requirements.txt -python3 transport.py +python3 transport.py "$@" close_log diff --git a/channel-transport/transport-nutils/transport.py b/channel-transport/transport-nutils/transport.py index 83aaebd82..53c19a0c8 100644 --- a/channel-transport/transport-nutils/transport.py +++ b/channel-transport/transport-nutils/transport.py @@ -6,27 +6,92 @@ from nutils import function, mesh, cli, solver, export import treelog as log +import argparse import numpy as np +import json import precice from mpi4py import MPI -def main(): +def reinitialize_namespace(domain, geom): + # cloud of Gauss points + gauss = domain.sample("bezier", degree=2) + + # Nutils namespace + ns = function.Namespace(fallback_length=2) + ns.x = geom + ns.basis = domain.basis("h-std", degree=1) # linear finite elements + ns.u = "basis_n ?lhs_n" # solution + ns.projectedu = "basis_n ?projectedlhs_n" + ns.gradu = "u_,i" # gradient of lhstion + ns.dudt = "basis_n (?lhs_n - ?lhs0_n) / ?dt" # time derivative + ns.vbasis = gauss.basis() + ns.velocity_i = "vbasis_n ?velocity_ni" + ns.k = 0.1 # diffusivity + ns.xblob = 1, 1 + ns.uinit = ".5 - .5 tanh(((x_i - xblob_i) (x_i - xblob_i) - .5) / .1)" # blob + + # define the weak form + res = gauss.integral("(basis_n (dudt + (velocity_i u)_,i) + k basis_n,i u_,i) d:x" @ ns) + + # define Dirichlet boundary condition + sqr = domain.boundary["inflow"].integral("u^2 d:x" @ ns, degree=2) + cons = solver.optimize("lhs", sqr, droptol=1e-15) + + return ns, res, cons, gauss + + +def refine_mesh(limits, ns, domain_coarse, domain_nm1, lhs_nm1): + """ + At the time of the calling of this function a predicted lhstion exists in ns.phi + """ + # ----- Refine the coarse mesh according to the projected lhstion to get a predicted refined topology ---- + domain_ref = domain_coarse + for level, limit in enumerate(limits): + print("refinement level = {}".format(level)) + domain_union1 = domain_nm1 & domain_ref + smpl = domain_union1.sample('uniform', 5) + ielem, criterion = smpl.eval([domain_ref.f_index, function.sqrt(ns.gradu[0]**2 + ns.gradu[1]**2) > limit], + lhs=lhs_nm1) + + # Refine the elements for which at least one point tests true. + domain_ref = domain_ref.refined_by(np.unique(ielem[criterion])) + # ---------------------------------------------------------------------------------------------------- + + # Create a new projection mesh which is the union of the previous refined mesh and the predicted mesh + domain_union = domain_nm1 & domain_ref + + # ----- Project the lhstion of the last time step on the projection mesh ----- + ns.projectedu = function.dotarg('projectedlhs', domain_ref.basis('h-std', degree=1)) + sqru = domain_union.integral((ns.projectedu - ns.u) ** 2, degree=2) + lhs = solver.optimize('projectedlhs', sqru, droptol=1E-12, arguments=dict(lhs=lhs_nm1)) + + return domain_ref, lhs + + +def main(remesh: bool, frequency: int, limits: [float], visualize: bool): print("Running Nutils") - # define the Nutils mesh - nx = 120 - ny = 32 + if frequency > 0: + # Use a coarser grid for remeshing cases + nx = 60 + ny = 16 + else: + # define the Nutils mesh + nx = 120 + ny = 32 + step_start = nx // 3 step_end = nx // 2 - step_hight = ny // 2 + step_height = ny // 2 grid = np.linspace(0, 6, nx + 1), np.linspace(0, 2, ny + 1) domain, geom = mesh.rectilinear(grid) domain = domain.withboundary(inflow="left", outflow="right", wall="top,bottom") - domain[ - step_start:step_end, :step_hight + step_start:step_end, :step_height ].withboundary(wall="left,top,right") + domain_coarse = domain # Retain the original coarse domain for mesh refinement later on # cloud of Gauss points gauss = domain.sample("gauss", degree=2) @@ -36,6 +101,8 @@ def main(): ns.x = geom ns.basis = domain.basis("std", degree=1) # linear finite elements ns.u = "basis_n ?lhs_n" # solution + ns.projectedu = "basis_n ?projectedlhs_n" + ns.gradu = "u_,i" # gradient of lhstion ns.dudt = "basis_n (?lhs_n - ?lhs0_n) / ?dt" # time derivative ns.vbasis = gauss.basis() ns.velocity_i = "vbasis_n ?velocity_ni" @@ -50,6 +117,30 @@ def main(): sqr = domain.boundary["inflow"].integral("u^2 d:x" @ ns, degree=2) cons = solver.optimize("lhs", sqr, droptol=1e-15) + # set blob as initial condition + sqr = domain.integral("(u - uinit)^2" @ ns, degree=2) + lhs0 = solver.optimize("lhs", sqr) + + if remesh: + # Initial refinement according to initial condition + print("Performing initial mesh refinement") + for level, limit in enumerate(limits): + print("refinement level = {}".format(level)) + smpl = domain.sample('uniform', 5) + ielem, criterion = smpl.eval([domain.f_index, function.sqrt( + ns.gradu[0]**2 + ns.gradu[1]**2) > limit], lhs=lhs0) + + # Refine the elements for which at least one point tests true. + domain = domain.refined_by(np.unique(ielem[criterion])) + + ns, res, cons, gauss = reinitialize_namespace(domain, geom) + + # set blob as initial condition after each refinement + sqr = domain.integral("(u - uinit)^2" @ ns, degree=2) + lhs0 = solver.optimize("lhs", sqr) + + lhs = lhs0 + # preCICE setup participant = precice.Participant("Transport", "../precice-config.xml", 0, 1) @@ -63,25 +154,18 @@ def main(): participant.initialize() - timestep = 0 + timestep = 1 solver_dt = 0.005 precice_dt = participant.get_max_time_step_size() dt = min(precice_dt, solver_dt) - - # set blob as initial condition - sqr = domain.integral("(u - uinit)^2" @ ns, degree=2) - lhs0 = solver.optimize("lhs", sqr) + t = 0 # initialize the velocity values velocity_values = np.zeros_like(vertices) - while participant.is_coupling_ongoing(): + series = {"file-series-version": "1.0", "files": []} - if timestep % 1 == 0: # visualize - bezier = domain.sample("bezier", 2) - x, u = bezier.eval(["x_i", "u"] @ ns, lhs=lhs0) - with log.add(log.DataLog()): - export.vtk("Transport_" + str(timestep), bezier.tri, x, T=u) + while participant.is_coupling_ongoing(): precice_dt = participant.get_max_time_step_size() @@ -96,15 +180,40 @@ def main(): "lhs", res, constrain=cons, arguments=dict(lhs0=lhs0, dt=dt, velocity=velocity_values) ) + if remesh and timestep % frequency == 0: + domain, lhs = refine_mesh(limits, ns, domain_coarse, domain, lhs) + ns, res, cons, gauss = reinitialize_namespace(domain, geom) + + vertices = gauss.eval(ns.x) + participant.reset_mesh(mesh_name) # Throws away the entire mesh + vertex_ids = participant.set_mesh_vertices(mesh_name, vertices) # Redefine the mesh + + if visualize: + bezier = domain.sample("bezier", 2) + x, u = bezier.eval(["x_i", "u"] @ ns, lhs=lhs) + with log.add(log.DataLog()): + export.vtk(f"Transport_{timestep}", bezier.tri, x, T=u) + series["files"].append({"name": f"Transport_{timestep}.vtk", "time": t}) + # do the coupling participant.advance(dt) # advance variables + t += dt timestep += 1 lhs0 = lhs + if visualize: + with open("Transport.vtk.series", "w") as f: + json.dump(series, f) + participant.finalize() if __name__ == "__main__": - cli.run(main) + + def run(remesh: bool = False, frequency: int = 2, visualize: bool = True): + # we use 1 and 2 as default refinement limits. + main(remesh, frequency, [1.0, 2.0], visualize) + + cli.run(run)