Growing mesh case

.gitignore

@@ -33,6 +33,7 @@ build/
 .venv/
 __pycache__/
 *.pyc
+*.egg-info
 
 # Rust
 Cargo.lock

growing-mesh/A/run.sh

@@ -0,0 +1,18 @@
+#!/usr/bin/env bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+python3 -m venv .venv
+. .venv/bin/activate
+pip install ../solver-python
+pip freeze ../solver-python/ > ../solver-python/pip-installed-packages.log
+
+if [ $# -eq 0 ]; then
+  growing A
+else
+  mpirun -n "$@" growing A
+fi
+
+close_log

growing-mesh/B/run.sh

@@ -0,0 +1,18 @@
+#!/usr/bin/env bash
+set -e -u
+
+. ../../tools/log.sh
+exec > >(tee --append "$LOGFILE") 2>&1
+
+python3 -m venv .venv
+. .venv/bin/activate
+pip install ../solver-python
+pip freeze ../solver-python/ > ../solver-python/pip-installed-packages.log
+
+if [ $# -eq 0 ]; then
+  growing B
+else
+  mpirun -n "$@" growing B
+fi
+
+close_log

growing-mesh/README.md

@@ -0,0 +1,56 @@
+---
+title: Growing Mesh
+permalink: tutorials-growing-mesh.html
+keywords: python, remeshing
+summary: The growing mesh case is a showcase example of two solvers which grow their mesh at predefined points in time.
+---
+
+{% note %}
+Get the [case files of this tutorial](https://github.com/precice/tutorials/tree/master/growing-mesh). Read how in the [tutorials introduction](https://precice.org/tutorials.html).
+{% endnote %}
+
+## Motivation
+
+This case is an extremely simplified version of the scenario presented by Jean-Marc Gratien during Coupled 2023 *Implementing a Comprehensive Hydromechanical Model for Sedimentary Basins by Coupling a 3D Mechanical Code to a Classic Basin Fluid Flow Code with the PreCICE Library*.
+
+In the original scenario a hydrogeological and a geomechanical model are coupled via preCICE using an implicit coupling schemes.
+The scheme iterates until both model agree on effective stress.
+At fixed points in time, geological events deposit new layers onto the mesh, effectively growing the mesh into z-direction.
+
+This tutorial aims to implement such a growing mesh scenario without the implicit coupling.
+
+## Setup
+
+The problem consists of a unit-square uniformly discretized by 512 x 512 nodes.
+The unit square is partitioned among the ranks of a solver according to an n x m grid, where 512 must be divisible by n and m.
+The mesh starts with 2 nodes in z direction and at a given frequency, 2 nodes are added to the mesh, changing only the load per rank, not the partitioning.
+
+## Configuration
+
+preCICE configuration (image generated using the [precice-config-visualizer](https://precice.org/tooling-config-visualization.html)):
+
+![preCICE configuration visualization](images/tutorials-growing-mesh-precice-config.png)
+
+## Available solvers
+
+There is a generic python solver that can be used for either A or B.
+
+The runs scripts in participant folders `A` and `B` accept the amount of ranks as an optional parameter.
+
+## Running the Simulation
+
+Pass the amount of ranks to the run script of the solvers.
+Not passing a number, runs the simulation on a single rank.
+To run both on a two rank each, use:
+
+```bash
+cd A
+./run.sh 2
+```
+
+and
+
+```bash
+cd B
+./run.sh 2
+```

growing-mesh/precice-config.xml

@@ -0,0 +1,50 @@
+<?xml version="1.0" encoding="UTF-8" ?>
+<precice-configuration experimental="True" allow-remeshing="True">
+  <log>
+    <sink
+      filter="%Severity% >= debug and %Rank% = 0"
+      format="---[precice] %ColorizedSeverity% %Message%"
+      enabled="true" />
+  </log>
+
+  <profiling mode="all" />
+
+  <data:scalar name="Data-A" />
+  <data:scalar name="Data-B" />
+
+  <mesh name="A-Mesh" dimensions="3">
+    <use-data name="Data-A" />
+    <use-data name="Data-B" />
+  </mesh>
+
+  <mesh name="B-Mesh" dimensions="3">
+    <use-data name="Data-A" />
+    <use-data name="Data-B" />
+  </mesh>
+
+  <participant name="A">
+    <provide-mesh name="A-Mesh" />
+    <write-data name="Data-A" mesh="A-Mesh" />
+    <read-data name="Data-B" mesh="A-Mesh" />
+    <receive-mesh name="B-Mesh" from="B" />
+    <mapping:nearest-neighbor direction="read" from="B-Mesh" to="A-Mesh" constraint="consistent" />
+  </participant>
+
+  <participant name="B">
+    <provide-mesh name="B-Mesh" />
+    <read-data name="Data-A" mesh="B-Mesh" />
+    <write-data name="Data-B" mesh="B-Mesh" />
+    <receive-mesh name="A-Mesh" from="A" />
+    <mapping:nearest-neighbor direction="read" from="A-Mesh" to="B-Mesh" constraint="consistent" />
+  </participant>
+
+  <m2n:sockets acceptor="A" connector="B" exchange-directory=".." />
+
+  <coupling-scheme:parallel-explicit>
+    <time-window-size value="1" />
+    <max-time-windows value="12" />
+    <participants first="A" second="B" />
+    <exchange data="Data-A" mesh="A-Mesh" from="A" to="B" />
+    <exchange data="Data-B" mesh="B-Mesh" from="B" to="A" />
+  </coupling-scheme:parallel-explicit>
+</precice-configuration>

growing-mesh/solver-python/pyproject.toml

@@ -0,0 +1,11 @@
+[project]
+name = "growing"
+version = "0"
+dependencies = [
+    "numpy",
+    "pyprecice @ git+https://github.com/precice/python-bindings.git@develop",
+    "mpi4py"
+]
+
+[project.scripts]
+growing = "solver:main"

growing-mesh/solver-python/solver.py

@@ -0,0 +1,137 @@
+#!/bin/python3
+
+import precice
+import numpy as np
+import math
+import sys
+from mpi4py import MPI
+
+import argparse
+
+
+def split(num):
+    for a in range(math.isqrt(num), 0, -1):
+        if num % a == 0:
+            return a, num // a
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("participant", choices=["A", "B"])
+    parser.add_argument("--config", "-c", default="../precice-config.xml")
+    parser.add_argument("--no-remesh", dest="remesh", action="store_false")
+    parser.add_argument("--nx", type=int, default=512)
+    parser.add_argument("--ny", type=int, default=512)
+    args = parser.parse_args()
+
+    participant_name = args.participant
+    remote_name = "A" if participant_name == "B" else "B"
+
+    # x is partitioned per rank and doesn't change
+    nx = args.nx
+    ny = args.ny
+    x = 0.0, 1.0
+    y = 0.0, 1.0
+
+    # y grows over time
+    newNodesPerEvent = 2
+    eventFrequency = 3  # time windows
+    dz = 0.1
+
+    # Handle partitioning
+    world = MPI.COMM_WORLD
+    size: int = world.size
+    rank: int = world.rank
+    xr, yr = split(size)
+
+    assert nx % xr == 0, f"Cannot split {nx=} by {xr=} ranks"
+    assert ny % yr == 0, f"Cannot split {ny=} by {yr=} ranks"
+
+    pnx = nx // xr
+    pny = ny // yr
+
+    px = rank % xr
+    py = rank // xr
+    print(f"Rank {rank}/{size} has partition ({px}, {py})/({xr}, {yr})")
+    if rank == 0:
+        print(
+            f"Each of {size} partitions has node size {pnx}x{pny} = {pnx * pny}"
+            "for a total of {nx * ny} nodes on the base"
+        )
+
+    def getMesh(nz):
+        basex = np.linspace(0, 1, nx)[px * pnx: (px + 1) * pnx]
+        basey = np.linspace(0, 1, ny)[py * pny: (py + 1) * pny]
+        z = np.array(range(nz)) * dz
+        return np.stack(np.meshgrid(basex, basey, z, indexing="ij"), axis=-1).reshape(
+            -1, 3
+        )
+
+    def requiresEvent(tw):
+        return tw % eventFrequency == 0
+
+    assert not requiresEvent(eventFrequency - 1)
+    assert requiresEvent(eventFrequency)
+    assert not requiresEvent(eventFrequency + 1)
+
+    def eventsAt(tw):
+        # First event block at tw=0, second at eventFrequency
+        return 1 + math.floor(tw / eventFrequency)
+
+    assert eventsAt(0) == 1
+    assert eventsAt(eventFrequency - 1) == 1
+    assert eventsAt(eventFrequency) == 2
+    assert eventsAt(eventFrequency + 1) == 2
+
+    def getMeshAtTimeWindow(tw):
+        znodes = eventsAt(tw) * newNodesPerEvent
+        return getMesh(znodes)
+
+    def dataAtTimeWindow(tw):
+        znodes = eventsAt(tw) * newNodesPerEvent
+        total = pnx * pny * znodes
+        return np.full(total, tw)
+
+    participant = precice.Participant(participant_name, args.config, rank, size)
+
+    mesh_name = participant_name + "-Mesh"
+    read_data_name = "Data-" + remote_name
+    write_data_name = "Data-" + participant_name
+
+    coords = getMeshAtTimeWindow(0)
+    vertex_ids = participant.set_mesh_vertices(mesh_name, coords)
+    participant.initialize()
+
+    tw = 1
+    while participant.is_coupling_ongoing():
+        dt = participant.get_max_time_step_size()
+
+        data = participant.read_data(mesh_name, read_data_name, vertex_ids, dt)
+        if rank == 0:
+            print(f"{participant_name}: data: {data[0]} * {len(data)}")
+
+        if args.remesh and requiresEvent(tw):
+            oldCount = len(coords)
+            coords = getMeshAtTimeWindow(tw)
+            if rank == 0:
+                print(
+                    f"{participant_name}: Event grows local mesh from {oldCount} to {len(coords)}"
+                    "and global mesh from {oldCount * size} to {len(coords) * size}"
+                )
+            participant.reset_mesh(mesh_name)
+            vertex_ids = participant.set_mesh_vertices(mesh_name, coords)
+
+        participant.write_data(
+            mesh_name, write_data_name, vertex_ids, dataAtTimeWindow(tw)
+        )
+
+        participant.advance(dt)
+        tw += 1
+
+
+if __name__ == "__main__":
+    try:
+        main()
+    except Exception as e:
+        print(e)
+        sys.exit(1)