Rename PreconditionerBase to SystemJacobian

This commit is contained in:
Ray Speth 2025-01-01 21:53:51 -05:00 committed by Ingmar Schoegl
parent a3b1765714
commit 020c690c7e
22 changed files with 369 additions and 328 deletions

View File

@ -1,6 +1,6 @@
/**
* @file AdaptivePreconditioner.h Declarations for the class
* AdaptivePreconditioner which is a child class of PreconditionerBase
* AdaptivePreconditioner which is a child class of SystemJacobian
* for preconditioners used by sundials
*/
@ -10,7 +10,7 @@
#ifndef ADAPTIVEPRECONDITIONER_H
#define ADAPTIVEPRECONDITIONER_H
#include "cantera/numerics/PreconditionerBase.h"
#include "cantera/numerics/SystemJacobian.h"
#include "cantera/numerics/eigen_sparse.h"
#include "cantera/base/global.h"
#include <iostream>
@ -23,7 +23,7 @@ namespace Cantera
//! the preconditioner by a threshold value. It also neglects pressure
//! dependence and third body contributions in its formation and has a
//! finite difference approximation for temperature.
class AdaptivePreconditioner : public PreconditionerBase
class AdaptivePreconditioner : public SystemJacobian
{
public:
AdaptivePreconditioner();

View File

@ -8,7 +8,7 @@
#ifndef CT_INTEGRATOR_H
#define CT_INTEGRATOR_H
#include "FuncEval.h"
#include "PreconditionerBase.h"
#include "SystemJacobian.h"
#include "cantera/base/global.h"
#include "cantera/base/AnyMap.h"
@ -91,7 +91,7 @@ public:
/*!
* @param preconditioner preconditioner object used for the linear solver
*/
virtual void setPreconditioner(shared_ptr<PreconditionerBase> preconditioner) {
virtual void setPreconditioner(shared_ptr<SystemJacobian> preconditioner) {
m_preconditioner = preconditioner;
if (preconditioner->preconditionerSide() == "none") {
m_prec_side = PreconditionerSide::NO_PRECONDITION;
@ -123,7 +123,7 @@ public:
}
//! Return preconditioner reference to object
virtual shared_ptr<PreconditionerBase> preconditioner() {
virtual shared_ptr<SystemJacobian> preconditioner() {
return m_preconditioner;
}
@ -296,7 +296,7 @@ protected:
//! Pointer to preconditioner object used in integration which is
//! set by setPreconditioner and initialized inside of
//! ReactorNet::initialize()
shared_ptr<PreconditionerBase> m_preconditioner;
shared_ptr<SystemJacobian> m_preconditioner;
//! Type of preconditioning used in applyOptions
PreconditionerSide m_prec_side = PreconditionerSide::NO_PRECONDITION;
// methods for DAE solvers

View File

@ -1,195 +1,15 @@
/**
* @file PreconditionerBase.h Declarations for the class
* PreconditionerBase which is a virtual base class for
* preconditioning systems.
*/
//! @file PreconditionerBase.h
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#ifndef PRECONDITIONERBASE_H
#define PRECONDITIONERBASE_H
#ifndef PRECONDITIONERBASE_FACTORY_H
#define PRECONDITIONERBASE_FACTORY_H
#include "cantera/base/ctexceptions.h"
#pragma message("warning: PreconditionerBase.h and class PreconditionerBase are " \
"deprecated and will be removed after Cantera 3.2. Use " \
"SystemJacobian.h and class SystemJacobian.")
namespace Cantera
{
#include "SystemJacobian.h"
/**
* Specifies the side of the system on which the preconditioner is applied. Not all
* methods are supported by all integrators.
*/
enum class PreconditionerSide {
NO_PRECONDITION, //! No preconditioning
LEFT_PRECONDITION, //! Left side preconditioning
RIGHT_PRECONDITION, //! Right side preconditioning
BOTH_PRECONDITION //! Left and right side preconditioning
};
//! PreconditionerBase serves as an abstract type to extend different preconditioners
class PreconditionerBase
{
public:
PreconditionerBase() {}
virtual ~PreconditionerBase() {}
virtual const string type() const = 0;
//! Set a value at the specified row and column of the jacobian triplet vector
//! @param row row in the jacobian matrix
//! @param col column in the jacobian matrix
//! @param value value of the element at row and col
virtual void setValue(size_t row, size_t col, double value) {
throw NotImplementedError("PreconditionerBase::setValue");
}
//! Adjust the state vector based on the preconditioner, e.g., Adaptive
//! preconditioning uses a strictly positive composition when preconditioning which
//! is handled by this function
//! @param state a vector containing the state to be updated
virtual void stateAdjustment(vector<double>& state) {
throw NotImplementedError("PreconditionerBase::stateAdjustment");
}
//! Get preconditioner application side for CVODES
string preconditionerSide() const {
return m_precon_side;
}
virtual void setPreconditionerSide(const string& preconSide) {
m_precon_side = preconSide;
}
//! Update the diagonal terms in the Jacobian by using the transient mask.
//! @param rdt Reciprocal of the time step [1/s]
//! @param mask Mask for transient terms: 1 if transient, 0 if algebraic.
virtual void updateTransient(double rdt, int* mask) {
throw NotImplementedError("PreconditionerBase::updateTransient");
}
//! Solve a linear system Ax=b where A is the preconditioner
//! @param[in] stateSize length of the rhs and output vectors
//! @param[in] rhs_vector right hand side vector used in linear system
//! @param[out] output output vector for solution
virtual void solve(const size_t stateSize, double* rhs_vector, double* output) {
throw NotImplementedError("PreconditionerBase::solve");
};
//! Perform preconditioner specific post-reactor setup operations such as factorize.
//! @deprecated To be removed after %Cantera 3.2. Use updatePreconditioner() and
//! factorize() instead.
virtual void setup() {
throw NotImplementedError("PreconditionerBase::setup");
};
//! Reset preconditioner parameters as needed
virtual void reset() {
throw NotImplementedError("PreconditionerBase::reset");
};
//! Called during setup for any processes that need
//! to be completed prior to setup functions used in sundials.
//! @param networkSize the number of variables in the associated reactor network
virtual void initialize(size_t networkSize) {
throw NotImplementedError("PreconditionerBase::initialize");
};
//! Used to provide system bandwidth for implementations that use banded matrix
//! storage. Ignored if not needed.
virtual void setBandwidth(size_t bw) {}
//! Print preconditioner contents
virtual void printPreconditioner() {
throw NotImplementedError("PreconditionerBase::printPreconditioner");
};
//! Transform Jacobian vector and write into
//! preconditioner, P = (I - gamma * J)
virtual void updatePreconditioner() {
throw NotImplementedError("PreconditionerBase::updatePreconditioner");
}
//! Set gamma used in preconditioning
//! @param gamma used in M = I - gamma*J
virtual void setGamma(double gamma) {
m_gamma = gamma;
};
//! Get gamma used in preconditioning
virtual double gamma() {
return m_gamma;
};
//! Set the absolute tolerance in the solver outside of the network initialization
//! @param atol the specified tolerance
virtual void setAbsoluteTolerance(double atol) {
m_atol = atol;
}
//! Get latest status of linear solver. Zero indicates success. Meaning of non-zero
//! values is solver dependent.
virtual int info() const {
throw NotImplementedError("PreconditionerBase::info");
}
//! Elapsed CPU time spent computing the Jacobian elements.
double elapsedTime() const {
return m_elapsed;
}
void updateElapsed(double evalTime) {
m_elapsed += evalTime;
}
//! Number of Jacobian evaluations.
int nEvals() const {
return m_nevals;
}
void incrementEvals() {
m_nevals++;
}
//! Number of times 'incrementAge' has been called since the last evaluation
int age() const {
return m_age;
}
//! Increment the Jacobian age.
void incrementAge() {
m_age++;
}
//! Set the Jacobian age.
void setAge(int age) {
m_age = age;
}
protected:
//! Factorize the system matrix. This method should be called at the end of
//! implementations of updatePreconditioner() and updateTransient().
virtual void factorize() {
throw NotImplementedError("preconditionerBase::factorize");
}
//! Dimension of the preconditioner
size_t m_dim;
//! gamma value used in M = I - gamma*J
double m_gamma = 1.0;
//! bool saying whether or not the preconditioner is initialized
bool m_init = false;
//! Absolute tolerance of the ODE solver
double m_atol = 0;
double m_elapsed = 0.0; //!< Elapsed CPU time taken to compute the Jacobian
int m_nevals = 0; //!< Number of Jacobian evaluations.
int m_age = 100000; //!< Age of the Jacobian (times incrementAge() has been called)
string m_precon_side = "none";
};
}
#endif

View File

@ -3,34 +3,13 @@
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#ifndef PRECONDITIONER_FACTORY_H
#define PRECONDITIONER_FACTORY_H
#ifndef PRECONDITIONERFACTORY_FACTORY_H
#define PRECONDITIONERFACTORY_FACTORY_H
#include "cantera/base/FactoryBase.h"
#pragma message("warning: PreconditionerFactory.h and class PreconditionerFactory " \
"are deprecated and will be removed after Cantera 3.2. Use " \
"SystemJacobianFactory.h and class SystemJacobianFactory.")
namespace Cantera
{
class PreconditionerBase;
//! Factory class to create preconditioner objects
class PreconditionerFactory : public Factory<PreconditionerBase>
{
public:
static PreconditionerFactory* factory();
//! Delete preconditioner factory
void deleteFactory() override;
private:
static PreconditionerFactory* s_factory;
static std::mutex precon_mutex;
PreconditionerFactory();
};
//! Create a Preconditioner object of the specified type
shared_ptr<PreconditionerBase> newPreconditioner(const string& precon);
}
#include "SystemJacobianFactory.h"
#endif

View File

@ -0,0 +1,202 @@
//! @file SystemJacobian.h Declarations for class SystemJacobian
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#ifndef SYSTEMJACOBIAN_H
#define SYSTEMJACOBIAN_H
#include "cantera/base/ctexceptions.h"
#include "cantera/base/global.h"
namespace Cantera
{
/**
* Specifies the side of the system on which the preconditioner is applied. Not all
* methods are supported by all integrators.
*/
enum class PreconditionerSide {
NO_PRECONDITION, //! No preconditioning
LEFT_PRECONDITION, //! Left side preconditioning
RIGHT_PRECONDITION, //! Right side preconditioning
BOTH_PRECONDITION //! Left and right side preconditioning
};
//! Abstract base class representing Jacobian matrices and preconditioners used in
//! nonlinear solvers.
class SystemJacobian
{
public:
SystemJacobian() {}
virtual ~SystemJacobian() {}
virtual const string type() const = 0;
//! Set a value at the specified row and column of the jacobian triplet vector
//! @param row row in the jacobian matrix
//! @param col column in the jacobian matrix
//! @param value value of the element at row and col
virtual void setValue(size_t row, size_t col, double value) {
throw NotImplementedError("SystemJacobian::setValue");
}
//! Adjust the state vector based on the preconditioner, e.g., Adaptive
//! preconditioning uses a strictly positive composition when preconditioning which
//! is handled by this function
//! @param state a vector containing the state to be updated
virtual void stateAdjustment(vector<double>& state) {
throw NotImplementedError("SystemJacobian::stateAdjustment");
}
//! Get preconditioner application side for CVODES
string preconditionerSide() const {
return m_precon_side;
}
virtual void setPreconditionerSide(const string& preconSide) {
m_precon_side = preconSide;
}
//! Update the diagonal terms in the Jacobian by using the transient mask.
//! @param rdt Reciprocal of the time step [1/s]
//! @param mask Mask for transient terms: 1 if transient, 0 if algebraic.
virtual void updateTransient(double rdt, int* mask) {
throw NotImplementedError("SystemJacobian::updateTransient");
}
//! Solve a linear system Ax=b where A is the preconditioner
//! @param[in] stateSize length of the rhs and output vectors
//! @param[in] rhs_vector right hand side vector used in linear system
//! @param[out] output output vector for solution
virtual void solve(const size_t stateSize, double* rhs_vector, double* output) {
throw NotImplementedError("SystemJacobian::solve");
};
//! Perform preconditioner specific post-reactor setup operations such as factorize.
//! @deprecated To be removed after %Cantera 3.2. Use updatePreconditioner() and
//! factorize() instead.
virtual void setup() {
throw NotImplementedError("SystemJacobian::setup");
};
//! Reset parameters as needed
virtual void reset() {
throw NotImplementedError("SystemJacobian::reset");
};
//! Called during setup for any processes that need
//! to be completed prior to setup functions used in sundials.
//! @param networkSize the number of variables in the associated reactor network
virtual void initialize(size_t networkSize) {
throw NotImplementedError("SystemJacobian::initialize");
};
//! Used to provide system bandwidth for implementations that use banded matrix
//! storage. Ignored if not needed.
virtual void setBandwidth(size_t bw) {}
//! Print preconditioner contents
virtual void printPreconditioner() {
throw NotImplementedError("SystemJacobian::printPreconditioner");
};
//! Transform Jacobian vector and write into
//! preconditioner, P = (I - gamma * J)
virtual void updatePreconditioner() {
throw NotImplementedError("SystemJacobian::updatePreconditioner");
}
//! Set gamma used in preconditioning
//! @param gamma used in M = I - gamma*J
virtual void setGamma(double gamma) {
m_gamma = gamma;
};
//! Get gamma used in preconditioning
virtual double gamma() {
return m_gamma;
};
//! Set the absolute tolerance in the solver outside of the network initialization
//! @param atol the specified tolerance
virtual void setAbsoluteTolerance(double atol) {
m_atol = atol;
}
//! Get latest status of linear solver. Zero indicates success. Meaning of non-zero
//! values is solver dependent.
virtual int info() const {
throw NotImplementedError("SystemJacobian::info");
}
//! Elapsed CPU time spent computing the Jacobian elements.
double elapsedTime() const {
return m_elapsed;
}
void updateElapsed(double evalTime) {
m_elapsed += evalTime;
}
//! Number of Jacobian evaluations.
int nEvals() const {
return m_nevals;
}
void incrementEvals() {
m_nevals++;
}
//! Number of times 'incrementAge' has been called since the last evaluation
int age() const {
return m_age;
}
//! Increment the Jacobian age.
void incrementAge() {
m_age++;
}
//! Set the Jacobian age.
void setAge(int age) {
m_age = age;
}
protected:
//! Factorize the system matrix. This method should be called at the end of
//! implementations of updatePreconditioner() and updateTransient().
virtual void factorize() {
throw NotImplementedError("SystemJacobian::factorize");
}
//! Dimension of the system
size_t m_dim;
//! gamma value used in M = I - gamma*J
double m_gamma = 1.0;
//! bool saying whether or not the system is initialized
bool m_init = false;
//! Absolute tolerance of the ODE solver
double m_atol = 0;
double m_elapsed = 0.0; //!< Elapsed CPU time taken to compute the Jacobian
int m_nevals = 0; //!< Number of Jacobian evaluations.
int m_age = 100000; //!< Age of the Jacobian (times incrementAge() has been called)
string m_precon_side = "none";
};
//! @deprecated To be removed after Cantera 3.2. Renamed to SystemJacobian.
class PreconditionerBase : public SystemJacobian
{
PreconditionerBase() {
warn_deprecated("PreconditionerBase::PreconditionerBase",
"To be removed after Cantera 3.2. Renamed to SystemJacobian.");
}
};
}
#endif

View File

@ -0,0 +1,42 @@
//! @file SystemJacobianFactory.h
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#ifndef SYSTEMJACOBIANFACTORY_FACTORY_H
#define SYSTEMJACOBIANFACTORY_FACTORY_H
#include "cantera/base/FactoryBase.h"
#include "cantera/base/global.h"
namespace Cantera
{
class SystemJacobian;
//! Factory class to create Jacobian objects for use by linear solvers
class SystemJacobianFactory : public Factory<SystemJacobian>
{
public:
static SystemJacobianFactory* factory();
void deleteFactory() override;
private:
static SystemJacobianFactory* s_factory;
static std::mutex jac_mutex;
SystemJacobianFactory();
};
//! Create a SystemJacobian object of the specified type
shared_ptr<SystemJacobian> newSystemJacobian(const string& type);
//! @deprecated To be removed after %Cantera 3.2. Renamed to newSystemJacobian()
inline shared_ptr<SystemJacobian> newPreconditioner(const string& type) {
warn_deprecated("newPreconditioner",
"To be removed after Cantera 3.2. Renamed to newSystemJacobian");
return newSystemJacobian(type);
}
}
#endif

View File

@ -7,7 +7,7 @@
#define CT_MULTIJAC_H
#include "cantera/numerics/BandMatrix.h"
#include "cantera/numerics/PreconditionerBase.h"
#include "cantera/numerics/SystemJacobian.h"
#include "OneDim.h"
namespace Cantera
@ -21,7 +21,7 @@ namespace Cantera
* @ingroup onedUtilsGroup
* @ingroup derivGroup
*/
class MultiJac : public PreconditionerBase
class MultiJac : public SystemJacobian
{
public:
//! Constructor.

View File

@ -10,7 +10,7 @@
#include "Domain1D.h"
#include "MultiJac.h"
#include "cantera/numerics/PreconditionerBase.h"
#include "cantera/numerics/SystemJacobian.h"
namespace Cantera
{
@ -44,7 +44,7 @@ public:
//! @ingroup derivGroup
MultiJac& jacobian();
shared_ptr<PreconditionerBase> getJacobian() {
shared_ptr<SystemJacobian> getJacobian() {
return m_jac;
}
@ -53,10 +53,10 @@ public:
//! Set the linear solver used to hold the Jacobian matrix and solve linear systems
//! as part of each Newton iteration. The default is a direct, banded solver.
void setLinearSolver(shared_ptr<PreconditionerBase> solver);
void setLinearSolver(shared_ptr<SystemJacobian> solver);
//! Get the type of the linear solver being used.
shared_ptr<PreconditionerBase> linearSolver() const { return m_jac; }
shared_ptr<SystemJacobian> linearSolver() const { return m_jac; }
/**
* Solve F(x) = 0, where F(x) is the multi-domain residual function.
@ -386,7 +386,7 @@ protected:
shared_ptr<vector<double>> m_state; //!< Solution vector
shared_ptr<PreconditionerBase> m_jac; //!< Jacobian evaluator
shared_ptr<SystemJacobian> m_jac; //!< Jacobian evaluator
unique_ptr<MultiNewton> m_newt; //!< Newton iterator
double m_rdt = 0.0; //!< reciprocal of time step
bool m_jac_ok = false; //!< if true, Jacobian is current

View File

@ -15,7 +15,7 @@ namespace Cantera
class Array2D;
class Integrator;
class PreconditionerBase;
class SystemJacobian;
//! A class representing a network of connected reactors.
/*!
@ -44,7 +44,7 @@ public:
//! Set preconditioner used by the linear solver
//! @param preconditioner preconditioner object used for the linear solver
void setPreconditioner(shared_ptr<PreconditionerBase> preconditioner);
void setPreconditioner(shared_ptr<SystemJacobian> preconditioner);
//! Set the initial value of the independent variable (typically time).
//! Default = 0.0 s. Restarts integration from this value using the current mixture
@ -341,7 +341,7 @@ protected:
double m_rtolsens = 1.0e-4;
double m_atols = 1.0e-15;
double m_atolsens = 1.0e-6;
shared_ptr<PreconditionerBase> m_precon;
shared_ptr<SystemJacobian> m_precon;
string m_linearSolverType;
//! Maximum integrator internal timestep. Default of 0.0 means infinity.

View File

@ -10,7 +10,7 @@ from cantera.delegator cimport *
from cantera.func1 cimport *
from cantera.kinetics cimport *
from cantera.mixture cimport *
from cantera.preconditioners cimport *
from cantera.jacobians cimport *
from cantera.reaction cimport *
from cantera.reactionpath cimport *
from cantera.reactor cimport *

View File

@ -43,7 +43,7 @@ from .transport import *
from .units import *
from .yamlwriter import *
from .constants import *
from .preconditioners import *
from .jacobians import *
# Custom finder/loader no longer needed, so remove it
sys.meta_path.pop()

View File

@ -9,7 +9,7 @@ from .solutionbase cimport *
from .kinetics cimport *
from .func1 cimport *
from .thermo cimport *
from .preconditioners cimport *
from .jacobians cimport *
cdef extern from "cantera/oneD/DomainFactory.h" namespace "Cantera":
@ -115,8 +115,8 @@ cdef extern from "cantera/oneD/Sim1D.h":
void restoreSteadySolution() except +translate_exception
void setMaxTimeStepCount(int)
int maxTimeStepCount()
void setLinearSolver(shared_ptr[CxxPreconditionerBase]) except +translate_exception
shared_ptr[CxxPreconditionerBase] linearSolver()
void setLinearSolver(shared_ptr[CxxSystemJacobian]) except +translate_exception
shared_ptr[CxxSystemJacobian] linearSolver()
void getInitialSoln() except +translate_exception
void solve(int, cbool) except +translate_exception
void refine(int) except +translate_exception

View File

@ -1019,10 +1019,10 @@ cdef class Sim1D:
systems as part of each Newton iteration. The default is a banded, direct
solver.
"""
return PreconditionerBase.wrap(self.sim.linearSolver())
return SystemJacobian.wrap(self.sim.linearSolver())
@linear_solver.setter
def linear_solver(self, PreconditionerBase precon):
def linear_solver(self, SystemJacobian precon):
self.sim.setLinearSolver(precon.pbase)
def set_initial_guess(self, *args, **kwargs):

View File

@ -7,16 +7,16 @@
from .ctcxx cimport *
from .kinetics cimport CxxSparseMatrix
cdef extern from "cantera/numerics/PreconditionerBase.h" namespace "Cantera":
cdef cppclass CxxPreconditionerBase "Cantera::PreconditionerBase":
CxxPreconditionerBase()
cdef extern from "cantera/numerics/SystemJacobian.h" namespace "Cantera":
cdef cppclass CxxSystemJacobian "Cantera::SystemJacobian":
CxxSystemJacobian()
string preconditionerSide()
string type()
void setPreconditionerSide(string) except +translate_exception
cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera":
cdef cppclass CxxAdaptivePreconditioner "Cantera::AdaptivePreconditioner" \
(CxxPreconditionerBase):
(CxxSystemJacobian):
CxxAdaptivePreconditioner() except +translate_exception
void setThreshold(double threshold)
double threshold()
@ -28,24 +28,23 @@ cdef extern from "cantera/numerics/AdaptivePreconditioner.h" namespace "Cantera"
CxxSparseMatrix matrix() except +translate_exception
cdef extern from "cantera/oneD/MultiJac.h" namespace "Cantera":
cdef cppclass CxxMultiJac "Cantera::MultiJac" (CxxPreconditionerBase):
cdef cppclass CxxMultiJac "Cantera::MultiJac" (CxxSystemJacobian):
CxxMultiJac() except +translate_exception
cdef extern from "cantera/numerics/PreconditionerFactory.h" namespace "Cantera":
cdef shared_ptr[CxxPreconditionerBase] newPreconditioner(string) except\
cdef extern from "cantera/numerics/SystemJacobianFactory.h" namespace "Cantera":
cdef shared_ptr[CxxSystemJacobian] newSystemJacobian(string) except\
+translate_exception
cdef class PreconditionerBase:
cdef class SystemJacobian:
@staticmethod
cdef wrap(shared_ptr[CxxPreconditionerBase])
cdef wrap(shared_ptr[CxxSystemJacobian])
cdef set_cxx_object(self)
cdef shared_ptr[CxxPreconditionerBase] pbase
cdef shared_ptr[CxxSystemJacobian] pbase
cdef class AdaptivePreconditioner(PreconditionerBase):
cdef class AdaptivePreconditioner(SystemJacobian):
cdef set_cxx_object(self)
cdef CxxAdaptivePreconditioner* preconditioner
cdef class BandedJacobian(PreconditionerBase):
cdef class BandedJacobian(SystemJacobian):
cdef set_cxx_object(self)
cdef CxxMultiJac* preconditioner
cdef CxxMultiJac* jacobian

View File

@ -5,43 +5,43 @@ from ._utils cimport stringify, pystr, py_to_anymap, anymap_to_py
from .kinetics cimport get_from_sparse
# dictionary to store reaction rate classes
cdef dict _preconditioner_class_registry = {}
cdef dict _class_registry = {}
cdef class PreconditionerBase:
cdef class SystemJacobian:
"""
Common base class for preconditioners.
Common base class for Jacobian matrices used in the solution of nonlinear systems.
"""
precon_type = "PreconditionerBase"
precon_linear_solver_type = "GMRES"
_type = "SystemJacobian"
linear_solver_type = "GMRES"
def __cinit__(self, *args, init=True, **kwargs):
if init:
self.pbase = newPreconditioner(stringify(self.precon_type))
self.pbase = newSystemJacobian(stringify(self._type))
@staticmethod
cdef wrap(shared_ptr[CxxPreconditionerBase] base):
cdef wrap(shared_ptr[CxxSystemJacobian] base):
"""
Wrap a C++ PreconditionerBase object with a Python object of the correct
Wrap a C++ SystemJacobian object with a Python object of the correct
derived type.
"""
# Ensure all derived types are registered
if not _preconditioner_class_registry:
if not _class_registry:
def register_subclasses(cls):
for c in cls.__subclasses__():
_preconditioner_class_registry[getattr(c, "precon_type")] = c
_class_registry[getattr(c, "_type")] = c
register_subclasses(c)
# Update global class registry
register_subclasses(PreconditionerBase)
register_subclasses(SystemJacobian)
cxx_type = pystr(base.get().type())
cls = _preconditioner_class_registry.get(cxx_type, PreconditionerBase)
cdef PreconditionerBase precon
precon = cls(init=False)
precon.pbase = base
precon.set_cxx_object()
return precon
cls = _class_registry.get(cxx_type, SystemJacobian)
cdef SystemJacobian jac
jac = cls(init=False)
jac.pbase = base
jac.set_cxx_object()
return jac
cdef set_cxx_object(self):
pass
@ -58,9 +58,9 @@ cdef class PreconditionerBase:
def __set__(self, side):
self.pbase.get().setPreconditionerSide(stringify(side))
cdef class AdaptivePreconditioner(PreconditionerBase):
precon_type = "Adaptive"
precon_linear_solver_type = "GMRES"
cdef class AdaptivePreconditioner(SystemJacobian):
_type = "Adaptive"
linear_solver_type = "GMRES"
def __cinit__(self, *args, **kwargs):
self.preconditioner = <CxxAdaptivePreconditioner*>(self.pbase.get())
@ -139,9 +139,9 @@ cdef class AdaptivePreconditioner(PreconditionerBase):
return get_from_sparse(smat, smat.rows(), smat.cols())
cdef class BandedJacobian(PreconditionerBase):
precon_type = "banded-direct"
precon_linear_solver_type = "direct"
cdef class BandedJacobian(SystemJacobian):
_type = "banded-direct"
linear_solver_type = "direct"
cdef set_cxx_object(self):
self.preconditioner = <CxxMultiJac*>self.pbase.get()
self.jacobian = <CxxMultiJac*>self.pbase.get()

View File

@ -7,7 +7,7 @@
from .ctcxx cimport *
from .kinetics cimport *
from .func1 cimport *
from .preconditioners cimport *
from .jacobians cimport *
cdef extern from "cantera/numerics/Integrator.h" namespace "Cantera":
# SUNDIALS integrator
@ -204,7 +204,7 @@ cdef extern from "cantera/zerodim.h" namespace "Cantera":
string sensitivityParameterName(size_t) except +translate_exception
void setLinearSolverType(string& integratorType) except +translate_exception
string linearSolverType()
void setPreconditioner(shared_ptr[CxxPreconditionerBase] preconditioner)
void setPreconditioner(shared_ptr[CxxSystemJacobian] preconditioner)
void setDerivativeSettings(CxxAnyMap&)
CxxAnyMap solverStats() except +translate_exception

View File

@ -1982,11 +1982,11 @@ cdef class ReactorNet:
property preconditioner:
"""Preconditioner associated with integrator"""
def __set__(self, PreconditionerBase precon):
def __set__(self, SystemJacobian precon):
# set preconditioner
self.net.setPreconditioner(precon.pbase)
# set problem type as default of preconditioner
self.linear_solver_type = precon.precon_linear_solver_type
self.linear_solver_type = precon.linear_solver_type
property linear_solver_type:
"""

View File

@ -1,42 +0,0 @@
//! @file PreconditionerFactory.cpp
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#include "cantera/numerics/PreconditionerFactory.h"
#include "cantera/numerics/AdaptivePreconditioner.h"
#include "cantera/oneD/MultiJac.h"
namespace Cantera
{
PreconditionerFactory* PreconditionerFactory::factory() {
std::unique_lock<std::mutex> lock(precon_mutex);
if (!s_factory) {
s_factory = new PreconditionerFactory;
}
return s_factory;
};
//! Delete preconditioner factory
void PreconditionerFactory::deleteFactory() {
std::unique_lock<std::mutex> lock(precon_mutex);
delete s_factory;
s_factory = 0;
};
PreconditionerFactory* PreconditionerFactory::s_factory = 0;
std::mutex PreconditionerFactory::precon_mutex;
PreconditionerFactory::PreconditionerFactory()
{
reg("Adaptive", []() { return new AdaptivePreconditioner(); });
reg("banded-direct", []() { return new MultiJac(); });
}
shared_ptr<PreconditionerBase> newPreconditioner(const string& precon)
{
return shared_ptr<PreconditionerBase>(PreconditionerFactory::factory()->create(precon));
};
}

View File

@ -0,0 +1,41 @@
//! @file SystemJacobianFactory.cpp
// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.
#include "cantera/numerics/SystemJacobianFactory.h"
#include "cantera/numerics/AdaptivePreconditioner.h"
#include "cantera/oneD/MultiJac.h"
namespace Cantera
{
SystemJacobianFactory* SystemJacobianFactory::factory() {
std::unique_lock<std::mutex> lock(jac_mutex);
if (!s_factory) {
s_factory = new SystemJacobianFactory;
}
return s_factory;
};
void SystemJacobianFactory::deleteFactory() {
std::unique_lock<std::mutex> lock(jac_mutex);
delete s_factory;
s_factory = 0;
};
SystemJacobianFactory* SystemJacobianFactory::s_factory = 0;
std::mutex SystemJacobianFactory::jac_mutex;
SystemJacobianFactory::SystemJacobianFactory()
{
reg("Adaptive", []() { return new AdaptivePreconditioner(); });
reg("banded-direct", []() { return new MultiJac(); });
}
shared_ptr<SystemJacobian> newSystemJacobian(const string& type)
{
return shared_ptr<SystemJacobian>(SystemJacobianFactory::factory()->create(type));
};
}

View File

@ -7,7 +7,7 @@
#include "cantera/numerics/Func1.h"
#include "cantera/oneD/MultiNewton.h"
#include "cantera/base/AnyMap.h"
#include "cantera/numerics/PreconditionerFactory.h"
#include "cantera/numerics/SystemJacobianFactory.h"
#include <fstream>
#include <ctime>
@ -102,7 +102,7 @@ MultiNewton& OneDim::newton()
return *m_newt;
}
void OneDim::setLinearSolver(shared_ptr<PreconditionerBase> solver)
void OneDim::setLinearSolver(shared_ptr<SystemJacobian> solver)
{
m_jac = solver;
m_jac->initialize(size());
@ -219,7 +219,7 @@ void OneDim::resize()
// delete the current Jacobian evaluator and create a new one
if (!m_jac) {
m_jac = newPreconditioner("banded-direct");
m_jac = newSystemJacobian("banded-direct");
}
m_jac->initialize(size());
m_jac->setBandwidth(bandwidth());

View File

@ -155,7 +155,7 @@ void ReactorNet::setLinearSolverType(const string& linSolverType)
m_integrator_init = false;
}
void ReactorNet::setPreconditioner(shared_ptr<PreconditionerBase> preconditioner)
void ReactorNet::setPreconditioner(shared_ptr<SystemJacobian> preconditioner)
{
m_precon = preconditioner;
m_integrator_init = false;

View File

@ -4,7 +4,7 @@
#include "cantera/zerodim.h"
#include "cantera/base/Interface.h"
#include "cantera/numerics/eigen_sparse.h"
#include "cantera/numerics/PreconditionerFactory.h"
#include "cantera/numerics/SystemJacobianFactory.h"
#include "cantera/numerics/AdaptivePreconditioner.h"
using namespace Cantera;
@ -188,7 +188,7 @@ TEST(AdaptivePreconditionerTests, test_precon_solver_stats)
ReactorNet network;
network.addReactor(reactor);
// setup preconditioner
shared_ptr<PreconditionerBase> precon_ptr = newPreconditioner("Adaptive");
shared_ptr<SystemJacobian> precon_ptr = newSystemJacobian("Adaptive");
network.setPreconditioner(precon_ptr);
EXPECT_THROW(network.step(), CanteraError);
// take a step