research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Changed the interface for ResidJacEval.

Browse source 
Updated NonlinearSolver with a lot of heuristic algorithm changes.
This commit is contained in:
Harry Moffat 2010-10-08 17:16:25 +00:00
parent 7905ea977f
commit 1d0ba2d3b8
5 changed files with 1067 additions and 424 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
Cantera/src/numerics/BEulerInt.cpp
View file

@ -272,7 +272,7 @@ namespace Cantera {
/*
* Get the initial conditions.
*/
func.getInitialConditionsDot(m_t0, m_neq, m_y_n, m_ydot_n);
func.getInitialConditions(m_t0, m_y_n, m_ydot_n);
// Store a pointer to the residual routine in the object
m_func = &func;
@ -672,7 +672,7 @@ namespace Cantera {
* current conditions.
*/
m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f);
m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f, JacBase_ResidEval);
m_nfe++;
m_nJacEval++;
@ -732,7 +732,7 @@ namespace Cantera {
m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, m_wksp,
true, j, dy);
JacDelta_ResidEval, j, dy);
m_nfe++;
double diff;
for (i = 0; i < m_neq; i++) {
@ -1669,7 +1669,7 @@ namespace Cantera {
int irow, jcol;
m_func->evalResidNJ(time_curr, delta_t_n, y_curr,
ydot_curr, delta_y);
ydot_curr, delta_y, Base_ResidEval);
m_nfe++;
int sz = m_func->nEquations();
for (int n = 0; n < sz; n++) {

794
Cantera/src/numerics/NonlinearSolver.cpp
View file

File diff suppressed because it is too large Load diff

164
Cantera/src/numerics/NonlinearSolver.h
View file

@ -87,12 +87,18 @@ namespace Cantera {
* calculate the norm of the solution vector. This will
* involve the column scaling of the matrix
*
* The second argument has a default of false. However,
* The third argument has a default of false. However,
* if true, then a table of the largest values is printed
* out to standard output.
*
* @param delta_y Vector to take the norm of
* @param title Optional title to be printed out
* @param printLargest int indicating how many specific lines should be printed out
* @param dampFactor Current value of the damping factor. Defaults to 1.
* only used for printout out a table.
*/
double solnErrorNorm(const double * const delta_y,
bool printLargest = false);
double solnErrorNorm(const double * const delta_y, const char * title = 0, int printLargest = 0,
const double dampFactor = 1.0);
//! L2 norm of the residual of the equation system
/*!
@ -102,33 +108,42 @@ namespace Cantera {
* The second argument has a default of false. However,
* if true, then a table of the largest values is printed
* out to standard output.
*
* @param resid Vector of the residuals
* @param title Optional title to be printed out
* @param printLargest Number of specific entries to be printed
* @param y Current value of y - only used for printouts
*/
double residErrorNorm(const double * const resid,
bool printLargest = false);
double residErrorNorm(const double * const resid, const char * title = 0, const int printLargest = 0,
const double * const y = 0);
//! Compute the current Residual
/*!
* Compute the time dependent residual of
* the set of equations.
*/
void doTDResidualCalc(const double time_curr, const int typeCalc,
const double * const y_curr,
const double * const ydot_curr, double* const residual,
int loglevel);
// void doTDResidualCalc(const double time_curr, const int typeCalc,
// const double * const y_curr, const double * const ydot_curr, int loglevel);
//! Compute the current Residual
/*!
* Compute the steady state residual of
* the set of equations.
*/
void doSteadyResidualCalc(const double time_curr, const int typeCalc,
const double * const y_curr,
double* const residual, int loglevel);
// void doSteadyResidualCalc(const double time_curr, const int typeCalc,
// const double * const y_curr, int loglevel);
void doResidualCalc(const double time_curr, const int typeCalc,
const double * const y_curr,
const double * const ydot_curr, double* const residual,
int loglevel);
//! Compute the current residual
/*!
* The current value of the residual is storred in the internal work array m_resid.
*
* @param time_curr Value of the time
* @param typeCalc Type of the calculation
* @param y_curr Current value of the solution vector
* @param ydot_curr Current value of the time derivative of the solution vector
*/
void doResidualCalc(const double time_curr, const int typeCalc, const double * const y_curr,
const double * const ydot_curr);
//! Compute the undamped Newton step
/*!
@ -316,6 +331,12 @@ namespace Cantera {
//! Set the column scales
void setColumnScales();
//! Scale the matrix
/*!
*
*/
void scaleMatrix(SquareMatrix& jac, double* y_comm, double* ydot_comm, double time_curr);
//! Print solution norm contribution
void
@ -329,6 +350,41 @@ namespace Cantera {
double damp,
int num_entries);
//! Compute the Residual Weights
/*!
* The residual weights are defined here to be equal to the inverse of the row scaling factors used to
* row scale the matrix, after column scaling is used. They are multiplied by 10-3 because the column
* weights are also multiplied by that same quantity.
*
* The basic idea is that a change in the solution vector on the order of the convergence tolerance
* multiplied by [RJC] which is of order one after row scaling should give you the relative weight
* of the row. Values of the residual for that row can then be normalized by the value of this weight.
* When the tolerance in delta x is achieved, the tolerance in the residual is also achieved.
*/
void computeResidWts();
//! Return the residual weights
/*!
* @param residWts Vector of length neq_
*/
void getResidWts(double * const residWts) const;
//! Check to see if the nonlinear problem has converged
/*!
*
* @return integer is returned. If positive, then the problem has converged
* 1 Successful step was taken: Next step's norm is less than 1.0.
* The final residual norm is less than 1.0.
* 2 Successful step: Next step's norm is less than 0.8.
* This step's norm is less than 1.0.
* The residual norm can be anything.
* 3 Success: The final residual is less than 1.0
* The predicted deltaSoln is below 1.0.
* 0 Not converged yet
*/
int convergenceCheck(int dampCode, double s1);
private:
//! Pointer to the residual and jacobian evaluator for the
@ -347,9 +403,13 @@ namespace Cantera {
//! Soln error weights
std::vector<doublereal> m_ewt;
//! Boolean indicating whether a manual delta bounds has been input.
//! Boolean indicating whether a manual delta bounds has been input.
int m_manualDeltaBoundsSet;
//! Soln Delta bounds magnitudes
std::vector<doublereal> m_deltaBoundsMagnitudes;
@ -359,28 +419,79 @@ namespace Cantera {
std::vector<doublereal> ydot_new;
//! Vector of column scaling factors
std::vector<doublereal> m_colScales;
//! Weights for normalizing the values of the residuals
/*!
* These are computed if row scaling, m_rowScaling, is turned on. They are calculated currently as the
* sum of the absolute values of the rows of the jacobian.
*/
std::vector<doublereal> m_rowScales;
//! Value of the residual for the nonlinear problem
std::vector<doublereal> m_resid;
//! Workspace of length neq_
std::vector<doublereal> m_wksp;
/*****************************************************************************************
* INTERNAL WEIGHTS FOR TAKING SOLUTION NORMS
******************************************************************************************/
//! Vector of residual weights
/*!
* These are used to establish useful and informative weighted norms of the residual vector.
*/
std::vector<doublereal> m_residWts;
//! Norm of the residual at the start of each nonlinear iteration
double m_normResid0;
//! Norm of the residual before damping
double m_normResidFRaw;
//! Norm of the solution update created by the iteration in its raw, undamped form.
double m_normSolnFRaw;
//! Norm of the residual for a trial calculation which may or may not be used
double m_normResidTrial;
//! Vector of the norm
double m_normResidPoints[15];
bool m_resid_scaled;
/*****************************************************************************************
* INTERNAL BOUNDARY INFO FOR SOLUTIONS
*****************************************************************************************/
//! Bounds vector for each species
std::vector<doublereal> m_y_high_bounds;
//! Lower bounds vector for each species
std::vector<doublereal> m_y_low_bounds;
//! Damping factor imposed by hard bounds and by delta bounds
double m_dampBound;
//! Additional damping factor due to bounds on the residual and solution norms
double m_dampRes;
//! Delta t for the current step
double delta_t_n;
//! Counter for the total number of function evaluations
int m_nfe;
//! The type of column scaled used in the solution of the problem
/*!
* If true then colScaling = m_ewt[]
* if false then colScaling = 1.0
* Currently, this is not part of the interface
*/
bool m_colScaling;
//! int indicating whether row scaling is turned on (1) or not (0)
@ -417,11 +528,26 @@ namespace Cantera {
doublereal rtol_;
//! Base value of the absolute tolerance
doublereal atolBase_;
double * m_ydot_nm1;
std::vector<doublereal> atolk_;
//! Determines the level of printing for each time step.
/*!
* 0 -> absolutely nothing is printed for a single time step.
* 1 -> One line summary per solve_nonlinear call
* 2 -> short description, points of interest: Table of nonlinear solve - one line per iteration
* 3 -> Table is included -> More printing per nonlinear iteration (default) that occurs during the table
* 4 -> Summaries of the nonlinear solve iteration as they are occurring -> table no longer printed
* 5 -> Algorithm information on the nonlinear iterates are printed out
* 6 -> Additional info on the nonlinear iterates are printed out
* 7 -> Additional info on the linear solve is printed out.
* 8 -> Info on a per iterate of the linear solve is printed out.
*/
int m_print_flag;
};
}

281
Cantera/src/numerics/ResidJacEval.cpp
View file

@ -23,21 +23,14 @@
using namespace std;
namespace Cantera {
/*************************************************************************
*
* ResidJacEval():
*
* Default constructor for the ResidJacEval class.
*
* atol has a default of 1.0E-13.
*/
//====================================================================================================================
ResidJacEval::ResidJacEval(doublereal atol) :
ResidEval(),
m_atol(atol)
{
}
//====================================================================================================================
// Copy Constructor for the %ResidJacEval object
/*
*/
@ -46,14 +39,11 @@ namespace Cantera {
{
*this = operator=(right);
}
/*
*
*/
//====================================================================================================================
ResidJacEval::~ResidJacEval()
{
}
//====================================================================================================================
ResidJacEval& ResidJacEval::operator=(const ResidJacEval &right) {
if (this == &right) {
return *this;
@ -66,9 +56,8 @@ namespace Cantera {
return *this;
}
// Duplication routine for objects which inherit from
// %ResidJacEval
//====================================================================================================================
// Duplication routine for objects which inherit from %ResidJacEval
/*
* This virtual routine can be used to duplicate %ResidJacEval objects
* inherited from %ResidJacEval even if the application only has
@ -81,16 +70,14 @@ namespace Cantera {
ResidJacEval *ff = new ResidJacEval(*this);
return ff;
}
//====================================================================================================================
int ResidJacEval::nEquations() const {
return neq_;
}
//====================================================================================================================
// Set a global value of the absolute tolerance
/*
*
* setAtol():
*
* Set the absolute tolerance value
* @param atol Value of atol
*/
void ResidJacEval::setAtol(doublereal atol)
{
@ -100,17 +87,17 @@ namespace Cantera {
"atol must be greater than zero");
}
}
/**************************************************************************
//====================================================================================================================
//! Fill in the initial conditions
/*!
* Values for both the solution and the value of ydot may be provided.
*
*
*
* Fill the solution vector with the initial conditions
* at initial time t0.
* @param t0 Time (input)
* @param y Solution vector (output)
* @param ydot Rate of change of solution vector. (output)
*/
void ResidJacEval::
getInitialConditionsDot(const doublereal t0, const size_t leny,
doublereal * const y, doublereal * const ydot) {
getInitialConditions(doublereal t0, doublereal * const y, doublereal * const ydot) {
for (int i = 0; i < neq_; i++) {
y[i] = 0.0;
}
@ -120,65 +107,75 @@ namespace Cantera {
}
}
}
/**************************************************************************
//====================================================================================================================
// This function may be used to create output at various points in the execution of an application.
/*
*
* @param ifunc identity of the call
* 0 Initial call
* 1 Called at the end of every successful time step
* -1 Called at the end of every unsuccessful time step
* 2 Called at the end of every call to integrateRJE()
*
*
* Fill the solution vector with the initial conditions
* at initial time t0.
*
*/
void ResidJacEval::
getInitialConditions(doublereal t0,
doublereal * const y, doublereal * const ydot) {
size_t leny = neq_;
getInitialConditionsDot(t0, leny, y, 0);
}
/**************************************************************************
*
* user_out():
*
* This function may be used to create output at various points in the
* execution of an application.
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input)
*/
void ResidJacEval::
user_out2(const int ifunc, const doublereal t, const doublereal deltaT,
const doublereal *y, const doublereal *ydot) {
}
//====================================================================================================================
// This function may be used to create output at various points in the execution of an application.
/*
* This routine calls user_out2().
*
* @param ifunc identity of the call
* @param t Time (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input)
*/
void ResidJacEval::
user_out(const int ifunc, const doublereal t,
const doublereal *y, const doublereal *ydot) {
user_out2(ifunc, t, 0.0, y, ydot);
}
/**************************************************************************
*
//====================================================================================================================
//! Evaluate the time tracking equations, if any
/*!
* Evaluate time integrated quantities that are calculated at the
* end of every successful time step. This call is made once at the end of every successful
* time step that advances the time. It's also made once at the start of the time stepping.
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
*/
void ResidJacEval::
evalTimeTrackingEqns(const doublereal t, const doublereal deltaT,
const doublereal *y,
const doublereal *ydot) {
evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal *y,
const doublereal *ydot)
{
}
/********************************************************************
//====================================================================================================================
// Return a vector of delta y's for calculation of the numerical Jacobian
/*
* There is a default algorithm provided.
*
* delta_y[i] = atol[i] + 1.0E-6 ysoln[i]
* delta_y[i] = atol[i] + MAX(1.0E-6 ysoln[i] * 0.01 * solnWeights[i])
*
*
* Return a vector of delta y's for calculation of the
* numerical Jacobian
* @param t Time (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param delta_y Value of the delta to be used in calculating the numerical jacobian
* @param solnWeights Value of the solution weights that are used in determining convergence (default = 0)
*/
void ResidJacEval::
calcDeltaSolnVariables(const doublereal t,
const doublereal * const ySoln,
const doublereal * const ySolnDot,
doublereal * const deltaYSoln,
calcDeltaSolnVariables(const doublereal t, const doublereal * const ySoln,
const doublereal * const ySolnDot, doublereal * const deltaYSoln,
const doublereal *const solnWeights)
{
if (!solnWeights) {
@ -192,95 +189,123 @@ namespace Cantera {
}
}
}
/******************************************************************
//====================================================================================================================
// Returns a vector of column scale factors that can be used to column scale Jacobians.
/*
* Default to yScales[] = 1.0
*
* calcSolnScales():
*
* Returns a vector of ysolnScales[] that can be used to column scale
* Jacobians.
* @param t Time (input)
* @param y Solution vector (input, do not modify)
* @param y_old Old Solution vector (input, do not modify)
* @param yScales Value of the column scales
*/
void ResidJacEval::
calcSolnScales(const doublereal t,
const doublereal * const ysoln,
const doublereal * const ysolnOld,
calcSolnScales(const doublereal t, const doublereal * const ysoln, const doublereal * const ysolnOld,
doublereal * const ysolnScales)
{
for (int i = 0; i < neq_; i++) {
ysolnScales[i] = 1.0;
if (ysolnScales[0] == 0.0) {
for (int i = 0; i < neq_; i++) {
ysolnScales[i] = 1.0;
}
}
}
void ResidJacEval::filterSolnPrediction(doublereal t,
doublereal * const y) {
//====================================================================================================================
// Filter the solution predictions
/*
* Codes might provide a predicted solution vector. This routine filters the predicted
* solution vector.
*
* @param t Time (input)
* @param y Solution vector (input, output)
*/
void ResidJacEval::filterSolnPrediction(doublereal t, doublereal * const y)
{
}
/**************************************************************************
//====================================================================================================================
// Evalulate any stopping criteria other than a final time limit
/*
* If we are to stop the time integration for any reason other than reaching a final time limit, tout,
* provide a test here. This call is made at the end of every succesful time step iteration
*
* evalStoppingCriteria()
*
* If there is a stopping critera other than time set it here.
* @return If true, the the time stepping is stopped. If false, then time stepping is stopped if t >= tout
* Defaults to false.
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
*/
bool ResidJacEval::
evalStoppingCritera(doublereal &time_current,
doublereal &delta_t_n,
doublereal *y_n,
doublereal *ydot_n)
evalStoppingCritera(const doublereal t,
const doublereal delta_t,
const doublereal * const y,
const doublereal * const ydot)
{
return false;
}
/**************************************************************************
//====================================================================================================================
// Multiply the matrix by another matrix that leads to better conditioning
/*
* Provide a left sided matrix that will multiply the current jacobian, after scaling
* and lead to a better conditioned system.
* This routine is called just before the matrix is factored.
*
* Original Problem:
* J delta_x = - Resid
*
* matrixConditioning()
* New problem:
* M (J delta_x) = - M Resid
*
* Multiply the matrix by the inverse of a matrix which lead to a
* better conditioned system. The default, specified here, is to
* do nothing.
* @param matrix Pointer to the current jacobian (if zero, it's already been factored)
* @param nrows offsets for the matrix
* @param rhs residual vector. This also needs to be lhs multiplied by M
*/
void ResidJacEval::
matrixConditioning(doublereal * const matrix, const int nrows,
doublereal * const rhs)
matrixConditioning(doublereal * const matrix, const int nrows, doublereal * const rhs)
{
}
/**************************************************************************
*
*/
//====================================================================================================================
// Evaluate the residual function
/*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param resid Value of the residual that is computed (output)
* @param evalType Type of the residual being computed (defaults to Base_ResidEval)
* @param id_x Index of the variable that is being numerically differenced to find
* the jacobian (defaults to -1, which indicates that no variable is being
* differenced or that the residual doesn't take this issue into account)
* @param delta_x Value of the delta used in the numerical differencing
*/
void ResidJacEval::
evalResidNJ(doublereal t, const doublereal deltaT,
const doublereal * y,
const doublereal * ydot,
doublereal * resid,
bool NJevaluation,
int id_x,
doublereal delta_x)
{
printf("Not implemented\n");
std::exit(-1);
evalResidNJ(const doublereal t, const doublereal deltaT, const doublereal * y,
const doublereal * ydot, doublereal * const resid, const ResidEval_Type_Enum evalType,
const int id_x, const doublereal delta_x)
{
throw CanteraError("ResidJacEval::evalResidNJ()", "Not implemented\n");
}
/**************************************************************************
//====================================================================================================================
// Calculate an analytical jacobian and the residual at the current time and values.
/*
* Only called if the jacFormation method is set to analytical
*
* evalJacobian()
*
* Calculate the jacobian and the residual at the current
* time and values.
* Backwards Euler is assumed.
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param J Reference to the SquareMatrix object to be calculated (output)
* @param resid Value of the residual that is computed (output)
*/
void ResidJacEval::
evalJacobian(const doublereal t, const doublereal deltaT,
evalJacobian(const doublereal t, const doublereal delta_t,
const doublereal * const y,
const doublereal * const ydot,
SquareMatrix &J,
doublereal * const resid)
{
printf("Not implemented\n");
std::exit(-1);
throw CanteraError("ResidJacEval::evalJacobian()", "Not implemented\n");
}
//====================================================================================================================
}

244
Cantera/src/numerics/ResidJacEval.h
View file

@ -25,16 +25,38 @@
namespace Cantera {
/**
* A class for full (non-sparse) matrices with Fortran-compatible
* data storage. Adds matrix operations to class Array2D.
//! Differentiates the type of residual evaluations according to functionality
enum ResidEval_Type_Enum
{
//! Base residual calculation for the time-stepping function
Base_ResidEval = 0,
//! Base residual calculation for the Jacobian calculation
JacBase_ResidEval,
//! Delta residual calculation for the Jacbobian calculation
JacDelta_ResidEval,
//! Base residual calculation for the showSolution routine
/*!
* We calculate this when we want to display a solution
*/
Base_ShowSolution
};
//! Wrappers for the function evaluators for Nonlinear solvers and Time steppers
/*!
* A class for full (non-sparse dense matrices with Fortran-compatible data storage.
* The class adds support for identifying what types of calls are made to the residual
* evaluator by adding the ResidEval_Type_Enum class.
*
*
*/
class ResidJacEval : public ResidEval {
public:
/**
* Default constructor
//!Default constructor
/*!
* @param atol Initial value of the global tolerance (defaults to 1.0E-13)
*/
ResidJacEval(doublereal atol = 1.0e-13);
@ -71,106 +93,196 @@ namespace Cantera {
//! Return the number of equations in the equation system
virtual int nEquations() const;
/**
* Evaluate the residual function.
* @param t time (input, do not modify)
* @param y solution vector (input, do not modify)
* @param ydot rate of change of solution vector. (input, do
* not modify)
//! Evaluate the residual function
/*!
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param resid Value of the residual that is computed (output)
* @param evalType Type of the residual being computed (defaults to Base_ResidEval)
* @param id_x Index of the variable that is being numerically differenced to find
* the jacobian (defaults to -1, which indicates that no variable is being
* differenced or that the residual doesn't take this issue into account)
* @param delta_x Value of the delta used in the numerical differencing
*/
virtual void evalResidNJ(doublereal t, const doublereal deltaT,
virtual void evalResidNJ(const doublereal t, const doublereal delta_t,
const doublereal * const y,
const doublereal * const ydot,
doublereal * const resid,
bool NJevaluation = false,
int id_x = 0,
doublereal delta_x = 0.0);
const ResidEval_Type_Enum evalType = Base_ResidEval,
const int id_x = -1,
const doublereal delta_x = 0.0);
/**
* Fill the solution vector with the initial conditions
* at initial time t0.
//! Fill in the initial conditions
/*!
* Values for both the solution and the value of ydot may be provided.
*
* @param t0 Time (input)
* @param y Solution vector (output)
* @param ydot Rate of change of solution vector. (output)
*/
virtual void getInitialConditionsDot(const doublereal t0, size_t leny,
doublereal * const y,
doublereal * const ydot);
virtual void getInitialConditions(const doublereal t0,
doublereal * const y,
doublereal * const ydot);
virtual void getInitialConditions(const doublereal t0, doublereal * const y, doublereal * const ydot);
virtual void filterSolnPrediction(doublereal t,
doublereal * const y);
//! Filter the solution predictions
/*!
* Codes might provide a predicted solution vector. This routine filters the predicted
* solution vector.
*
* @param t Time (input)
* @param y Solution vector (input, output)
*/
virtual void filterSolnPrediction(const doublereal t, doublereal * const y);
//! Set a global value of the absolute tolerance
/*!
* @param atol Value of atol
*/
void setAtol(doublereal atol);
virtual void evalTimeTrackingEqns(const doublereal t, const doublereal deltaT,
const doublereal * const y,
//! Evaluate the time tracking equations, if any
/*!
* Evaluate time integrated quantities that are calculated at the
* end of every successful time step. This call is made once at the end of every successful
* time step that advances the time. It's also made once at the start of the time stepping.
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
*/
virtual void evalTimeTrackingEqns(const doublereal t, const doublereal delta_t, const doublereal * const y,
const doublereal * const ydot);
virtual bool evalStoppingCritera(doublereal &time_current,
doublereal &delta_t_n,
doublereal *y_n,
doublereal *ydot_n);
/**
* Return a vector of delta y's for calculation of the
* numerical Jacobian
//! Evalulate any stopping criteria other than a final time limit
/*!
* If we are to stop the time integration for any reason other than reaching a final time limit, tout,
* provide a test here. This call is made at the end of every succesful time step iteration
*
* @return If true, the the time stepping is stopped. If false, then time stepping is stopped if t >= tout
* Defaults to false.
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
*/
virtual bool evalStoppingCritera(const doublereal t,
const doublereal delta_t,
const doublereal * const y,
const doublereal * const ydot);
//! Return a vector of delta y's for calculation of the numerical Jacobian
/*!
* There is a default algorithm provided.
*
* delta_y[i] = atol[i] + 1.0E-6 ysoln[i]
* delta_y[i] = atol[i] + MAX(1.0E-6 ysoln[i] * 0.01 * solnWeights[i])
*
* @param t Time (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param delta_y Value of the delta to be used in calculating the numerical jacobian
* @param solnWeights Value of the solution weights that are used in determining convergence (default = 0)
*/
virtual void
calcDeltaSolnVariables(const doublereal t,
const doublereal * const ysoln,
const doublereal * const ysolnDot,
doublereal * const deltaYsoln,
const doublereal * const solnWeights=0);
calcDeltaSolnVariables(const doublereal t,
const doublereal * const y,
const doublereal * const ydot,
doublereal * const delta_y,
const doublereal * const solnWeights = 0);
/**
* Returns a vector of ysolnScales[] that can be used to column
* scale Jacobians.
*/
virtual void calcSolnScales(const doublereal t,
const doublereal * const ysoln,
const doublereal * const ysolnOld,
doublereal * const ysolnScales);
/**
* This function may be used to create output at various points in the
* execution of an application.
//! Returns a vector of column scale factors that can be used to column scale Jacobians.
/*!
* Default to yScales[] = 1.0
*
* @param t Time (input)
* @param y Solution vector (input, do not modify)
* @param y_old Old Solution vector (input, do not modify)
* @param yScales Value of the column scales
*/
virtual void calcSolnScales(const doublereal t, const doublereal * const y,
const doublereal * const y_old, doublereal * const yScales);
//! This function may be used to create output at various points in the execution of an application.
/*!
*
* @param ifunc identity of the call
* 0 Initial call
* 1 Called at the end of every successful time step
* -1 Called at the end of every unsuccessful time step
* 2 Called at the end of every call to integrateRJE()
*
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input)
*/
virtual void user_out2(const int ifunc, const doublereal t,
const doublereal deltaT,
const doublereal delta_t,
const doublereal * const y,
const doublereal * const ydot);
//! This function may be used to create output at various points in the execution of an application.
/*!
* This routine calls user_out2().
*
* @param ifunc identity of the call
* @param t Time (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input)
*/
virtual void user_out(const int ifunc, const doublereal t,
const doublereal *y,
const doublereal *ydot);
//! Multiply the matrix by another matrix that leads to better conditioning
/*!
* Provide a left sided matrix that will multiply the current jacobian, after scaling
* and lead to a better conditioned system.
* This routine is called just before the matrix is factored.
*
* Original Problem:
* J delta_x = - Resid
*
* New problem:
* M (J delta_x) = - M Resid
*
* @param matrix Pointer to the current jacobian (if zero, it's already been factored)
* @param nrows offsets for the matrix
* @param rhs residual vector. This also needs to be lhs multiplied by M
*/
virtual void matrixConditioning(doublereal * const matrix, const int nrows,
doublereal * const rhs);
/*********************************************************************
//! Calculate an analytical jacobian and the residual at the current time and values.
/*!
* Only called if the jacFormation method is set to analytical
*
* evalJacobian()
*
* Calculate the jacobian and the residual at the current
* time and values.
* Backwards Euler is assumed.
* @param t Time (input)
* @param delta_t The current value of the time step (input)
* @param y Solution vector (input, do not modify)
* @param ydot Rate of change of solution vector. (input, do not modify)
* @param J Reference to the SquareMatrix object to be calculated (output)
* @param resid Value of the residual that is computed (output)
*/
virtual void evalJacobian(const doublereal t, const doublereal deltaT,
const double* const y,
const double* const ydot,
virtual void evalJacobian(const doublereal t, const doublereal delta_t,
const doublereal* const y,
const doublereal* const ydot,
SquareMatrix &J,
doublereal * const resid);
protected:
//! constant value of atol
doublereal m_atol;
//! Number of equations
int neq_;
};
}