diff --git a/.gitignore b/.gitignore index ac8ae3c1f..9526bada0 100644 --- a/.gitignore +++ b/.gitignore @@ -33,6 +33,7 @@ build/ .venv/ __pycache__/ *.pyc +*.egg-info # Rust Cargo.lock diff --git a/growing-mesh/A/clean.sh b/growing-mesh/A/clean.sh new file mode 100755 index 000000000..494c80414 --- /dev/null +++ b/growing-mesh/A/clean.sh @@ -0,0 +1,7 @@ +#!/usr/bin/env sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . diff --git a/growing-mesh/A/run.sh b/growing-mesh/A/run.sh new file mode 100755 index 000000000..d5ea61494 --- /dev/null +++ b/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 \ No newline at end of file diff --git a/growing-mesh/B/clean.sh b/growing-mesh/B/clean.sh new file mode 100755 index 000000000..494c80414 --- /dev/null +++ b/growing-mesh/B/clean.sh @@ -0,0 +1,7 @@ +#!/usr/bin/env sh +set -e -u + +. ../../tools/cleaning-tools.sh + +clean_precice_logs . +clean_case_logs . diff --git a/growing-mesh/B/run.sh b/growing-mesh/B/run.sh new file mode 100755 index 000000000..fb3fbe2bd --- /dev/null +++ b/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 \ No newline at end of file diff --git a/growing-mesh/README.md b/growing-mesh/README.md new file mode 100644 index 000000000..b4f842dbb --- /dev/null +++ b/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 +``` diff --git a/growing-mesh/clean-tutorial.sh b/growing-mesh/clean-tutorial.sh new file mode 120000 index 000000000..4713f5092 --- /dev/null +++ b/growing-mesh/clean-tutorial.sh @@ -0,0 +1 @@ +../tools/clean-tutorial-base.sh \ No newline at end of file diff --git a/growing-mesh/images/tutorials-growing-mesh-precice-config.png b/growing-mesh/images/tutorials-growing-mesh-precice-config.png new file mode 100644 index 000000000..f46f4ab57 Binary files /dev/null and b/growing-mesh/images/tutorials-growing-mesh-precice-config.png differ diff --git a/growing-mesh/precice-config.xml b/growing-mesh/precice-config.xml new file mode 100644 index 000000000..1ffc743e8 --- /dev/null +++ b/growing-mesh/precice-config.xml @@ -0,0 +1,50 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/growing-mesh/solver-python/pyproject.toml b/growing-mesh/solver-python/pyproject.toml new file mode 100644 index 000000000..f3969b34f --- /dev/null +++ b/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" diff --git a/growing-mesh/solver-python/solver.py b/growing-mesh/solver-python/solver.py new file mode 100644 index 000000000..34f429d84 --- /dev/null +++ b/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} " + f"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)} " + f"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)