Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Added solve_system to all DAE problem classes #371

Merged
merged 6 commits into from
Nov 3, 2023
+143 −103
Merged

Added solve_system to all DAE problem classes #371

Changes from 1 commit
Commits
Show all changes
6 commits
Select commit Hold shift + click to select a range
6ae4779
Added solve_system method to all DAE classes
lisawim Oct 31, 2023
c325087
Removed comment
lisawim Oct 31, 2023
76efa6a
Moved solve_system to generic ProblemDAE class
lisawim Oct 31, 2023
6fd2a87
Used possibility to update documentation :D
lisawim Oct 31, 2023
076607f
Moved work_counters for rhs also to ProblemDAE
lisawim Oct 31, 2023
ff5805f
Inheritance from generic_implicit
lisawim Oct 31, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Prev Previous commit
Next Next commit
Moved solve_system to generic ProblemDAE class
@lisawim
lisawim committed Oct 31, 2023
commit 76efa6aaf9a6e022d831f86867ed5825085eca31
61 changes: 50 additions & 11 deletions pySDC/projects/DAE/misc/ProblemDAE.py
View file
Original file line number Diff line number Diff line change
@@ -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
103 changes: 3 additions & 100 deletions pySDC/projects/DAE/problems/simple_DAE.py
View file
Original file line number Diff line number Diff line change
@@ -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()
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks like at least this one could go to the generic parent class, too?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

True! I added this for all problem classes. I'll do that.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks better and more convenient (like before), thanks a lot!


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.
35 changes: 1 addition & 34 deletions pySDC/projects/DAE/problems/synchronous_machine.py
View file
Original file line number Diff line number Diff line change
@@ -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.
69 changes: 2 additions & 67 deletions pySDC/projects/DAE/problems/transistor_amplifier.py
View file
Original file line number Diff line number Diff line change
@@ -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