Coverage for pySDC/projects/DAE/sweepers/fullyImplicitDAEMPI.py: 97%
58 statements
« prev ^ index » next coverage.py v7.6.7, created at 2024-11-16 14:51 +0000
1from mpi4py import MPI
3from pySDC.core.errors import ParameterError
4from pySDC.implementations.sweeper_classes.generic_implicit_MPI import SweeperMPI, generic_implicit_MPI
5from pySDC.projects.DAE.sweepers.fullyImplicitDAE import FullyImplicitDAE
8class SweeperDAEMPI(SweeperMPI):
9 r"""
10 MPI base class for DAE sweepers where each rank administers one collocation node.
12 This class provides all the stuff to be done in parallel during the simulation. When implementing a new MPI-DAE sweeper multiple inheritance is used here
14 >>> class FullyImplicitDAEMPI(SweeperDAEMPI, generic_implicit_MPI, FullyImplicitDAE):
16 Be careful with ordering of classes where the child class should inherit from! Due to multiple inheritance several methods are overwritten here. In the example above (which is the class below) the class first inherits the methods
18 - compute_residual()
19 - predict()
21 from ``SweeperDAEMPI``. ``compute_end_point()`` is inherited from ``SweeperMPI``. From second inherited class ``generic_implicit_MPI`` the methods
23 - integrate()
24 - update_nodes()
26 are inherited, which then overwritten by the child class.
28 Parameters
29 ----------
30 params : dict
31 Parameters passed to the sweeper.
32 """
34 def compute_residual(self, stage=None):
35 r"""
36 Uses the absolute value of the DAE system
38 .. math::
39 ||F(t, u, u')||
41 for computing the residual in a chosen norm. If norm computes residual by using all values on collocation nodes, result is broadcasted to all processes; when norm only uses the value on last node, the result is collected on last process!
43 Parameters
44 ----------
45 stage : str, optional
46 The current stage of the step the level belongs to.
47 """
49 L = self.level
50 P = L.prob
52 # Check if we want to skip the residual computation to gain performance
53 # Keep in mind that skipping any residual computation is likely to give incorrect outputs of the residual!
54 if stage in self.params.skip_residual_computation:
55 L.status.residual = 0.0 if L.status.residual is None else L.status.residual
56 return None
58 res = P.eval_f(L.u[self.rank + 1], L.f[self.rank + 1], L.time + L.dt * self.coll.nodes[self.rank])
59 res_norm = abs(res) # use abs function from data type here
61 # find maximal residual over the nodes
62 if L.params.residual_type == 'full_abs':
63 L.status.residual = self.comm.allreduce(res_norm, op=MPI.MAX)
64 elif L.params.residual_type == 'last_abs':
65 L.status.residual = self.comm.bcast(res_norm, root=self.comm.size - 1)
66 elif L.params.residual_type == 'full_rel':
67 L.status.residual = self.comm.allreduce(res_norm / abs(L.u[0]), op=MPI.MAX)
68 elif L.params.residual_type == 'last_rel':
69 L.status.residual = self.comm.bcast(res_norm / abs(L.u[0]), root=self.comm.size - 1)
70 else:
71 raise NotImplementedError(f'residual type \"{L.params.residual_type}\" not implemented!')
73 # indicate that the residual has seen the new values
74 L.status.updated = False
76 return None
78 def predict(self):
79 r"""
80 Predictor to fill values at nodes before first sweep. It can decide whether the
82 - initial condition is spread to each node ('initial_guess' = 'spread'),
83 - or zero values are spread to each node ('initial_guess' = 'zero').
85 Default prediction for sweepers, only copies values to all collocation nodes. This function
86 overrides the base implementation by always initialising ``level.f`` to zero. This is necessary since
87 ``level.f`` stores the solution derivative in the fully implicit case, which is not initially known.
88 """
90 L = self.level
91 P = L.prob
93 # set initial guess for gradient to zero
94 L.f[0] = P.dtype_f(init=P.init, val=0.0)
96 if self.params.initial_guess == 'spread':
97 L.u[self.rank + 1] = P.dtype_u(L.u[0])
98 L.f[self.rank + 1] = P.dtype_f(init=P.init, val=0.0)
99 elif self.params.initial_guess == 'zero':
100 L.u[self.rank + 1] = P.dtype_u(init=P.init, val=0.0)
101 L.f[self.rank + 1] = P.dtype_f(init=P.init, val=0.0)
102 else:
103 raise ParameterError(f'initial_guess option {self.params.initial_guess} not implemented')
105 # indicate that this level is now ready for sweeps
106 L.status.unlocked = True
107 L.status.updated = True
110class FullyImplicitDAEMPI(SweeperDAEMPI, generic_implicit_MPI, FullyImplicitDAE):
111 r"""
112 Custom sweeper class to implement the fully-implicit SDC parallelized across the nodes for solving fully-implicit DAE problems of the form
114 .. math::
115 F(t, u, u') = 0.
117 More detailed description can be found in ``fullyImplicitDAE.py``.
119 This sweeper uses diagonal matrices :math:`\mathbf{Q}_\Delta` to
120 solve the implicit system on each node in parallel. First diagonal matrices were first developed and applied in [1]_. Years later intensive theory in [2]_ is developed to that topic. Ideas of parallelizing SDC across the method are basically applied for the DAE case here.
122 Note
123 ----
124 Please supply a communicator as `comm` to the parameters! The number of processes used to parallelize needs to be equal to the number of collocation nodes used. For example, three collocation nodes are used and passed to the sweeper via
126 >>> sweeper_params = {'num_nodes': 3}
128 running a script ``example_mpi.py`` using ``mpi4py`` that is used to enable parallelism have to be done by using the command
130 >>> mpiexec -n 3 python3 example.py
132 Reference
133 ---------
134 .. [1] R. Speck. Parallelizing spectral deferred corrections across the method. Comput. Vis. Sci. 19, No. 3-4, 75-83 (2018).
135 .. [2] G. Čaklović, T. Lunet, S. Götschel, D. Ruprecht. Improving Efficiency of Parallel Across the Method Spectral Deferred Corrections. Preprint, arXiv:2403.18641 [math.NA] (2024).
136 """
138 def integrate(self, last_only=False):
139 r"""
140 Integrates the gradient. ``me`` serves as buffer, and root process ``m`` stores
141 the result of integral at node :math:`\tau_m` there.
143 Parameters
144 ----------
145 last_only : bool, optional
146 Integrate only the last node for the residual or all of them.
148 Returns
149 -------
150 me : list of dtype_u
151 Containing the integral as values.
152 """
154 L = self.level
155 P = L.prob
157 me = P.dtype_u(P.init, val=0.0)
158 for m in [self.coll.num_nodes - 1] if last_only else range(self.coll.num_nodes):
159 integral = L.dt * self.coll.Qmat[m + 1, self.rank + 1] * L.f[self.rank + 1][:]
160 recvBuf = me[:] if m == self.rank else None
161 self.comm.Reduce(integral, recvBuf, root=m, op=MPI.SUM)
163 return me
165 def update_nodes(self):
166 r"""
167 Updates values of ``u`` and ``f`` at collocation nodes. This correspond to a single iteration of the
168 preconditioned Richardson iteration in **"ordinary"** SDC.
169 """
171 L = self.level
172 P = L.prob
174 # only if the level has been touched before
175 assert L.status.unlocked
177 integral = self.integrate()
178 integral -= L.dt * self.QI[self.rank + 1, self.rank + 1] * L.f[self.rank + 1]
179 integral += L.u[0]
181 u_approx = P.dtype_u(integral)
183 L.f[self.rank + 1][:] = P.solve_system(
184 FullyImplicitDAE.F,
185 u_approx,
186 L.dt * self.QI[self.rank + 1, self.rank + 1],
187 L.f[self.rank + 1],
188 L.time + L.dt * self.coll.nodes[self.rank],
189 )
191 integral = self.integrate()
192 L.u[self.rank + 1] = L.u[0] + integral
193 # indicate presence of new values at this level
194 L.status.updated = True
196 return None