From 6ae477946e9a0041ae77c69bcabaea4dae3184a2 Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 31 Oct 2023 08:03:27 +0100 Subject: [PATCH 1/6] Added solve_system method to all DAE classes --- pySDC/projects/DAE/problems/simple_DAE.py | 126 +++++++++++++++++- .../DAE/problems/synchronous_machine.py | 40 ++++++ .../DAE/problems/transistor_amplifier.py | 78 +++++++++++ .../DAE/sweepers/fully_implicit_DAE.py | 35 +++-- 4 files changed, 267 insertions(+), 12 deletions(-) diff --git a/pySDC/projects/DAE/problems/simple_DAE.py b/pySDC/projects/DAE/problems/simple_DAE.py index 78716b5343..859dbf7c5f 100644 --- a/pySDC/projects/DAE/problems/simple_DAE.py +++ b/pySDC/projects/DAE/problems/simple_DAE.py @@ -1,8 +1,10 @@ import warnings import numpy as np from scipy.interpolate import interp1d +from scipy.optimize import root from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae +from pySDC.core.Problem import WorkCounter class pendulum_2d(ptype_dae): @@ -32,6 +34,12 @@ class pendulum_2d(ptype_dae): ---------- t_end: float The end time at which the reference solution is determined. + Attributes + ---------- + work_counters : WorkCounter + Counts the work, i.e., number of function calls of right-hand side is called and stored in + ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in + ``work_counters['newton']``. References ---------- @@ -49,6 +57,8 @@ def __init__(self, nvars, newton_tol): # solution = data[:, 1:] # self.u_ref = interp1d(t, solution, kind='cubic', axis=0, fill_value='extrapolate') self.t_end = 0.0 + self.work_counters['newton'] = WorkCounter() + self.work_counters['rhs'] = WorkCounter() def eval_f(self, u, du, t): r""" @@ -73,8 +83,39 @@ def eval_f(self, u, du, t): # weight somehow f = self.dtype_f(self.init) f[:] = (du[0] - u[2], du[1] - u[3], du[2] + u[4] * u[0], du[3] + u[4] * u[1] + g, u[0] ** 2 + u[1] ** 2 - 1) + self.work_counters['rhs']() return f + def solve_system(self, impl_sys, u0, t): + r""" + Solver for nonlinear implicit system (defined in sweeper). + + Parameters + ---------- + impl_sys : callable + Implicit system to be solved. + u0 : dtype_u + Initial guess for solver. + t : float + Current time :math:`t`. + + Returns + ------- + me : dtype_u + Numerical solution. + """ + + me = self.dtype_u(self.init) + opt = root( + impl_sys, + u0, + method='hybr', + tol=self.newton_tol, + ) + me[:] = opt.x + self.work_counters['newton'].niter += opt.nfev + return me + def u_exact(self, t): """ Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical reference solution. @@ -132,12 +173,26 @@ class simple_dae_1(ptype_dae): newton_tol : float Tolerance for Newton solver. + Attributes + ---------- + work_counters : WorkCounter + Counts the work, i.e., number of function calls of right-hand side is called and stored in + ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in + ``work_counters['newton']``. + References ---------- .. [1] U. Ascher, L. R. Petzold. Computer method for ordinary differential equations and differential-algebraic equations. Society for Industrial and Applied Mathematics (1998). """ + def __init__(self, nvars, newton_tol): + """Initialization routine""" + super().__init__(nvars, newton_tol) + + self.work_counters['newton'] = WorkCounter() + self.work_counters['rhs'] = WorkCounter() + def eval_f(self, u, du, t): r""" Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. @@ -164,8 +219,39 @@ def eval_f(self, u, du, t): -du[1] + (1 - a) / (t - 2) * u[0] - u[1] + (a - 1) * u[2] + 2 * np.exp(t), (t + 2) * u[0] + (t**2 - 4) * u[1] - (t**2 + t - 2) * np.exp(t), ) + self.work_counters['rhs']() return f + def solve_system(self, impl_sys, u0, t): + r""" + Solver for nonlinear implicit system (defined in sweeper). + + Parameters + ---------- + impl_sys : callable + Implicit system to be solved. + u0 : dtype_u + Initial guess for solver. + t : float + Current time :math:`t`. + + Returns + ------- + me : dtype_u + Numerical solution. + """ + + me = self.dtype_u(self.init) + opt = root( + impl_sys, + u0, + method='hybr', + tol=self.newton_tol, + ) + me[:] = opt.x + self.work_counters['newton'].niter += opt.nfev + return me + def u_exact(self, t): """ Routine for the exact solution. @@ -207,8 +293,12 @@ class problematic_f(ptype_dae): Attributes ---------- - eta: float + eta : float Specific parameter of the problem. + work_counters : WorkCounter + Counts the work, i.e., number of function calls of right-hand side is called and stored in + ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in + ``work_counters['newton']``. References ---------- @@ -221,6 +311,9 @@ def __init__(self, nvars, newton_tol, eta=1): super().__init__(nvars, newton_tol) self._makeAttributeAndRegister('eta', localVars=locals()) + self.work_counters['rhs'] = WorkCounter() + self.work_counters['newton'] = WorkCounter() + def eval_f(self, u, du, t): r""" Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. @@ -244,8 +337,39 @@ def eval_f(self, u, du, t): u[0] + self.eta * t * u[1] - np.sin(t), du[0] + self.eta * t * du[1] + (1 + self.eta) * u[1] - np.cos(t), ) + self.work_counters['rhs']() return f + def solve_system(self, impl_sys, u0, t): + r""" + Solver for nonlinear implicit system (defined in sweeper). + + Parameters + ---------- + impl_sys : callable + Implicit system to be solved. + u0 : dtype_u + Initial guess for solver. + t : float + Current time :math:`t`. + + Returns + ------- + me : dtype_u + Numerical solution. + """ + + me = self.dtype_u(self.init) + opt = root( + impl_sys, + u0, + method='hybr', + tol=self.newton_tol, + ) + me[:] = opt.x + self.work_counters['newton'].niter += opt.nfev + return me + def u_exact(self, t): """ Routine for the exact solution. diff --git a/pySDC/projects/DAE/problems/synchronous_machine.py b/pySDC/projects/DAE/problems/synchronous_machine.py index 1cd1855455..1f99f5644a 100644 --- a/pySDC/projects/DAE/problems/synchronous_machine.py +++ b/pySDC/projects/DAE/problems/synchronous_machine.py @@ -1,8 +1,10 @@ import numpy as np import warnings from scipy.interpolate import interp1d +from scipy.optimize import root from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae +from pySDC.core.Problem import WorkCounter from pySDC.implementations.datatype_classes.mesh import mesh @@ -143,6 +145,10 @@ class synchronous_machine_infinite_bus(ptype_dae): Voltage at the field winding. T_m: float Defines the mechanical torque applied to the rotor shaft. + work_counters : WorkCounter + Counts the work, i.e., number of function calls of right-hand side is called and stored in + ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in + ``work_counters['newton']``. References ---------- @@ -184,6 +190,9 @@ def __init__(self, nvars, newton_tol): self.v_F = 8.736809687330562e-4 self.T_m = 0.854 + self.work_counters['rhs'] = WorkCounter() + self.work_counters['newton'] = WorkCounter() + def eval_f(self, u, du, t): r""" Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. @@ -257,8 +266,39 @@ def eval_f(self, u, du, t): -psi_Q1 + self.L_mq * i_q + self.L_Q1 * i_Q1 + self.L_mq * i_Q2, -psi_Q2 + self.L_mq * i_q + self.L_mq * i_Q1 + self.L_Q2 * i_Q2, ) + self.work_counters['rhs']() return f + def solve_system(self, impl_sys, u0, t): + r""" + Solver for nonlinear implicit system (defined in sweeper). + + Parameters + ---------- + impl_sys : callable + Implicit system to be solved. + u0 : dtype_u + Initial guess for solver. + t : float + Current time :math:`t`. + + Returns + ------- + me : dtype_u + Numerical solution. + """ + + me = self.dtype_u(self.init) + opt = root( + impl_sys, + u0, + method='hybr', + tol=self.newton_tol, + ) + me[:] = opt.x + self.work_counters['newton'].niter += opt.nfev + return me + def u_exact(self, t): """ Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical reference solution. diff --git a/pySDC/projects/DAE/problems/transistor_amplifier.py b/pySDC/projects/DAE/problems/transistor_amplifier.py index 2ea06b632b..b160be8ef1 100644 --- a/pySDC/projects/DAE/problems/transistor_amplifier.py +++ b/pySDC/projects/DAE/problems/transistor_amplifier.py @@ -1,8 +1,10 @@ import warnings import numpy as np +from scipy.optimize import root from scipy.interpolate import interp1d from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae +from pySDC.core.Problem import WorkCounter # Helper function @@ -54,6 +56,10 @@ class one_transistor_amplifier(ptype_dae): ---------- t_end: float The end time at which the reference solution is determined. + work_counters : WorkCounter + Counts the work, i.e., number of function calls of right-hand side is called and stored in + ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in + ``work_counters['newton']``. References ---------- @@ -72,6 +78,9 @@ def __init__(self, nvars, newton_tol): # self.u_ref = interp1d(x, y, kind='cubic', axis=0, fill_value='extrapolate') self.t_end = 0.0 + self.work_counters['rhs'] = WorkCounter() + self.work_counters['newton'] = WorkCounter() + def eval_f(self, u, du, t): r""" Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. @@ -104,8 +113,39 @@ def eval_f(self, u, du, t): (u_b - u[3]) / r_k + c_3 * (du[4] - du[3]) - alpha * _transistor(u[1] - u[2]), -u[4] / r_k + c_3 * (du[3] - du[4]), ) + self.work_counters['rhs']() return f + def solve_system(self, impl_sys, u0, t): + r""" + Solver for nonlinear implicit system (defined in sweeper). + + Parameters + ---------- + impl_sys : callable + Implicit system to be solved. + u0 : dtype_u + Initial guess for solver. + t : float + Current time :math:`t`. + + Returns + ------- + me : dtype_u + Numerical solution. + """ + + me = self.dtype_u(self.init) + opt = root( + impl_sys, + u0, + method='hybr', + tol=self.newton_tol, + ) + me[:] = opt.x + self.work_counters['newton'].niter += opt.nfev + return me + def u_exact(self, t): """ Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical @@ -187,6 +227,10 @@ class two_transistor_amplifier(ptype_dae): ---------- t_end: float The end time at which the reference solution is determined. + work_counters : WorkCounter + Counts the work, i.e., number of function calls of right-hand side is called and stored in + ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in + ``work_counters['newton']``. References ---------- @@ -205,6 +249,9 @@ def __init__(self, nvars, newton_tol): # self.u_ref = interp1d(x, y, kind='cubic', axis=0, fill_value='extrapolate') self.t_end = 0.0 + self.work_counters['rhs'] = WorkCounter() + self.work_counters['newton'] = WorkCounter() + def eval_f(self, u, du, t): r""" Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. @@ -240,8 +287,39 @@ def eval_f(self, u, du, t): (u_b - u[6]) / r_k - c_5 * (du[6] - du[7]) - alpha * _transistor(u[4] - u[5]), -u[7] / r_k + c_5 * (du[6] - du[7]), ) + self.work_counters['rhs']() return f + def solve_system(self, impl_sys, u0, t): + r""" + Solver for nonlinear implicit system (defined in sweeper). + + Parameters + ---------- + impl_sys : callable + Implicit system to be solved. + u0 : dtype_u + Initial guess for solver. + t : float + Current time :math:`t`. + + Returns + ------- + me : dtype_u + Numerical solution. + """ + + me = self.dtype_u(self.init) + opt = root( + impl_sys, + u0, + method='hybr', + tol=self.newton_tol, + ) + me[:] = opt.x + self.work_counters['newton'].niter += opt.nfev + return me + def u_exact(self, t): """ Dummy exact solution that should only be used to get initial conditions for the problem. This makes diff --git a/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py b/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py index 28461a5634..57f4b2fce1 100644 --- a/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py +++ b/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py @@ -104,16 +104,33 @@ def update_nodes(self): u_approx += L.dt * self.QI[m, j] * L.f[j] # params contains U = u' - def impl_fn(params): - # make params into a mesh object + def implSystem(params): + """ + Build implicit system to solve in order to find the unknowns. + + Parameters + ---------- + params : dtype_u + Unknowns of the system. + + Returns + ------- + sys : + System to be solved as implicit function. + """ + params_mesh = P.dtype_f(P.init) params_mesh[:] = params + # build parameters to pass to implicit function local_u_approx = u_approx + # note that derivatives of algebraic variables are taken into account here too # these do not directly affect the output of eval_f but rather indirectly via QI local_u_approx += L.dt * self.QI[m, m] * params_mesh - return P.eval_f(local_u_approx, params_mesh, L.time + L.dt * self.coll.nodes[m - 1]) + + sys = P.eval_f(local_u_approx, params_mesh, L.time + L.dt * self.coll.nodes[m - 1]) + return sys # get U_k+1 # note: not using solve_system here because this solve step is the same for any problem @@ -121,15 +138,11 @@ def impl_fn(params): # https://github.com/scipy/scipy/blob/8a6f1a0621542f059a532953661cd43b8167fce0/scipy/optimize/_root.py#L220 # options['xtol'] = P.params.newton_tol # options['eps'] = 1e-16 - opt = optimize.root( - impl_fn, - L.f[m], - method='hybr', - tol=P.newton_tol - # callback= lambda x, f: print("solution:", x, " residual: ", f) - ) + + u_new = P.solve_system(implSystem, L.f[m], L.time + L.dt * self.coll.nodes[m - 1]) + # update gradient (recall L.f is being used to store the gradient) - L.f[m][:] = opt.x + L.f[m][:] = u_new # opt.x # Update solution approximation integral = self.integrate() From c325087ca5a890ffb1c9ced2606fa1776dd2f686 Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 31 Oct 2023 08:14:19 +0100 Subject: [PATCH 2/6] Removed comment --- pySDC/projects/DAE/sweepers/fully_implicit_DAE.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py b/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py index 57f4b2fce1..414fde4545 100644 --- a/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py +++ b/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py @@ -142,7 +142,7 @@ def implSystem(params): u_new = P.solve_system(implSystem, L.f[m], L.time + L.dt * self.coll.nodes[m - 1]) # update gradient (recall L.f is being used to store the gradient) - L.f[m][:] = u_new # opt.x + L.f[m][:] = u_new # Update solution approximation integral = self.integrate() From 76efa6aaf9a6e022d831f86867ed5825085eca31 Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 31 Oct 2023 08:44:52 +0100 Subject: [PATCH 3/6] Moved solve_system to generic ProblemDAE class --- pySDC/projects/DAE/misc/ProblemDAE.py | 61 +++++++++-- pySDC/projects/DAE/problems/simple_DAE.py | 103 +----------------- .../DAE/problems/synchronous_machine.py | 35 +----- .../DAE/problems/transistor_amplifier.py | 69 +----------- 4 files changed, 56 insertions(+), 212 deletions(-) diff --git a/pySDC/projects/DAE/misc/ProblemDAE.py b/pySDC/projects/DAE/misc/ProblemDAE.py index 97c64b9d87..91cc6bc148 100644 --- a/pySDC/projects/DAE/misc/ProblemDAE.py +++ b/pySDC/projects/DAE/misc/ProblemDAE.py @@ -1,25 +1,64 @@ import numpy as np +from scipy.optimize import root -from pySDC.core.Problem import ptype +from pySDC.core.Problem import ptype, WorkCounter from pySDC.implementations.datatype_classes.mesh import mesh class ptype_dae(ptype): - """ - Interface class for DAE problems. Ensures that all parameters are passed that are needed by DAE sweepers + r""" + This class implements a generic DAE class and illustrates the interface class for DAE problems. + It ensures that all parameters are passed that are needed by DAE sweepers. + + Parameters + ---------- + nvars : int + Number of unknowns of the problem class. + newton_tol : float + Tolerance for the nonlinear solver. + + Attributes + ---------- + work_counters : WorkCounter + Counts the work, here the number of function calls during the nonlinear solve is logged. """ dtype_u = mesh dtype_f = mesh def __init__(self, nvars, newton_tol): - """ - Initialization routine - - Args: - problem_params (dict): custom parameters for the example - dtype_u: mesh data type (will be passed parent class) - dtype_f: mesh data type (will be passed parent class) - """ + """Initialization routine""" super().__init__((nvars, None, np.dtype('float64'))) self._makeAttributeAndRegister('nvars', 'newton_tol', localVars=locals(), readOnly=True) + + self.work_counters['newton'] = WorkCounter() + + def solve_system(self, impl_sys, u0, t): + r""" + Solver for nonlinear implicit system (defined in sweeper). + + Parameters + ---------- + impl_sys : callable + Implicit system to be solved. + u0 : dtype_u + Initial guess for solver. + t : float + Current time :math:`t`. + + Returns + ------- + me : dtype_u + Numerical solution. + """ + + me = self.dtype_u(self.init) + opt = root( + impl_sys, + u0, + method='hybr', + tol=self.newton_tol, + ) + me[:] = opt.x + self.work_counters['newton'].niter += opt.nfev + return me diff --git a/pySDC/projects/DAE/problems/simple_DAE.py b/pySDC/projects/DAE/problems/simple_DAE.py index 859dbf7c5f..21ff17be26 100644 --- a/pySDC/projects/DAE/problems/simple_DAE.py +++ b/pySDC/projects/DAE/problems/simple_DAE.py @@ -1,7 +1,6 @@ import warnings import numpy as np from scipy.interpolate import interp1d -from scipy.optimize import root from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae from pySDC.core.Problem import WorkCounter @@ -38,8 +37,7 @@ class pendulum_2d(ptype_dae): ---------- work_counters : WorkCounter Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in - ``work_counters['newton']``. + ``work_counters['rhs']``. References ---------- @@ -57,7 +55,6 @@ def __init__(self, nvars, newton_tol): # solution = data[:, 1:] # self.u_ref = interp1d(t, solution, kind='cubic', axis=0, fill_value='extrapolate') self.t_end = 0.0 - self.work_counters['newton'] = WorkCounter() self.work_counters['rhs'] = WorkCounter() def eval_f(self, u, du, t): @@ -86,36 +83,6 @@ def eval_f(self, u, du, t): self.work_counters['rhs']() return f - def solve_system(self, impl_sys, u0, t): - r""" - Solver for nonlinear implicit system (defined in sweeper). - - Parameters - ---------- - impl_sys : callable - Implicit system to be solved. - u0 : dtype_u - Initial guess for solver. - t : float - Current time :math:`t`. - - Returns - ------- - me : dtype_u - Numerical solution. - """ - - me = self.dtype_u(self.init) - opt = root( - impl_sys, - u0, - method='hybr', - tol=self.newton_tol, - ) - me[:] = opt.x - self.work_counters['newton'].niter += opt.nfev - return me - def u_exact(self, t): """ Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical reference solution. @@ -177,8 +144,7 @@ class simple_dae_1(ptype_dae): ---------- work_counters : WorkCounter Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in - ``work_counters['newton']``. + ``work_counters['rhs']``. References ---------- @@ -190,7 +156,6 @@ def __init__(self, nvars, newton_tol): """Initialization routine""" super().__init__(nvars, newton_tol) - self.work_counters['newton'] = WorkCounter() self.work_counters['rhs'] = WorkCounter() def eval_f(self, u, du, t): @@ -222,36 +187,6 @@ def eval_f(self, u, du, t): self.work_counters['rhs']() return f - def solve_system(self, impl_sys, u0, t): - r""" - Solver for nonlinear implicit system (defined in sweeper). - - Parameters - ---------- - impl_sys : callable - Implicit system to be solved. - u0 : dtype_u - Initial guess for solver. - t : float - Current time :math:`t`. - - Returns - ------- - me : dtype_u - Numerical solution. - """ - - me = self.dtype_u(self.init) - opt = root( - impl_sys, - u0, - method='hybr', - tol=self.newton_tol, - ) - me[:] = opt.x - self.work_counters['newton'].niter += opt.nfev - return me - def u_exact(self, t): """ Routine for the exact solution. @@ -297,8 +232,7 @@ class problematic_f(ptype_dae): Specific parameter of the problem. work_counters : WorkCounter Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in - ``work_counters['newton']``. + ``work_counters['rhs']`` References ---------- @@ -312,7 +246,6 @@ def __init__(self, nvars, newton_tol, eta=1): self._makeAttributeAndRegister('eta', localVars=locals()) self.work_counters['rhs'] = WorkCounter() - self.work_counters['newton'] = WorkCounter() def eval_f(self, u, du, t): r""" @@ -340,36 +273,6 @@ def eval_f(self, u, du, t): self.work_counters['rhs']() return f - def solve_system(self, impl_sys, u0, t): - r""" - Solver for nonlinear implicit system (defined in sweeper). - - Parameters - ---------- - impl_sys : callable - Implicit system to be solved. - u0 : dtype_u - Initial guess for solver. - t : float - Current time :math:`t`. - - Returns - ------- - me : dtype_u - Numerical solution. - """ - - me = self.dtype_u(self.init) - opt = root( - impl_sys, - u0, - method='hybr', - tol=self.newton_tol, - ) - me[:] = opt.x - self.work_counters['newton'].niter += opt.nfev - return me - def u_exact(self, t): """ Routine for the exact solution. diff --git a/pySDC/projects/DAE/problems/synchronous_machine.py b/pySDC/projects/DAE/problems/synchronous_machine.py index 1f99f5644a..a3309ab369 100644 --- a/pySDC/projects/DAE/problems/synchronous_machine.py +++ b/pySDC/projects/DAE/problems/synchronous_machine.py @@ -1,7 +1,6 @@ import numpy as np import warnings from scipy.interpolate import interp1d -from scipy.optimize import root from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae from pySDC.core.Problem import WorkCounter @@ -147,8 +146,7 @@ class synchronous_machine_infinite_bus(ptype_dae): Defines the mechanical torque applied to the rotor shaft. work_counters : WorkCounter Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in - ``work_counters['newton']``. + ``work_counters['rhs']``. References ---------- @@ -191,7 +189,6 @@ def __init__(self, nvars, newton_tol): self.T_m = 0.854 self.work_counters['rhs'] = WorkCounter() - self.work_counters['newton'] = WorkCounter() def eval_f(self, u, du, t): r""" @@ -269,36 +266,6 @@ def eval_f(self, u, du, t): self.work_counters['rhs']() return f - def solve_system(self, impl_sys, u0, t): - r""" - Solver for nonlinear implicit system (defined in sweeper). - - Parameters - ---------- - impl_sys : callable - Implicit system to be solved. - u0 : dtype_u - Initial guess for solver. - t : float - Current time :math:`t`. - - Returns - ------- - me : dtype_u - Numerical solution. - """ - - me = self.dtype_u(self.init) - opt = root( - impl_sys, - u0, - method='hybr', - tol=self.newton_tol, - ) - me[:] = opt.x - self.work_counters['newton'].niter += opt.nfev - return me - def u_exact(self, t): """ Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical reference solution. diff --git a/pySDC/projects/DAE/problems/transistor_amplifier.py b/pySDC/projects/DAE/problems/transistor_amplifier.py index b160be8ef1..039678ebbe 100644 --- a/pySDC/projects/DAE/problems/transistor_amplifier.py +++ b/pySDC/projects/DAE/problems/transistor_amplifier.py @@ -1,6 +1,5 @@ import warnings import numpy as np -from scipy.optimize import root from scipy.interpolate import interp1d from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae @@ -58,8 +57,7 @@ class one_transistor_amplifier(ptype_dae): The end time at which the reference solution is determined. work_counters : WorkCounter Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in - ``work_counters['newton']``. + ``work_counters['rhs']``. References ---------- @@ -79,7 +77,6 @@ def __init__(self, nvars, newton_tol): self.t_end = 0.0 self.work_counters['rhs'] = WorkCounter() - self.work_counters['newton'] = WorkCounter() def eval_f(self, u, du, t): r""" @@ -116,36 +113,6 @@ def eval_f(self, u, du, t): self.work_counters['rhs']() return f - def solve_system(self, impl_sys, u0, t): - r""" - Solver for nonlinear implicit system (defined in sweeper). - - Parameters - ---------- - impl_sys : callable - Implicit system to be solved. - u0 : dtype_u - Initial guess for solver. - t : float - Current time :math:`t`. - - Returns - ------- - me : dtype_u - Numerical solution. - """ - - me = self.dtype_u(self.init) - opt = root( - impl_sys, - u0, - method='hybr', - tol=self.newton_tol, - ) - me[:] = opt.x - self.work_counters['newton'].niter += opt.nfev - return me - def u_exact(self, t): """ Approximation of the exact solution generated by spline interpolation of an extremely accurate numerical @@ -229,8 +196,7 @@ class two_transistor_amplifier(ptype_dae): The end time at which the reference solution is determined. work_counters : WorkCounter Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``, the number of function calls of the Newton-like solver is stored in - ``work_counters['newton']``. + ``work_counters['rhs']``. References ---------- @@ -250,7 +216,6 @@ def __init__(self, nvars, newton_tol): self.t_end = 0.0 self.work_counters['rhs'] = WorkCounter() - self.work_counters['newton'] = WorkCounter() def eval_f(self, u, du, t): r""" @@ -290,36 +255,6 @@ def eval_f(self, u, du, t): self.work_counters['rhs']() return f - def solve_system(self, impl_sys, u0, t): - r""" - Solver for nonlinear implicit system (defined in sweeper). - - Parameters - ---------- - impl_sys : callable - Implicit system to be solved. - u0 : dtype_u - Initial guess for solver. - t : float - Current time :math:`t`. - - Returns - ------- - me : dtype_u - Numerical solution. - """ - - me = self.dtype_u(self.init) - opt = root( - impl_sys, - u0, - method='hybr', - tol=self.newton_tol, - ) - me[:] = opt.x - self.work_counters['newton'].niter += opt.nfev - return me - def u_exact(self, t): """ Dummy exact solution that should only be used to get initial conditions for the problem. This makes From 6fd2a876ef8c86bb088176c7174ddd268da2cb31 Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 31 Oct 2023 08:45:18 +0100 Subject: [PATCH 4/6] Used possibility to update documentation :D --- .../DAE/sweepers/fully_implicit_DAE.py | 108 ++++++++++++------ 1 file changed, 71 insertions(+), 37 deletions(-) diff --git a/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py b/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py index 414fde4545..e073eabde0 100644 --- a/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py +++ b/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py @@ -6,15 +6,43 @@ class fully_implicit_DAE(sweeper): - """ - Custom sweeper class, implements Sweeper.py + r""" + Custom sweeper class to implement the fully-implicit SDC for solving DAEs. It solves fully-implicit DAE problems + of the form + + .. math:: + F(t, u, u') = 0. + + It solves a collocation problem of the form + + .. math:: + F(\tau, \vec{U}_0 + \Delta t (\mathbf{Q} \otimes \mathbf{I}_n) \vec{U}, \vec{U}) = 0, + + where - Sweeper to solve first order differential equations in fully implicit form - Primarily implemented to be used with differential algebraic equations - Based on the concepts outlined in "Arbitrary order Krylov deferred correction methods for differential algebraic equations" by Huang et al. + - :math:`\tau=(\tau_1,..,\tau_M) in \mathbb{R}^M` the vector of collocation nodes, + - :math:`\vec{U}_0 = (u_0,..,u_0) \in \mathbb{R}^{Mn}` the vector of initial condition spread to each node, + - spectral integration matrix :math:`\mathbf{Q} \in \mathbb{R}^{M \times M}`, + - :math:`\vec{U}=(U_1,..,U_M) \in \mathbb{R}^{Mn}` the vector of unknown derivatives + :math:`U_m \approx U(\tau_m) = u'(\tau_m) \in \mathbb{R}^n`, + - and identity matrix :math:`\mathbf{I}_n \in \mathbb{R}^{n \times n}`. - Attributes: - QI: implicit Euler integration matrix + The construction of this sweeper is based on the concepts outlined in [1]_. + + Parameters + ---------- + params : dict + Parameters passed to the sweeper. + + Attributes + ---------- + QI : np.2darray + Implicit Euler integration matrix. + + References + ---------- + .. [1] J. Huang, J. Jun, M. L. Minion. Arbitrary order Krylov deferred correction methods for differential algebraic equation. + J. Comput. Phys. Vol. 221 No. 2 (2007). """ def __init__(self, params): @@ -40,12 +68,15 @@ def __init__(self, params): # TODO: hijacking this function to return solution from its gradient i.e. fundamental theorem of calculus. # This works well since (ab)using level.f to store the gradient. Might need to change this for release? def integrate(self): - """ - Returns the solution by integrating its gradient (fundamental theorem of calculus) - Note that level.f stores the gradient values in the fully implicit case, rather than the evaluation of the rhs as in the ODE case - - Returns: - list of dtype_u: containing the integral as values + r""" + Returns the solution by integrating its gradient (fundamental theorem of calculus) at each collocation node. + Note that ``level.f`` stores the gradient values in the fully implicit case, rather than the evaluation of + the right-hand side as in the ODE case. + + Returns + ------- + me : list of lists + Integral of the gradient at each collocation node. """ # get current level and problem description @@ -65,11 +96,9 @@ def integrate(self): return me def update_nodes(self): - """ - Update the u- and f-values at the collocation nodes -> corresponds to a single iteration of the preconditioned Richardson iteration in "ordinary" SDC - - Returns: - None + r""" + Updates values of ``u`` and ``f`` at collocation nodes. This correspond to a single iteration of the + preconditioned Richardson iteration in **"ordinary"** SDC. """ # get current level and problem description @@ -154,12 +183,16 @@ def implSystem(params): return None def predict(self): - """ - Predictor to fill values at nodes before first sweep + r""" + Predictor to fill values at nodes before first sweep. It can decides whether the - Default prediction for the sweepers, only copies the values to all collocation nodes - This function overrides the base implementation by always initialising level.f to zero - This is necessary since level.f stores the solution derivative in the fully implicit case, which is not initially known + - initial condition is spread to each node ('initial_guess' = 'spread'), + - zero values are spread to each node ('initial_guess' = 'zero'), + - or random values are spread to each collocation node ('initial_guess' = 'random'). + + Default prediction for the sweepers, only copies the values to all collocation nodes. This function + overrides the base implementation by always initialising ``level.f`` to zero. This is necessary since + ``level.f`` stores the solution derivative in the fully implicit case, which is not initially known. """ # get current level and problem description L = self.level @@ -187,15 +220,18 @@ def predict(self): L.status.updated = True def compute_residual(self, stage=None): - """ - Overrides the base implementation - Uses the absolute value of the implicit function ||F(u', u, t)|| as the residual + r""" + Uses the absolute value of the DAE system - Args: - stage (str): The current stage of the step the level belongs to + .. math:: + ||F(t, u, u')|| + + for computing the residual. - Returns: - None + Parameters + ---------- + stage : str, optional + The current stage of the step the level belongs to. """ # get current level and problem description @@ -239,13 +275,11 @@ def compute_residual(self, stage=None): 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 + r""" + Computes the solution ``u`` at the right-hand point. For ``quad_type='RADAU-LEFT'`` a collocation update + has to be done, which is the full evaluation of the Picard formulation. In cases of + ``quad_type='RADAU-RIGHT'`` or ``quad_type='LOBATTO'`` the value at last collocation node is the new value + for the next step. """ # get current level and problem description From 076607f0ab1066e7ab92d15f0e578a2890e7ee66 Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 31 Oct 2023 08:57:57 +0100 Subject: [PATCH 5/6] Moved work_counters for rhs also to ProblemDAE --- pySDC/projects/DAE/misc/ProblemDAE.py | 5 +++- pySDC/projects/DAE/problems/simple_DAE.py | 24 ------------------- .../DAE/problems/synchronous_machine.py | 6 ----- .../DAE/problems/transistor_amplifier.py | 11 --------- 4 files changed, 4 insertions(+), 42 deletions(-) diff --git a/pySDC/projects/DAE/misc/ProblemDAE.py b/pySDC/projects/DAE/misc/ProblemDAE.py index 91cc6bc148..08cc16a3a7 100644 --- a/pySDC/projects/DAE/misc/ProblemDAE.py +++ b/pySDC/projects/DAE/misc/ProblemDAE.py @@ -20,7 +20,9 @@ class ptype_dae(ptype): Attributes ---------- work_counters : WorkCounter - Counts the work, here the number of function calls during the nonlinear solve is logged. + Counts the work, here the number of function calls during the nonlinear solve is logged and stored + in work_counters['newton']. The number of each function class of the right-hand side is then stored + in work_counters['rhs'] """ dtype_u = mesh @@ -32,6 +34,7 @@ def __init__(self, nvars, newton_tol): self._makeAttributeAndRegister('nvars', 'newton_tol', localVars=locals(), readOnly=True) self.work_counters['newton'] = WorkCounter() + self.work_counters['rhs'] = WorkCounter() def solve_system(self, impl_sys, u0, t): r""" diff --git a/pySDC/projects/DAE/problems/simple_DAE.py b/pySDC/projects/DAE/problems/simple_DAE.py index 21ff17be26..343128624a 100644 --- a/pySDC/projects/DAE/problems/simple_DAE.py +++ b/pySDC/projects/DAE/problems/simple_DAE.py @@ -3,7 +3,6 @@ from scipy.interpolate import interp1d from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae -from pySDC.core.Problem import WorkCounter class pendulum_2d(ptype_dae): @@ -33,11 +32,6 @@ class pendulum_2d(ptype_dae): ---------- t_end: float The end time at which the reference solution is determined. - Attributes - ---------- - work_counters : WorkCounter - Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``. References ---------- @@ -55,7 +49,6 @@ def __init__(self, nvars, newton_tol): # solution = data[:, 1:] # self.u_ref = interp1d(t, solution, kind='cubic', axis=0, fill_value='extrapolate') self.t_end = 0.0 - self.work_counters['rhs'] = WorkCounter() def eval_f(self, u, du, t): r""" @@ -140,24 +133,12 @@ class simple_dae_1(ptype_dae): newton_tol : float Tolerance for Newton solver. - Attributes - ---------- - work_counters : WorkCounter - Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``. - References ---------- .. [1] U. Ascher, L. R. Petzold. Computer method for ordinary differential equations and differential-algebraic equations. Society for Industrial and Applied Mathematics (1998). """ - def __init__(self, nvars, newton_tol): - """Initialization routine""" - super().__init__(nvars, newton_tol) - - self.work_counters['rhs'] = WorkCounter() - def eval_f(self, u, du, t): r""" Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. @@ -230,9 +211,6 @@ class problematic_f(ptype_dae): ---------- eta : float Specific parameter of the problem. - work_counters : WorkCounter - Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']`` References ---------- @@ -245,8 +223,6 @@ def __init__(self, nvars, newton_tol, eta=1): super().__init__(nvars, newton_tol) self._makeAttributeAndRegister('eta', localVars=locals()) - self.work_counters['rhs'] = WorkCounter() - def eval_f(self, u, du, t): r""" Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. diff --git a/pySDC/projects/DAE/problems/synchronous_machine.py b/pySDC/projects/DAE/problems/synchronous_machine.py index a3309ab369..5510f0dd46 100644 --- a/pySDC/projects/DAE/problems/synchronous_machine.py +++ b/pySDC/projects/DAE/problems/synchronous_machine.py @@ -3,7 +3,6 @@ from scipy.interpolate import interp1d from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae -from pySDC.core.Problem import WorkCounter from pySDC.implementations.datatype_classes.mesh import mesh @@ -144,9 +143,6 @@ class synchronous_machine_infinite_bus(ptype_dae): Voltage at the field winding. T_m: float Defines the mechanical torque applied to the rotor shaft. - work_counters : WorkCounter - Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``. References ---------- @@ -188,8 +184,6 @@ def __init__(self, nvars, newton_tol): self.v_F = 8.736809687330562e-4 self.T_m = 0.854 - self.work_counters['rhs'] = WorkCounter() - def eval_f(self, u, du, t): r""" Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. diff --git a/pySDC/projects/DAE/problems/transistor_amplifier.py b/pySDC/projects/DAE/problems/transistor_amplifier.py index 039678ebbe..ae2adeb9cc 100644 --- a/pySDC/projects/DAE/problems/transistor_amplifier.py +++ b/pySDC/projects/DAE/problems/transistor_amplifier.py @@ -3,7 +3,6 @@ from scipy.interpolate import interp1d from pySDC.projects.DAE.misc.ProblemDAE import ptype_dae -from pySDC.core.Problem import WorkCounter # Helper function @@ -55,9 +54,6 @@ class one_transistor_amplifier(ptype_dae): ---------- t_end: float The end time at which the reference solution is determined. - work_counters : WorkCounter - Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``. References ---------- @@ -76,8 +72,6 @@ def __init__(self, nvars, newton_tol): # self.u_ref = interp1d(x, y, kind='cubic', axis=0, fill_value='extrapolate') self.t_end = 0.0 - self.work_counters['rhs'] = WorkCounter() - def eval_f(self, u, du, t): r""" Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. @@ -194,9 +188,6 @@ class two_transistor_amplifier(ptype_dae): ---------- t_end: float The end time at which the reference solution is determined. - work_counters : WorkCounter - Counts the work, i.e., number of function calls of right-hand side is called and stored in - ``work_counters['rhs']``. References ---------- @@ -215,8 +206,6 @@ def __init__(self, nvars, newton_tol): # self.u_ref = interp1d(x, y, kind='cubic', axis=0, fill_value='extrapolate') self.t_end = 0.0 - self.work_counters['rhs'] = WorkCounter() - def eval_f(self, u, du, t): r""" Routine to evaluate the implicit representation of the problem, i.e., :math:`F(u, u', t)`. From ff5805f8da3c37bb760131ae1259596d57beaa0c Mon Sep 17 00:00:00 2001 From: lisawim Date: Tue, 31 Oct 2023 10:42:52 +0100 Subject: [PATCH 6/6] Inheritance from generic_implicit --- .../DAE/sweepers/fully_implicit_DAE.py | 61 +------------------ 1 file changed, 3 insertions(+), 58 deletions(-) diff --git a/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py b/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py index e073eabde0..aab811ffb9 100644 --- a/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py +++ b/pySDC/projects/DAE/sweepers/fully_implicit_DAE.py @@ -2,10 +2,10 @@ from scipy import optimize from pySDC.core.Errors import ParameterError -from pySDC.core.Sweeper import sweeper +from pySDC.implementations.sweeper_classes.generic_implicit import generic_implicit -class fully_implicit_DAE(sweeper): +class fully_implicit_DAE(generic_implicit): r""" Custom sweeper class to implement the fully-implicit SDC for solving DAEs. It solves fully-implicit DAE problems of the form @@ -65,36 +65,6 @@ def __init__(self, params): self.QI = self.get_Qdelta_implicit(coll=self.coll, qd_type=self.params.QI) - # TODO: hijacking this function to return solution from its gradient i.e. fundamental theorem of calculus. - # This works well since (ab)using level.f to store the gradient. Might need to change this for release? - def integrate(self): - r""" - Returns the solution by integrating its gradient (fundamental theorem of calculus) at each collocation node. - Note that ``level.f`` stores the gradient values in the fully implicit case, rather than the evaluation of - the right-hand side as in the ODE case. - - Returns - ------- - me : list of lists - Integral of the gradient at each collocation node. - """ - - # get current level and problem description - L = self.level - P = L.prob - M = self.coll.num_nodes - - me = [] - - # integrate gradient over all collocation nodes - for m in range(1, M + 1): - # new instance of dtype_u, initialize values with 0 - me.append(P.dtype_u(P.init, val=0.0)) - for j in range(1, M + 1): - me[-1] += L.dt * self.coll.Qmat[m, j] * L.f[j] - - return me - def update_nodes(self): r""" Updates values of ``u`` and ``f`` at collocation nodes. This correspond to a single iteration of the @@ -184,7 +154,7 @@ def implSystem(params): def predict(self): r""" - Predictor to fill values at nodes before first sweep. It can decides whether the + Predictor to fill values at nodes before first sweep. It can decide whether the - initial condition is spread to each node ('initial_guess' = 'spread'), - zero values are spread to each node ('initial_guess' = 'zero'), @@ -204,7 +174,6 @@ def predict(self): if self.params.initial_guess == 'spread': L.u[m] = P.dtype_u(L.u[0]) L.f[m] = P.dtype_f(init=P.init, val=0.0) - # start with zero everywhere elif self.params.initial_guess == 'zero': L.u[m] = P.dtype_u(init=P.init, val=0.0) L.f[m] = P.dtype_f(init=P.init, val=0.0) @@ -273,27 +242,3 @@ def compute_residual(self, stage=None): L.status.updated = False return None - - def compute_end_point(self): - r""" - Computes the solution ``u`` at the right-hand point. For ``quad_type='RADAU-LEFT'`` a collocation update - has to be done, which is the full evaluation of the Picard formulation. In cases of - ``quad_type='RADAU-RIGHT'`` or ``quad_type='LOBATTO'`` the value at last collocation node is the new value - for the next step. - """ - - # get current level and problem description - L = self.level - P = L.prob - - # 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: - # a copy is sufficient - L.uend = P.dtype_u(L.u[-1]) - else: - # start with u0 and add integral over the full interval (using coll.weights) - L.uend = P.dtype_u(L.u[0]) - for m in range(self.coll.num_nodes): - L.uend += L.dt * self.coll.weights[m] * L.f[m + 1] - - return None