Skip to content

Commit 9f17050

Browse files
bangerthbangerth
authored
Merge pull request #24 from tjhei/lla-landlab2
Landlab: correctly support MPI on a subset
2 parents 0ae6a07 + a192272 commit 9f17050

File tree

6 files changed

+315
-31
lines changed

6 files changed

+315
-31
lines changed

‎include/aspect/mesh_deformation/landlab.h‎

Lines changed: 10 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -96,6 +96,16 @@ namespace aspect
9696
*/
9797
unsigned int n_landlab_ranks;
9898

99+
/**
100+
* The MPI communicator for the Landlab simulation.
101+
*/
102+
MPI_Comm landlab_communicator;
103+
104+
/**
105+
* Whether this rank is one of the ranks that runs the Landlab simulation.
106+
*/
107+
bool this_rank_runs_landlab;
108+
99109
/**
100110
* The path to the Landlab Python module.
101111
*/

‎source/mesh_deformation/landlab.cc‎

Lines changed: 42 additions & 31 deletions
Original file line numberDiff line numberDiff line change
@@ -82,42 +82,48 @@ namespace aspect
8282
ExcMessage("The surface diffusion mesh deformation plugin only works for Box geometries."));
8383

8484
unsigned int rank = Utilities::MPI::this_mpi_process(this->get_mpi_communicator());
85-
if (rank >= n_landlab_ranks)
86-
return;
8785

86+
this_rank_runs_landlab = (rank < n_landlab_ranks);
87+
const int color = this_rank_runs_landlab?1:0;
88+
int ierr = MPI_Comm_split(this->get_mpi_communicator(), color, rank, &landlab_communicator);
89+
AssertThrow(ierr == MPI_SUCCESS, ExcMessage("Failed to split MPI communicator for Landlab simulation"));
90+
91+
if (this_rank_runs_landlab)
92+
{
8893
#ifdef ASPECT_WITH_PYTHON
89-
// Append script dirs so env packages (venv site-packages, PYTHONPATH) are found first
90-
// for "import landlab":
91-
PyRun_SimpleString("import sys");
92-
PyRun_SimpleString("sys.path.append(\"" ASPECT_SOURCE_DIR "/tests\")");
93-
PyRun_SimpleString("sys.path.append(\".\")");
94+
// Append script dirs so env packages (venv site-packages, PYTHONPATH) are found first
95+
// for "import landlab":
96+
PyRun_SimpleString("import sys");
97+
PyRun_SimpleString("sys.path.append(\"" ASPECT_SOURCE_DIR "/tests\")");
98+
PyRun_SimpleString("sys.path.append(\".\")");
9499

95-
// avoid floating point exceptions in Landlab Python code:
100+
// avoid floating point exceptions in Landlab Python code:
96101
#ifdef ASPECT_USE_FP_EXCEPTIONS
97-
fedisableexcept(FE_DIVBYZERO|FE_INVALID);
102+
fedisableexcept(FE_DIVBYZERO|FE_INVALID);
98103
#endif
99104

100-
std::cout << "importing '" << script_module_name << "' ..." << std::endl;
101-
pModule = PyImport_ImportModule(script_module_name.c_str());
102-
if (PyErr_Occurred())
103-
PyErr_Print();
104-
AssertThrow(pModule, ExcMessage("Failed to load Python module"));
105+
std::cout << "importing '" << script_module_name << "' ..." << std::endl;
106+
pModule = PyImport_ImportModule(script_module_name.c_str());
107+
if (PyErr_Occurred())
108+
PyErr_Print();
109+
AssertThrow(pModule, ExcMessage("Failed to load Python module"));
105110

106-
// Call Python initialize() function with communicator handle
107-
PyObject *pArgs;
108-
if (n_landlab_ranks == 1)
109-
pArgs = PyTuple_Pack(1, Py_None);
110-
else
111-
pArgs = PyTuple_Pack(1, PyLong_FromLong(MPI_Comm_c2f(this->get_mpi_communicator())));
112-
PyObject *pValue = call_python_function(pModule, "initialize", pArgs);
111+
// Call Python initialize() function with communicator handle
112+
PyObject *pArgs;
113+
if (n_landlab_ranks == 1)
114+
pArgs = PyTuple_Pack(1, Py_None);
115+
else
116+
pArgs = PyTuple_Pack(1, PyLong_FromLong(MPI_Comm_c2f(landlab_communicator)));
117+
PyObject *pValue = call_python_function(pModule, "initialize", pArgs);
113118

114-
Py_DECREF(pArgs);
115-
Py_DECREF(pValue);
119+
Py_DECREF(pArgs);
120+
Py_DECREF(pValue);
116121

117122
#else
118-
AssertThrow(false, ExcMessage("ASPECT needs to be configure with Python support "
119-
"(ASPECT_WITH_PYTHON=ON in CMake) to be able to use the Landlab mesh deformation model."));
123+
AssertThrow(false, ExcMessage("ASPECT needs to be configure with Python support "
124+
"(ASPECT_WITH_PYTHON=ON in CMake) to be able to use the Landlab mesh deformation model."));
120125
#endif
126+
}
121127
}
122128

