Skip to content

Commit 45f5ca0

Browse files
lisawimlisawim
authored
SDC-DAE MPI Sweepers (#439)
* Added FI-SDC-DAE-MPI sweeper * Changes in docu + SI-SDC-DAE-MPI sweeper with playground * Small changes for playgrounds * Test for MPI sweepers + adding MPI stuff to environment.yml * Removed print line * Added test to check uend's from MPI and non-MPI versions + requested changes * What the jellyfish did I do here.. * Changed parametrization in test * Update README.rst
1 parent e372a43 commit 45f5ca0

File tree

9 files changed

+715
-5
lines changed

9 files changed

+715
-5
lines changed

‎pySDC/playgrounds/DAE/DiscontinuousTestDAE.py‎

Lines changed: 2 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -74,7 +74,7 @@ def simulateDAE():
7474
plotSolution(stats)
7575

7676

77-
def plotSolution(stats):
77+
def plotSolution(stats, file_name='data/solution.png'):
7878
r"""
7979
Here, the solution of the DAE is plotted.
8080
@@ -96,7 +96,7 @@ def plotSolution(stats):
9696
plt.xlabel(r"Time $t$")
9797
plt.ylabel(r"Solution $y(t)$, $z(t)$")
9898

99-
plt.savefig('data/solution.png', dpi=300, bbox_inches='tight')
99+
plt.savefig(file_name, dpi=300, bbox_inches='tight')
100100
plt.close()
101101

102102

‎pySDC/playgrounds/DAE/playground_MPI.py‎

Lines changed: 90 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,90 @@
1+
from mpi4py import MPI
2+
from pathlib import Path
3+
4+
from pySDC.projects.DAE.sweepers.SemiImplicitDAEMPI import SemiImplicitDAEMPI
5+
from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE
6+
from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
7+
8+
from pySDC.playgrounds.DAE.DiscontinuousTestDAE import plotSolution
9+
10+
from pySDC.implementations.hooks.log_solution import LogSolution
11+
12+
13+
def run():
14+
r"""
15+
Routine to do a run using the semi-implicit SDC-DAE sweeper enabling parallelization across the nodes. The number of processes that is used to run this file is the number of collocation nodes used! When you run the script with the command
16+
17+
>>> mpiexec -n 3 python3 playground_MPI.py
18+
19+
then 3 collocation nodes are used for SDC.
20+
"""
21+
22+
# This communicator is needed for the SDC-DAE sweeper
23+
comm = MPI.COMM_WORLD
24+
25+
# initialize level parameters
26+
level_params = {
27+
'restol': 1e-12,
28+
'dt': 0.1,
29+
}
30+
31+
# initialize problem parameters
32+
problem_params = {
33+
'newton_tol': 1e-6,
34+
}
35+
36+
# initialize sweeper parameters
37+
sweeper_params = {
38+
'quad_type': 'RADAU-RIGHT',
39+
'num_nodes': comm.Get_size(),
40+
'QI': 'MIN-SR-S', # use a diagonal Q_Delta here!
41+
'initial_guess': 'spread',
42+
'comm': comm,
43+
}
44+
45+
# check if number of processes requested matches with number of nodes
46+
assert sweeper_params['num_nodes'] == comm.Get_size(), f"Number of nodes does not match with number of processes!"
47+
48+
# initialize step parameters
49+
step_params = {
50+
'maxiter': 20,
51+
}
52+
53+
# initialize controller parameters
54+
controller_params = {
55+
'logger_level': 30,
56+
'hook_class': [LogSolution],
57+
}
58+
59+
# fill description dictionary for easy step instantiation
60+
description = {
61+
'problem_class': DiscontinuousTestDAE,
62+
'problem_params': problem_params,
63+
'sweeper_class': SemiImplicitDAEMPI,
64+
'sweeper_params': sweeper_params,
65+
'level_params': level_params,
66+
'step_params': step_params,
67+
}
68+
69+
# instantiate controller
70+
controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
71+
P = controller.MS[0].levels[0].prob
72+
73+
t0 = 1.0
74+
Tend = 3.0
75+
76+
uinit = P.u_exact(t0)
77+
78+
# call main function to get things done...
79+
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
80+
81+
Path("data").mkdir(parents=True, exist_ok=True)
82+
83+
# only process with index 0 should plot
84+
if comm.Get_rank() == 0:
85+
file_name = 'data/solution_MPI.png'
86+
plotSolution(stats, file_name)
87+
88+
89+
if __name__ == "__main__":
90+
run()

‎pySDC/projects/DAE/README.rst‎

Lines changed: 9 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -35,11 +35,19 @@ Project overview
3535
- | ``fully_implicit_dae_playground.py``
3636
| Testing arena for the fully implicit sweeper.
3737
- | ``synchronous_machine_playground.py``
38-
| Testing arena for the synchronous machine model.
38+
| Testing arena for the synchronous machine model.
39+
- | ``accuracy_check_MPI.py``
40+
| Script checking the order of accuracy of MPI sweepers for DAEs of different indices.
3941
4042
- sweepers
4143
- | ``fully_implicit_DAE.py``
4244
| Sweeper that accepts a fully implicit formulation of a system of differential equations and applies to it a modified version of spectral deferred correction
45+
- | ``SemiImplicitDAE.py``
46+
| SDC sweeper especially for semi-explicit DAEs. This sweeper is based on ideas mentioned in `Huang et al. (2007) <https://www.sciencedirect.com/science/article/abs/pii/S0021999106003147>`_.
47+
- | ``fully_implicit_DAE_MPI.py``
48+
| MPI version of fully-implicit SDC-DAE sweeper.
49+
- | ``SemiImplicitDAEMPI.py``
50+
| MPI version of semi-implicit SDC-DAE sweeper.
4351
4452
Theoretical details
4553
----------------------

‎pySDC/projects/DAE/environment.yml‎

Lines changed: 2 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -8,3 +8,5 @@ dependencies:
88
- numpy
99
- scipy>=0.17.1
1010
- dill>=0.2.6
11+
- mpich
12+
- mpi4py>=3.0.0

‎pySDC/projects/DAE/run/accuracy_check_MPI.py‎

Lines changed: 187 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,187 @@
1+
import numpy as np
2+
from mpi4py import MPI
3+
4+
from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI
5+
from pySDC.projects.DAE.misc.HookClass_DAE import (
6+
LogGlobalErrorPostStepDifferentialVariable,
7+
LogGlobalErrorPostStepAlgebraicVariable,
8+
)
9+
from pySDC.helpers.stats_helper import get_sorted
10+
11+
12+
def run(dt, num_nodes, use_MPI, semi_implicit, residual_type, index_case, initial_guess='spread', comm=None):
13+
r"""
14+
Prepares the controller with all the description needed. Here, the function decides to choose the correct sweeper
15+
for the test.
16+
17+
Parameters
18+
----------
19+
dt : float
20+
Time step size chosen for simulation.
21+
num_nodes : int
22+
Number of collocation nodes.
23+
use_MPI : bool
24+
If True, the MPI sweeper classes are used.
25+
semi_implicit : bool
26+
Modules are loaded either for fully-implicit case or semi-implicit case.
27+
residual_type : str
28+
Choose how to compute the residual.
29+
index_case : int
30+
Denotes the index case of a DAE to be tested here, can be either ``1`` or ``2``.
31+
initial_guess : str, optional
32+
Type of initial guess for simulation.
33+
comm : mpi4py.MPI.COMM_WORLD
34+
Communicator.
35+
"""
36+
37+
if not semi_implicit:
38+
if use_MPI:
39+
from pySDC.projects.DAE.sweepers.fully_implicit_DAE_MPI import fully_implicit_DAE_MPI as sweeper
40+
41+
else:
42+
from pySDC.projects.DAE.sweepers.fully_implicit_DAE import fully_implicit_DAE as sweeper
43+
44+
else:
45+
if use_MPI:
46+
from pySDC.projects.DAE.sweepers.SemiImplicitDAEMPI import SemiImplicitDAEMPI as sweeper
47+
48+
else:
49+
from pySDC.projects.DAE.sweepers.SemiImplicitDAE import SemiImplicitDAE as sweeper
50+
51+
if index_case == 1:
52+
from pySDC.projects.DAE.problems.DiscontinuousTestDAE import DiscontinuousTestDAE as problem
53+
54+
t0 = 1.0
55+
Tend = 1.5
56+
57+
elif index_case == 2:
58+
from pySDC.projects.DAE.problems.simple_DAE import simple_dae_1 as problem
59+
60+
t0 = 0.0
61+
Tend = 0.4
62+
63+
else:
64+
raise NotImplementedError(f"DAE case of index {index_case} is not implemented!")
65+
66+
# initialize level parameters
67+
level_params = {
68+
'restol': 1e-12,
69+
'residual_type': residual_type,
70+
'dt': dt,
71+
}
72+
73+
# initialize problem parameters
74+
problem_params = {
75+
'newton_tol': 1e-6,
76+
}
77+
78+
# initialize sweeper parameters
79+
sweeper_params = {
80+
'quad_type': 'RADAU-RIGHT',
81+
'num_nodes': num_nodes,
82+
'QI': 'MIN-SR-S', # use a diagonal Q_Delta here!
83+
'initial_guess': initial_guess,
84+
}
85+
86+
# check if number of processes requested matches with number of nodes
87+
if comm is not None:
88+
sweeper_params.update({'comm': comm})
89+
assert (
90+
sweeper_params['num_nodes'] == comm.Get_size()
91+
), f"Number of nodes does not match with number of processes! Expected {sweeper_params['num_nodes']}, got {comm.Get_size()}!"
92+
93+
# initialize step parameters
94+
step_params = {
95+
'maxiter': 20,
96+
}
97+
98+
# initialize controller parameters
99+
controller_params = {
100+
'logger_level': 30,
101+
'hook_class': [LogGlobalErrorPostStepDifferentialVariable, LogGlobalErrorPostStepAlgebraicVariable],
102+
}
103+
104+
# fill description dictionary for easy step instantiation
105+
description = {
106+
'problem_class': problem,
107+
'problem_params': problem_params,
108+
'sweeper_class': sweeper,
109+
'sweeper_params': sweeper_params,
110+
'level_params': level_params,
111+
'step_params': step_params,
112+
}
113+
114+
# instantiate controller
115+
controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description)
116+
P = controller.MS[0].levels[0].prob
117+
118+
uinit = P.u_exact(t0)
119+
120+
# call main function to get things done...
121+
uend, stats = controller.run(u0=uinit, t0=t0, Tend=Tend)
122+
controller.MS[0].levels[0].sweep.compute_end_point()
123+
124+
residual = controller.MS[0].levels[0].status.residual
125+
126+
return uend, residual, stats
127+
128+
129+
def check_order(comm):
130+
num_nodes = comm.Get_size()
131+
use_MPI = True
132+
residual_type = 'full_abs'
133+
for semi_implicit in [False, True]:
134+
for index_case in [1, 2]:
135+
dt_list = np.logspace(-1.7, -1.0, num=5)
136+
137+
errorsDiff, errorsAlg = np.zeros(len(dt_list)), np.zeros(len(dt_list))
138+
for i, dt in enumerate(dt_list):
139+
_, _, stats = run(
140+
dt=dt,
141+
num_nodes=num_nodes,
142+
use_MPI=use_MPI,
143+
semi_implicit=semi_implicit,
144+
residual_type=residual_type,
145+
index_case=index_case,
146+
comm=comm,
147+
)
148+
149+
errorsDiff[i] = max(
150+
np.array(
151+
get_sorted(stats, type='e_global_differential_post_step', sortby='time', recomputed=False)
152+
)[:, 1]
153+
)
154+
errorsAlg[i] = max(
155+
np.array(get_sorted(stats, type='e_global_algebraic_post_step', sortby='time', recomputed=False))[
156+
:, 1
157+
]
158+
)
159+
160+
# only process with index 0 should plot
161+
if comm.Get_rank() == 0:
162+
orderDiff = np.mean(
163+
[
164+
np.log(errorsDiff[i] / errorsDiff[i - 1]) / np.log(dt_list[i] / dt_list[i - 1])
165+
for i in range(1, len(dt_list))
166+
]
167+
)
168+
orderAlg = np.mean(
169+
[
170+
np.log(errorsAlg[i] / errorsAlg[i - 1]) / np.log(dt_list[i] / dt_list[i - 1])
171+
for i in range(1, len(dt_list))
172+
]
173+
)
174+
175+
refOrderDiff = 2 * comm.Get_size() - 1
176+
refOrderAlg = 2 * comm.Get_size() - 1 if index_case == 1 else comm.Get_size()
177+
assert np.isclose(
178+
orderDiff, refOrderDiff, atol=1e0
179+
), f"Expected order {refOrderDiff} in differential variable, got {orderDiff}"
180+
assert np.isclose(
181+
orderAlg, refOrderAlg, atol=1e0
182+
), f"Expected order {refOrderAlg} in algebraic variable, got {orderAlg}"
183+
184+
185+
if __name__ == "__main__":
186+
comm = MPI.COMM_WORLD
187+
check_order(comm)

0 commit comments

Comments
 (0)