diff --git a/pySDC/implementations/convergence_controller_classes/adaptive_collocation.py b/pySDC/implementations/convergence_controller_classes/adaptive_collocation.py index 743335bc51..e914273ee0 100644 --- a/pySDC/implementations/convergence_controller_classes/adaptive_collocation.py +++ b/pySDC/implementations/convergence_controller_classes/adaptive_collocation.py @@ -78,8 +78,39 @@ def setup(self, controller, params, description, **kwargs): defaults['num_colls'] = max([defaults['num_colls'], len(params[key])]) + self.comm = description['sweeper_params'].get('comm', None) + if self.comm: + from mpi4py import MPI + + self.MPI_SUM = MPI.SUM + return {**defaults, **super().setup(controller, params, description, **kwargs)} + def matmul(self, A, b): + """ + Matrix vector multiplication, possibly MPI parallel. + The parallel implementation performs a reduce operation in every row of the matrix. While communicating the + entire vector once could reduce the number of communications, this way we never need to store the entire vector + on any specific rank. + + Args: + A (2d np.ndarray): Matrix + b (list): Vector + + Returns: + List: Axb + """ + if self.comm: + res = [A[i, 0] * b[0] if b[i] is not None else None for i in range(A.shape[0])] + buf = b[0] * 0.0 + for i in range(1, A.shape[1]): + self.comm.Reduce(A[i, self.comm.rank + 1] * b[self.comm.rank + 1], buf, op=self.MPI_SUM, root=i - 1) + if i == self.comm.rank + 1: + res[i] += buf + return res + else: + return A @ b + def switch_sweeper(self, S): """ Update to the next sweeper in line. @@ -105,7 +136,7 @@ def switch_sweeper(self, S): P = L.prob # store solution of current level which will be interpolated to new level - u_old = [me.flatten() for me in L.u] + u_old = [me.flatten() if me is not None else me for me in L.u] nodes_old = L.sweep.coll.nodes.copy() # change sweeper @@ -120,17 +151,19 @@ def switch_sweeper(self, S): nodes_new = L.sweep.coll.nodes.copy() interpolator = LagrangeApproximation(points=np.append(0, nodes_old)) - u_inter = interpolator.getInterpolationMatrix(np.append(0, nodes_new)) @ u_old + u_inter = self.matmul(interpolator.getInterpolationMatrix(np.append(0, nodes_new)), u_old) # assign the interpolated values to the nodes in the level for i in range(0, len(u_inter)): - me = P.dtype_u(P.init) - me[:] = np.reshape(u_inter[i], P.init[0]) - L.u[i] = me + if u_inter[i] is not None: + me = P.dtype_u(P.init) + me[:] = np.reshape(u_inter[i], P.init[0]) + L.u[i] = me # reevaluate rhs for i in range(L.sweep.coll.num_nodes + 1): - L.f[i] = L.prob.eval_f(L.u[i], L.time) + if L.u[i] is not None: + L.f[i] = L.prob.eval_f(L.u[i], L.time) # log the new parameters self.log(f'Switching to collocation {self.status.active_coll + 1} of {self.params.num_colls}', S, level=20) diff --git a/pySDC/implementations/problem_classes/polynomial_test_problem.py b/pySDC/implementations/problem_classes/polynomial_test_problem.py new file mode 100644 index 0000000000..3878227b24 --- /dev/null +++ b/pySDC/implementations/problem_classes/polynomial_test_problem.py @@ -0,0 +1,90 @@ +import numpy as np + +from pySDC.core.Problem import ptype +from pySDC.implementations.datatype_classes.mesh import mesh + + +class polynomial_testequation(ptype): + """ + Dummy problem for tests only! In particular, the `solve_system` function just returns the exact solution instead of + solving an appropriate system. This class is indented to be used for tests of operations that are exact on polynomials. + """ + + dtype_u = mesh + dtype_f = mesh + + def __init__(self, degree=1, seed=26266): + """Initialization routine""" + + # invoke super init, passing number of dofs, dtype_u and dtype_f + super().__init__(init=(1, None, np.dtype('float64'))) + + self.rng = np.random.RandomState(seed=seed) + self.poly = np.polynomial.Polynomial(self.rng.rand(degree)) + self._makeAttributeAndRegister('degree', 'seed', localVars=locals(), readOnly=True) + + def eval_f(self, u, t): + """ + Derivative of the polynomial. + + Parameters + ---------- + u : dtype_u + Current values of the numerical solution. + t : float + Current time of the numerical solution is computed. + + Returns + ------- + f : dtype_f + The right-hand side of the problem. + """ + + f = self.dtype_f(self.init) + f[:] = self.poly.deriv(m=1)(t) + return f + + def solve_system(self, rhs, factor, u0, t): + """ + Just return the exact solution... + + Parameters + ---------- + rhs : dtype_f + Right-hand side for the linear system. + factor : float + Abbrev. for the local stepsize (or any other factor required). + u0 : dtype_u + Initial guess for the iterative solver. + t : float + Current time (e.g. for time-dependent BCs). + + Returns + ------- + me : dtype_u + The solution as mesh. + """ + + return self.u_exact(t) + + def u_exact(self, t, **kwargs): + """ + Evaluate the polynomial. + + Parameters + ---------- + t : float + Time of the exact solution. + u_init : pySDC.problem.testequation0d.dtype_u + Initial solution. + t_init : float + The initial time. + + Returns + ------- + me : dtype_u + The exact solution. + """ + me = self.dtype_u(self.init) + me[:] = self.poly(t) + return me diff --git a/pySDC/implementations/sweeper_classes/generic_implicit.py b/pySDC/implementations/sweeper_classes/generic_implicit.py index 7d15896104..e93060aba5 100644 --- a/pySDC/implementations/sweeper_classes/generic_implicit.py +++ b/pySDC/implementations/sweeper_classes/generic_implicit.py @@ -21,7 +21,7 @@ def __init__(self, params): params['QI'] = 'IE' # call parent's initialization routine - super(generic_implicit, self).__init__(params) + super().__init__(params) # get QI matrix self.QI = self.get_Qdelta_implicit(self.coll, qd_type=self.params.QI) @@ -34,7 +34,6 @@ def integrate(self): list of dtype_u: containing the integral as values """ - # get current level and problem description L = self.level P = L.prob @@ -57,7 +56,6 @@ def update_nodes(self): None """ - # get current level and problem description L = self.level P = L.prob @@ -112,7 +110,6 @@ def compute_end_point(self): None """ - # get current level and problem description L = self.level P = L.prob diff --git a/pySDC/implementations/sweeper_classes/generic_implicit_MPI.py b/pySDC/implementations/sweeper_classes/generic_implicit_MPI.py index 4a661cc9e2..d3aeb38dd8 100644 --- a/pySDC/implementations/sweeper_classes/generic_implicit_MPI.py +++ b/pySDC/implementations/sweeper_classes/generic_implicit_MPI.py @@ -47,7 +47,6 @@ def compute_end_point(self): None """ - # get current level and problem description L = self.level P = L.prob L.uend = P.dtype_u(P.init, val=0.0) @@ -69,7 +68,6 @@ def compute_residual(self, stage=None): stage (str): The current stage of the step the level belongs to """ - # get current level and problem description L = self.level # Check if we want to skip the residual computation to gain performance @@ -108,7 +106,6 @@ def predict(self): and evaluates the RHS of the ODE there """ - # get current level and problem description L = self.level P = L.prob @@ -144,7 +141,6 @@ def integrate(self): list of dtype_u: containing the integral as values """ - # get current level and problem description L = self.level P = L.prob @@ -165,7 +161,6 @@ def update_nodes(self): None """ - # get current level and problem description L = self.level P = L.prob @@ -204,3 +199,30 @@ def update_nodes(self): L.status.updated = True return None + + def compute_end_point(self): + """ + Compute u at the right point of the interval + + The value uend computed here is a full evaluation of the Picard formulation unless do_full_update==False + + Returns: + None + """ + + L = self.level + P = L.prob + L.uend = P.dtype_u(P.init, val=0.0) + + # check if Mth node is equal to right point and do_coll_update is false, perform a simple copy + if self.coll.right_is_node and not self.params.do_coll_update: + super().compute_end_point() + else: + L.uend = P.dtype_u(L.u[0]) + self.params.comm.Allreduce(L.dt * self.coll.weights[self.rank] * L.f[self.rank + 1], L.uend, op=MPI.SUM) + L.uend += L.u[0] + + # add up tau correction of the full interval (last entry) + if L.tau[-1] is not None: + L.uend += L.tau[-1] + return None diff --git a/pySDC/tests/test_convergence_controllers/test_adaptive_collocation.py b/pySDC/tests/test_convergence_controllers/test_adaptive_collocation.py new file mode 100644 index 0000000000..41ab3a83ff --- /dev/null +++ b/pySDC/tests/test_convergence_controllers/test_adaptive_collocation.py @@ -0,0 +1,163 @@ +import pytest + + +def single_run(num_nodes, quad_type, QI, useMPI, params): + """ + Runs a single advection problem with certain parameters + + Args: + num_nodes (int): Number of nodes + quad_type (str): Type of quadrature + QI (str): Preconditioner + useMPI (bool): Whether or not to use MPI + params (dict): Parameters for adaptive collocation convergence controller + + Returns: + (dict): Stats object generated during the run + (pySDC.Controller.controller): Controller used in the run + """ + from pySDC.implementations.problem_classes.polynomial_test_problem import polynomial_testequation + from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + from pySDC.implementations.convergence_controller_classes.adaptive_collocation import AdaptiveCollocation + + if useMPI: + from pySDC.implementations.sweeper_classes.generic_implicit_MPI import generic_implicit_MPI as sweeper_class + from mpi4py import MPI + + comm = MPI.COMM_WORLD + else: + from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit as sweeper_class + + comm = None + + # initialize level parameters + level_params = {} + level_params['dt'] = 1.0 + level_params['restol'] = 1.0 + + # initialize sweeper parameters + sweeper_params = {} + sweeper_params['quad_type'] = quad_type + sweeper_params['num_nodes'] = num_nodes + sweeper_params['QI'] = QI + sweeper_params['comm'] = comm + + problem_params = {'degree': num_nodes} + + # initialize step parameters + step_params = {} + step_params['maxiter'] = 99 + + # initialize controller parameters + controller_params = {} + controller_params['logger_level'] = 30 + + # fill description dictionary for easy step instantiation + description = {} + description['problem_class'] = polynomial_testequation + description['problem_params'] = problem_params + description['sweeper_class'] = sweeper_class + description['sweeper_params'] = sweeper_params + description['level_params'] = level_params + description['step_params'] = step_params + description['convergence_controllers'] = {AdaptiveCollocation: params} + + controller = controller_nonMPI(num_procs=1, controller_params=controller_params, description=description) + return controller + + +def single_test(**kwargs): + """ + Run a single test where the solution is replaced by a polynomial and the nodes are changed. + Because we know the polynomial going in, we can check if the interpolation based change was + exact. If the solution is not a polynomial or a polynomial of higher degree then the number + of nodes, the change in nodes does add some error, of course, but here it is on the order of + machine precision. + """ + import numpy as np + + coll_params_type = { + 'quad_type': ['GAUSS', 'RADAU-RIGHT'], + } + + args = { + 'num_nodes': 3, + 'quad_type': 'RADAU-RIGHT', + 'QI': 'MIN', + 'useMPI': False, + 'params': coll_params_type, + **kwargs, + } + + # prepare variables + controller = single_run(**args) + step = controller.MS[0] + level = step.levels[0] + prob = level.prob + cont = controller.convergence_controllers[ + np.arange(len(controller.convergence_controllers))[ + [type(me).__name__ == 'AdaptiveCollocation' for me in controller.convergence_controllers] + ][0] + ] + nodes = np.append([0], level.sweep.coll.nodes) + + # initialize variables + cont.status.active_coll = 0 + step.status.slot = 0 + level.u[0] = prob.u_exact(t=0) + level.status.time = 0.0 + level.sweep.predict() + for i in range(len(level.u)): + if level.u[i] is not None: + level.u[i][:] = prob.u_exact(nodes[i]) + + # perform the interpolation + cont.switch_sweeper(controller.MS[0]) + cont.status.active_coll = 1 + cont.switch_sweeper(controller.MS[0]) + nodes = np.append([0], level.sweep.coll.nodes) + error = max([abs(level.u[i] - prob.u_exact(nodes[i])) for i in range(len(level.u)) if level.u[i] is not None]) + assert error < 1e-15, f'Interpolation not exact!, Got {error}' + print(f'Passed test with error {error}') + + diff = min([abs(level.u[0] - prob.u_exact(nodes[i])) for i in range(1, len(level.u)) if level.u[i] is not None]) + assert diff > 1e-15, 'Solution is constant!' + + +@pytest.mark.base +def test_adaptive_collocation(): + single_test() + + +@pytest.mark.mpi4py +def test_adaptive_collocation_MPI(): + import subprocess + import os + + num_nodes = 3 + + # Set python path once + my_env = os.environ.copy() + my_env['PYTHONPATH'] = '../../..:.' + my_env['COVERAGE_PROCESS_START'] = 'pyproject.toml' + + cmd = f"mpirun -np {num_nodes} python {__file__} MPI".split() + + p = subprocess.Popen(cmd, env=my_env, cwd=".") + + p.wait() + assert p.returncode == 0, 'ERROR: did not get return code 0, got %s with %2i processes' % ( + p.returncode, + num_nodes, + ) + + +if __name__ == "__main__": + import sys + + kwargs = {} + if len(sys.argv) > 1: + kwargs = { + 'useMPI': True, + } + single_test(**kwargs) diff --git a/pySDC/tests/test_Multistep_sweeper.py b/pySDC/tests/test_sweepers/test_Multistep_sweeper.py similarity index 100% rename from pySDC/tests/test_Multistep_sweeper.py rename to pySDC/tests/test_sweepers/test_Multistep_sweeper.py diff --git a/pySDC/tests/test_Runge_Kutta_sweeper.py b/pySDC/tests/test_sweepers/test_Runge_Kutta_sweeper.py similarity index 100% rename from pySDC/tests/test_Runge_Kutta_sweeper.py rename to pySDC/tests/test_sweepers/test_Runge_Kutta_sweeper.py diff --git a/pySDC/tests/test_sweepers/test_compute_end_point.py b/pySDC/tests/test_sweepers/test_compute_end_point.py new file mode 100644 index 0000000000..e78ccc0085 --- /dev/null +++ b/pySDC/tests/test_sweepers/test_compute_end_point.py @@ -0,0 +1,93 @@ +import pytest + + +def compute_end_point_generic_implicit(num_nodes, useMPI): + """ + Check if the end point is computed exactly for a polynomial of sufficiently low degree. + + Args: + num_nodes (int): Number of collocation nodes + useMPI (bool): Use MPI or not + """ + import numpy as np + from pySDC.implementations.problem_classes.polynomial_test_problem import polynomial_testequation + from pySDC.implementations.controller_classes.controller_nonMPI import controller_nonMPI + + if useMPI: + from pySDC.implementations.sweeper_classes.generic_implicit_MPI import generic_implicit_MPI as sweeper_class + from mpi4py import MPI + + comm = MPI.COMM_WORLD + else: + from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit as sweeper_class + + comm = None + + description = {} + description['sweeper_class'] = sweeper_class + description['sweeper_params'] = {'num_nodes': num_nodes, 'quad_type': 'GAUSS', 'do_coll_update': True, 'comm': comm} + description['problem_params'] = {'degree': num_nodes + 1} + description['level_params'] = {'dt': 1.0} + description['step_params'] = {} + description['problem_class'] = polynomial_testequation + + # prepare variables + controller = controller_nonMPI(1, {'logger_level': 30}, description) + step = controller.MS[0] + level = step.levels[0] + prob = level.prob + sweep = level.sweep + nodes = np.append([0], level.sweep.coll.nodes) + + # initialize variables + for i in range(len(level.u)): + level.u[i] = prob.u_exact(nodes[i]) + level.f[i] = prob.eval_f(level.u[i], nodes[i]) + + # check if the end point was computed exactly + sweep.compute_end_point() + error = abs(level.uend - prob.u_exact(1.0)) + assert error < 1e-15, f'Failed to compute end point! Error: {error}' + print(f'Passed with error to exact end point: {error}') + + assert abs(level.uend - level.u[0]) > 1e-15, 'Solution is constant!' + + +@pytest.mark.base +@pytest.mark.parametrize("num_nodes", [2, 3]) +def test_compute_end_point_generic_implicit(num_nodes): + compute_end_point_generic_implicit(num_nodes, False) + + +@pytest.mark.mpi4py +@pytest.mark.parametrize("num_nodes", [2, 3]) +def test_compute_end_point_generic_implicit_MPI(num_nodes): + import subprocess + import os + + # Set python path once + my_env = os.environ.copy() + my_env['PYTHONPATH'] = '../../..:.' + my_env['COVERAGE_PROCESS_START'] = 'pyproject.toml' + + cmd = f"mpirun -np {num_nodes} python {__file__} MPI {num_nodes}".split() + + p = subprocess.Popen(cmd, env=my_env, cwd=".") + + p.wait() + assert p.returncode == 0, 'ERROR: did not get return code 0, got %s with %2i processes' % ( + p.returncode, + num_nodes, + ) + + +if __name__ == "__main__": + import sys + + kwargs = {} + if len(sys.argv) > 1: + kwargs = { + 'useMPI': bool(sys.argv[1]), + 'num_nodes': int(sys.argv[2]), + } + compute_end_point_generic_implicit(**kwargs) diff --git a/pySDC/tests/test_imexsweeper.py b/pySDC/tests/test_sweepers/test_imexsweeper.py similarity index 100% rename from pySDC/tests/test_imexsweeper.py rename to pySDC/tests/test_sweepers/test_imexsweeper.py diff --git a/pySDC/tests/test_sweeper.py b/pySDC/tests/test_sweepers/test_preconditioners.py similarity index 100% rename from pySDC/tests/test_sweeper.py rename to pySDC/tests/test_sweepers/test_preconditioners.py