123129

@@ -130,8 +136,7 @@ namespace aspect
130136
#ifdef ASPECT_WITH_PYTHON
131137
if (!this->remote_point_evaluator)
132138
{
133-
unsigned int rank = Utilities::MPI::this_mpi_process(this->get_mpi_communicator());
134-
if (rank >= n_landlab_ranks)
139+
if (!this_rank_runs_landlab)
135140
{
136141
// This rank does not participate, so we don't own any evaluation points:
137142
std::vector<Point<dim>> surface_points;
@@ -200,7 +205,7 @@ namespace aspect
200205
Assert(current_solution_at_points.size() == this->evaluation_points.size(), ExcInternalError());
201206
std::vector<Tensor<1,dim>> velocities(current_solution_at_points.size(), Tensor<1,dim>());
202207

203-
if (pModule)
208+
if (this_rank_runs_landlab)
204209
{
205210
// Build a dictionary with solution values for each variable to pass to Python:
206211
PyObject *pDict = PyDict_New();
@@ -282,7 +287,7 @@ namespace aspect
282287

283288
// 1. Grab initial deformation from Landlab:
284289
std::vector<Tensor<1,dim>> initial_deformation(this->evaluation_points.size(), Tensor<1,dim>());
285-
if (pModule)
290+
if (this_rank_runs_landlab)
286291
{
287292
Tensor<1,dim> topography_direction;
288293
topography_direction[dim-1] = 1.0;
@@ -344,7 +349,8 @@ namespace aspect
344349
{
345350
prm.declare_entry("MPI ranks for Landlab", "1",
346351
Patterns::Integer(1),
347-
"Number of ranks to use for the Landlab simulation. If set to 1, the Landlab simulation will run sequentially without MPI.");
352+
"Number of ranks to use for the Landlab simulation. If set to 1, the Landlab simulation will run sequentially "
353+
"without MPI. If set to -1, the Landlab simulation will run on all ranks.");
348354
prm.declare_entry("Script path", "",
349355
Patterns::Anything(),
350356
"Path to the Python script to execute. Relative paths and the placeholders "
@@ -371,7 +377,12 @@ namespace aspect
371377
{
372378
prm.enter_subsection ("Landlab");
373379
{
374-
n_landlab_ranks = prm.get_integer("MPI ranks for Landlab");
380+
const int ranks = prm.get_integer("MPI ranks for Landlab");
381+
if (ranks == -1)
382+
n_landlab_ranks = Utilities::MPI::n_mpi_processes(this->get_mpi_communicator());
383+
else
384+
n_landlab_ranks = ranks;
385+
375386
script_path = prm.get("Script path");
376387
script_module_name = prm.get("Script name");
377388
script_argument = prm.get("Script argument");

‎tests/mesh_deformation_external_landlab_02.prm‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
# Test the "landlab" external mesh deformation making use of
2+
# the MPI settings. Note that we are not actually running
3+
# Landlab in parallel here.
4+
#
5+
# Enable if: ASPECT_WITH_PYTHON
6+
# MPI: 3
7+
8+
set Dimension = 3
9+
set Use years instead of seconds = false
10+
set End time = 0.001
11+
set Maximum time step = 0.0005
12+
set Nonlinear solver scheme = no Advection, no Stokes
13+
set Pressure normalization = surface
14+
set Surface pressure = 0
15+
16+
subsection Geometry model
17+
set Model name = box
18+
19+
subsection Box
20+
set X extent = 1
21+
set Y extent = 1
22+
set Z extent = 1
23+
end
24+
end
25+
26+
# Temperature effects are ignored
27+
subsection Initial temperature model
28+
set Model name = function
29+
30+
subsection Function
31+
set Function expression = 100*x+y
32+
end
33+
end
34+
35+
subsection Boundary temperature model
36+
set Fixed temperature boundary indicators = bottom, top
37+
set List of model names = initial temperature
38+
end
39+
40+
# Free slip on all boundaries
41+
subsection Boundary velocity model
42+
set Tangential velocity boundary indicators = left, right, bottom, top
43+
end
44+
45+
subsection Mesh deformation
46+
set Mesh deformation boundary indicators = top: Landlab
47+
set Additional tangential mesh velocity boundary indicators = left, right
48+
subsection Landlab
49+
set Script name = mesh_deformation_external_landlab_02
50+
set MPI ranks for Landlab = 2
51+
end
52+
end
53+
54+
# Vertical gravity
55+
subsection Gravity model
56+
set Model name = vertical
57+
58+
subsection Vertical
59+
set Magnitude = 1
60+
end
61+
end
62+
63+
# One material with unity properties
64+
subsection Material model
65+
set Model name = simple
66+
67+
subsection Simple model
68+
set Reference density = 1
69+
set Reference specific heat = 1
70+
set Reference temperature = 0
71+
set Thermal conductivity = 1
72+
set Thermal expansion coefficient = 1
73+
set Viscosity = 1
74+
end
75+
end
76+
77+
# We also have to specify that we want to use the Boussinesq
78+
# approximation (assuming the density in the temperature
79+
# equation to be constant, and incompressibility).
80+
subsection Formulation
81+
set Formulation = Boussinesq approximation
82+
end
83+
84+
# We here use a globally refined mesh without
85+
# adaptive mesh refinement.
86+
subsection Mesh refinement
87+
set Initial global refinement = 2
88+
set Initial adaptive refinement = 0
89+
set Time steps between mesh refinement = 0
90+
end
91+
92+
subsection Postprocess
93+
set List of postprocessors = visualization, topography
94+
95+
subsection Visualization
96+
set Time between graphical output = 0
97+
set Output mesh velocity = true
98+
set Output mesh displacement = true
99+
set Output undeformed mesh = false
100+
set Interpolate output = false
101+
end
102+
end

‎tests/mesh_deformation_external_landlab_02.py‎

Lines changed: 96 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,96 @@
1+
print("mesh_deformation_external_landlab_02.py")
2+
3+
from mpi4py import MPI
4+
import numpy as np
5+
6+
current_time = 0
7+
8+
comm = None
9+
10+
11+
px = None
12+
py = None
13+
elevation = None
14+
15+
def initialize(comm_handle):
16+
global px, py, elevation
17+
18+
if not comm_handle is None:
19+
# Convert the handle back to an MPI communicator
20+
global comm
21+
comm = MPI.Comm.f2py(comm_handle)
22+
23+
rank = comm.Get_rank()
24+
size = comm.Get_size()
25+
26+
27+
data = 1
28+
globalsum = comm.allreduce(data, op=MPI.SUM)
29+
if comm.rank == 0:
30+
print(f"Python: comm size = {size}, communication test = {globalsum}")
31+
32+
if comm.rank == 0:
33+
px = np.array([0.1, 0.3, 0.4])
34+
py = np.array([0.2, 0.4, 0.5])
35+
36+
if comm.rank == 1:
37+
px = np.array([0.5, 0.7])
38+
py = np.array([0.6, 0.8])
39+
40+
elevation = np.zeros(px.size)
41+
42+
def finalize():
43+
pass
44+
45+
46+
# Run the Landlab simulation from the current time to end_time and return
47+
# the new topographic elevation (in m) at each local node.
48+
# dict_variable_name_to_value_in_nodes is a dictionary mapping variables
49+
# (x velocity, y velocity, temperature, etc.) to an array of values in each
50+
# node.
51+
def update_until(end_time, dict_variable_name_to_value_in_nodes):
52+
53+
print(f"update_until: end_time = {end_time}")
54+
55+
return elevation
56+
57+
def set_mesh_information(dict_grid_information):
58+
pass
59+
60+
# Return the x coordinates of the locally owned nodes on this
61+
# MPI rank. grid_id is always 0.
62+
def get_grid_x(grid_id):
63+
global px
64+
return px
65+
66+
# Return the y coordinates of the locally owned nodes on this
67+
# MPI rank. grid_id is always 0.
68+
def get_grid_y(grid_id):
69+
global py
70+
return py
71+
72+
# Return the initial topography at the start of the simulation
73+
# in each node.
74+
def get_initial_topography(grid_id):
75+
global elevation
76+
return elevation
77+
78+
79+
def write_output():
80+
pass
81+
82+
if __name__ == "__main__":
83+
comm = MPI.COMM_WORLD
84+
initialize(MPI.Comm.py2f(comm))
85+
86+
set_mesh_information({})
87+
print("grid coordinates:", get_grid_x(0), get_grid_y(0))
88+
89+
dt = 0.1
90+
for n in range(3):
91+
data = {}
92+
data["x velocity"] = np.zeros(px.size)
93+
data["y velocity"] = np.zeros(px.size)
94+
data["z velocity"] = np.zeros(px.size)
95+
update_until(n*dt, data)
96+
write_output()

‎tests/mesh_deformation_external_landlab_02/screen-output‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
-----------------------------------------------------------------------------
2+
-----------------------------------------------------------------------------
3+
-----------------------------------------------------------------------------
4+
5+
importing 'mesh_deformation_external_landlab_02' ...
6+
importing 'mesh_deformation_external_landlab_02' ...
7+
mesh_deformation_external_landlab_02.py
8+
mesh_deformation_external_landlab_02.py
9+
Python: comm size = 2, communication test = 2
10+
-----------------------------------------------------------------------------
11+
-----------------------------------------------------------------------------
12+
Number of active cells: 64 (on 3 levels)
13+
Number of degrees of freedom: 3,041 (2,187+125+729)
14+
15+
Number of mesh deformation degrees of freedom: 375
16+
Solving mesh displacement system... 0 iterations.
17+
*** Timestep 0: t=0 seconds, dt=0 seconds
18+
update_until: end_time = 0.0
19+
update_until: end_time = 0.0
20+
Solving mesh displacement system... 0 iterations.
21+
22+
Postprocessing:
23+
Writing graphical output: output-mesh_deformation_external_landlab_02/solution/solution-00000
24+
Topography min/max: 0 m, 0 m
25+
26+
*** Timestep 1: t=0.0005 seconds, dt=0.0005 seconds
27+
update_until: end_time = 0.0005
28+
update_until: end_time = 0.0005
29+
Solving mesh displacement system... 0 iterations.
30+
31+
Postprocessing:
32+
Writing graphical output: output-mesh_deformation_external_landlab_02/solution/solution-00001
33+
Topography min/max: 0 m, 0 m
34+
35+
*** Timestep 2: t=0.001 seconds, dt=0.0005 seconds
36+
update_until: end_time = 0.001
37+
update_until: end_time = 0.001
38+
Solving mesh displacement system... 0 iterations.
39+
40+
Postprocessing:
41+
Writing graphical output: output-mesh_deformation_external_landlab_02/solution/solution-00002
42+
Topography min/max: 0 m, 0 m
43+
44+
Termination requested by criterion: end time
45+
46+
47+
+---------------------------------------------+------------+------------+
48+
+---------------------------------+-----------+------------+------------+
49+
+---------------------------------+-----------+------------+------------+
50+
51+
-----------------------------------------------------------------------------
52+
-----------------------------------------------------------------------------

0 commit comments

Comments
 (0)