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

Cleaned up Doxygen documentation for class BEulerInt

Browse source
This commit is contained in:
Ray Speth 2013-05-29 00:11:17 +00:00
parent 0dba3d47a6
commit 7e8276809f
2 changed files with 286 additions and 528 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

356
include/cantera/numerics/BEulerInt.h
View file

@ -16,7 +16,6 @@
#include "cantera/base/utilities.h"
#include "cantera/base/ctexceptions.h"
#include "cantera/numerics/Integrator.h"
#include "cantera/numerics/ResidJacEval.h"
@ -45,6 +44,10 @@ enum BEulerMethodType {
class BEulerErr : public CanteraError
{
public:
/**
* Exception thrown when a BEuler error is encountered. We just call the
* Cantera Error handler in the initialization list.
*/
explicit BEulerErr(const std::string& msg);
};
@ -52,27 +55,38 @@ public:
#define BEULER_JAC_ANAL 2
#define BEULER_JAC_NUM 1
/**
/*!
* Wrapper class for 'beuler' integrator
* We derive the class from the class Integrator
*/
class BEulerInt : public Integrator
{
public:
//! The default constructor doesn't take an argument.
/*!
* Constructor. Default settings: dense jacobian, no user-supplied
* Jacobian function, Newton iteration.
*/
BEulerInt();
//! Destructor
virtual ~BEulerInt();
virtual void setTolerances(double reltol, size_t n, double* abstol);
virtual void setTolerances(double reltol, double abstol);
virtual void setProblemType(int probtype);
//! Find the initial conditions for y and ydot.
virtual void initializeRJE(double t0, ResidJacEval& func);
virtual void reinitializeRJE(double t0, ResidJacEval& func);
virtual double integrateRJE(double tout, double tinit = 0.0);
// This routine advances the calculations one step using a predictor
// corrector approach. We use an implicit algorithm here.
virtual doublereal step(double tout);
//! Set the solution weights. This is a very important routine as it affects
//! quite a few operations involving convergence.
virtual void setSolnWeights();
virtual double& solution(size_t k) {
return m_y_n[k];
}
@ -82,6 +96,8 @@ public:
int nEquations() const {
return m_neq;
}
//! Return the total number of function evaluations
virtual int nEvals() const;
virtual void setMethodBEMT(BEulerMethodType t);
virtual void setIterator(IterType t);
@ -99,68 +115,201 @@ public:
double damp,
int num_entries);
//! This routine controls when the solution is printed
/*!
* @param printSolnStepInterval If greater than 0, then the soln is
* printed every printSolnStepInterval steps.
* @param printSolnNumberToTout The solution is printed at regular
* invervals a total of "printSolnNumberToTout" times.
* @param printSolnFirstSteps The solution is printed out the first
* "printSolnFirstSteps" steps. After these steps the
* other parameters determine the printing. default = 0
* @param dumpJacobians Dump jacobians to disk.
*/
virtual void setPrintSolnOptions(int printSolnStepInterval,
int printSolnNumberToTout,
int printSolnFirstSteps = 0,
bool dumpJacobians = false);
//! Set the options for the nonlinear method
/*!
* Defaults are set in the .h file. These are the defaults:
* min_newt_its = 0
* matrixConditioning = false
* colScaling = false
* rowScaling = true
*/
void setNonLinOptions(int min_newt_its = 0,
bool matrixConditioning = false,
bool colScaling = false,
bool rowScaling = true);
virtual void setPrintFlag(int print_flag);
//! Set the column scaling vector at the current time
virtual void setColumnScales();
/**
* calculate the solution error norm
* Calculate the solution error norm. if printLargest is true, then a table
* of the largest values is printed to standard output.
*/
virtual double soln_error_norm(const double* const,
bool printLargest = false);
virtual void setInitialTimeStep(double delta_t);
void beuler_jac(GeneralMatrix&, double* const,
/*!
* Function called by to evaluate the Jacobian matrix and the current
* residual at the current time step.
* @param J = Jacobian matrix to be filled in
* @param f = Right hand side. This routine returns the current
* value of the rhs (output), so that it does
* not have to be computed again.
*/
void beuler_jac(GeneralMatrix& J, double* const f,
double, double, double* const, double* const, int);
protected:
//! Internal routine that sets up the fixed length storage based on
//! the size of the problem to solve.
void internalMalloc();
/**
* Internal function to calculate the predicted solution
* at a time step.
/*!
* Function to calculate the predicted solution vector, m_y_pred_n for the
* (n+1)th time step. This routine can be used by a first order - forward
* Euler / backward Euler predictor / corrector method or for a second order
* Adams-Bashforth / Trapezoidal Rule predictor / corrector method. See
* Nachos documentation Sand86-1816 and Gresho, Lee, Sani LLNL report UCRL -
* 83282 for more information.
*
* on input:
*
* N - number of unknowns
* order - indicates order of method
* = 1 -> first order forward Euler/backward Euler
* predictor/corrector
* = 2 -> second order Adams-Bashforth/Trapezoidal Rule
* predictor/corrector
*
* delta_t_n - magnitude of time step at time n (i.e., = t_n+1 - t_n)
* delta_t_nm1 - magnitude of time step at time n - 1 (i.e., = t_n - t_n-1)
* y_n[] - solution vector at time n
* y_dot_n[] - acceleration vector from the predictor at time n
* y_dot_nm1[] - acceleration vector from the predictor at time n - 1
*
* on output:
*
* m_y_pred_n[] - predicted solution vector at time n + 1
*/
void calc_y_pred(int);
/**
* Internal function to calculate the time derivative at the
* new step
/*!
* Function to calculate the acceleration vector ydot for the first or
* second order predictor/corrector time integrator. This routine can be
* called by a first order - forward Euler / backward Euler predictor /
* corrector or for a second order Adams - Bashforth / Trapezoidal Rule
* predictor / corrector. See Nachos documentation Sand86-1816 and Gresho,
* Lee, Sani LLNL report UCRL - 83282 for more information.
*
* on input:
*
* N - number of local unknowns on the processor
* This is equal to internal plus border unknowns.
* order - indicates order of method
* = 1 -> first order forward Euler/backward Euler
* predictor/corrector
* = 2 -> second order Adams-Bashforth/Trapezoidal Rule
* predictor/corrector
*
* delta_t_n - Magnitude of the current time step at time n
* (i.e., = t_n - t_n-1)
* y_curr[] - Current Solution vector at time n
* y_nm1[] - Solution vector at time n-1
* ydot_nm1[] - Acceleration vector at time n-1
*
* on output:
*
* ydot_curr[] - Current acceleration vector at time n
*
* Note we use the current attribute to denote the possibility that
* y_curr[] may not be equal to m_y_n[] during the nonlinear solve
* because we may be using a look-ahead scheme.
*/
void calc_ydot(int, double*, double*);
/**
* Internal function to calculate the time step truncation
* error for a predictor corrector time step
/*!
* Calculates the time step truncation error estimate from a very simple
* formula based on Gresho et al. This routine can be called for a first
* order - forward Euler/backward Euler predictor/ corrector and for a
* second order Adams- Bashforth/Trapezoidal Rule predictor/corrector. See
* Nachos documentation Sand86-1816 and Gresho, Lee, LLNL report UCRL -
* 83282 for more information.
*
* on input:
*
* abs_error - Generic absolute error tolerance
* rel_error - Generic realtive error tolerance
* x_coor[] - Solution vector from the implicit corrector
* x_pred_n[] - Solution vector from the explicit predictor
*
* on output:
*
* delta_t_n - Magnitude of next time step at time t_n+1
* delta_t_nm1 - Magnitude of previous time step at time t_n
*/
double time_error_norm();
/**
* Internal function to calculate the time step for the
* next step based on the time-truncation error on the
* current time step
/*!
* Time step control function for the selection of the time step size based on
* a desired accuracy of time integration and on an estimate of the relative
* error of the time integration process. This routine can be called for a
* first order - forward Euler/backward Euler predictor/ corrector and for a
* second order Adams- Bashforth/Trapezoidal Rule predictor/corrector. See
* Nachos documentation Sand86-1816 and Gresho, Lee, Sani LLNL report UCRL -
* 83282 for more information.
*
* on input:
*
* order - indicates order of method
* = 1 -> first order forward Euler/backward Euler
* predictor/corrector
* = 2 -> second order forward Adams-Bashforth/Trapezoidal
* rule predictor/corrector
*
* delta_t_n - Magnitude of time step at time t_n
* delta_t_nm1 - Magnitude of time step at time t_n-1
* rel_error - Generic realtive error tolerance
* time_error_factor - Estimated value of the time step truncation error
* factor. This value is a ratio of the computed
* error norms. The premultiplying constants
* and the power are not yet applied to normalize the
* predictor/corrector ratio. (see output value)
*
* on output:
*
* return - delta_t for the next time step
* If delta_t is negative, then the current time step is
* rejected because the time-step truncation error is
* too large. The return value will contain the negative
* of the recommended next time step.
*
* time_error_factor - This output value is normalized so that
* values greater than one indicate the current time
* integration error is greater than the user
* specified magnitude.
*/
double time_step_control(int m_order, double time_error_factor);
//! Solve a nonlinear system
/*!
*
* Find the solution to F(X, xprime) = 0 by damped Newton iteration. On
* entry, y_comm[] contains an initial estimate of the solution and
* ydot_comm[] contains an estimate of the derivative.
* On successful return, y_comm[] contains the converged solution
* and ydot_comm[] contains the derivative
*
*
* @param y_comm[] Contains the input solution. On output y_comm[] contains
* the converged solution
* @param ydot_comm Contains the input derivative solution. On output y_comm[] contains
* the converged derivative solution
* @param ydot_comm Contains the input derivative solution. On output
* y_comm[] contains the converged derivative solution
* @param CJ Inverse of the time step
* @param time_curr Current value of the time
* @param jac Jacobian
@ -179,8 +328,10 @@ protected:
int loglevel);
/**
* Compute the undamped Newton step. The residual function is
* evaluated at x, but the Jacobian is not recomputed.
* Compute the undamped Newton step. The residual function is
* evaluated at the current time, t_n, at the current values of the
* solution vector, m_y_n, and the solution time derivative, m_ydot_n,
* but the Jacobian is not recomputed.
*/
void doNewtonSolve(double, double*, double*, double*,
GeneralMatrix&, int);
@ -206,10 +357,8 @@ protected:
* couldn't possibly be representative if the
* variable is changed by a lot. (true for
* nonlinear systems, false for linear systems)
* Maximum increase in variable in any one newton iteration:
* factor of 2
* Maximum decrease in variable in any one newton iteration:
* factor of 5
* Maximum increase in variable in any one newton iteration: factor of 2
* Maximum decrease in variable in any one newton iteration: factor of 5
*
* @param y Current value of the solution
* @param step0 Current raw step change in y[]
@ -220,26 +369,25 @@ protected:
*/
double boundStep(const double* const y, const double* const step0, int loglevel);
/*
* Damp step
/*!
* On entry, step0 must contain an undamped Newton step for the
* solution x0. This method attempts to find a damping coefficient
* such that the next undamped step would have a norm smaller than
* that of step0. If successful, the new solution after taking the
* damped step is returned in y1, and the undamped step at y1 is
* returned in step1.
*/
int dampStep(double, const double*, const double*,
const double*, double*, double*,
double*, double&, GeneralMatrix&, int&, bool, int&);
/*
* Compute Residual Weights
*/
//! Compute Residual Weights
void computeResidWts(GeneralMatrix& jac);
/*
* Filter a new step
*/
//! Filter a new step
double filterNewStep(double, double*, double*);
/*
* get the next time to print out
*/
//! Get the next time to print out
double getPrintTime(double time_current);
/********************** Member data ***************************/
@ -284,10 +432,10 @@ protected:
* value of atol. If m_itol = 0, the all atols are equal.
*/
int m_itol;
/**
* Relative time truncation error tolerances
*/
//! Relative time truncation error tolerances
double m_reltol;
/**
* Absolute time truncation error tolerances, when uniform
* for all variables.
@ -298,45 +446,35 @@ protected:
* when not uniform for all variables.
*/
vector_fp m_abstol;
/**
* Error Weights. This is a surprisingly important quantity.
*/
//! Error Weights. This is a surprisingly important quantity.
vector_fp m_ewt;
//! Maximum step size
double m_hmax;
/**
* Maximum integration order
*/
//! Maximum integration order
int m_maxord;
/**
* Current integration order
*/
//! Current integration order
int m_order;
/**
* Time step number
*/
//! Time step number
int m_time_step_num;
int m_time_step_attempts;
/**
* Max time steps allowed
*/
//! Max time steps allowed
int m_max_time_step_attempts;
/**
* Number of initial time steps to take where the
* time truncation error tolerances are not checked. Instead
* the delta T is uniform
* Number of initial time steps to take where the time truncation error
* tolerances are not checked. Instead the delta T is uniform
*/
int m_numInitialConstantDeltaTSteps;
/**
* Failure Counter -> keeps track of the number
* of consequetive failures
*/
//! Failure Counter -> keeps track of the number of consecutive failures
int m_failure_counter;
/**
* Minimum Number of Newton Iterations per nonlinear step
* default = 0
*/
//! Minimum Number of Newton Iterations per nonlinear step. default = 0
int m_min_newt_its;
/************************
* PRINTING OPTIONS
@ -355,25 +493,17 @@ protected:
*/
int m_printSolnNumberToTout;
/**
* Number of initial steps that the solution is
* printed out.
* default = 0
*/
//! Number of initial steps that the solution is printed out. default = 0
int m_printSolnFirstSteps;
/**
* Dump Jacobians to disk
* default false
*/
//! Dump Jacobians to disk. default false
bool m_dumpJacobians;
/*********************
* INTERNAL SOLUTION VALUES
*********************/
/**
* Number of equations in the ode integrator
*/
//! Number of equations in the ode integrator
int m_neq;
vector_fp m_y_n;
vector_fp m_y_nm1;
@ -383,17 +513,13 @@ protected:
/************************
* TIME VARIABLES
************************/
/**
* Initial time at the start of the integration
*/
//! Initial time at the start of the integration
double m_t0;
/**
* Final time
*/
//! Final time
double m_time_final;
/**
*
*/
double time_n;
double time_nm1;
double time_nm2;
@ -401,11 +527,9 @@ protected:
double delta_t_nm1;
double delta_t_nm2;
double delta_t_np1;
/**
* Maximum permissible time step
*/
double delta_t_max;
//! Maximum permissible time step
double delta_t_max;
vector_fp m_resid;
vector_fp m_residWts;
@ -414,11 +538,9 @@ protected:
vector_fp m_rowScales;
vector_fp m_colScales;
/**
* Pointer to the jacobian representing the
* time dependent problem
*/
//! Pointer to the jacobian representing the time dependent problem
GeneralMatrix* tdjac_ptr;
/**
* Determines the level of printing for each time
* step.
@ -432,34 +554,28 @@ protected:
/***************************************************************************
* COUNTERS OF VARIOUS KINDS
***************************************************************************/
/**
* Number of function evaluations
*/
//! Number of function evaluations
int m_nfe;
/**
* Number of Jacobian Evaluations and
* factorization steps (they are the same)
*/
int m_nJacEval;
/**
* Number of total newton iterations
*/
//! Number of total newton iterations
int m_numTotalNewtIts;
/**
* Total number of linear iterations
*/
//! Total number of linear iterations
int m_numTotalLinearSolves;
/**
* Total number of convergence failures.
*/
//! Total number of convergence failures.
int m_numTotalConvFails;
/**
* Total Number of time truncation error failures
*/
//! Total Number of time truncation error failures
int m_numTotalTruncFails;
/*
*
*/
int num_failures;
};

458
src/numerics/BEulerInt.cpp
View file

@ -1,6 +1,5 @@
/**
* @file BEulerInt.cpp
*
*/
/*
@ -19,21 +18,11 @@ using namespace std;
namespace Cantera
{
//================================================================================================
/*
* Exception thrown when a BEuler error is encountered. We just call the
* Cantera Error handler in the initialization list
*/
BEulerErr::BEulerErr(const std::string& msg) :
CanteraError("BEulerInt", msg)
{
}
//================================================================================================
/*
* Constructor. Default settings: dense jacobian, no user-supplied
* Jacobian function, Newton iteration.
*/
BEulerInt::BEulerInt() :
m_iter(Newton_Iter),
m_method(BEulerVarStep),
@ -92,15 +81,12 @@ BEulerInt::BEulerInt() :
{
}
//================================================================================================
/*
* Destructor
*/
BEulerInt::~BEulerInt()
{
delete tdjac_ptr;
}
//================================================================================================
void BEulerInt::setTolerances(double reltol, size_t n, double* abstol)
{
m_itol = 1;
@ -114,64 +100,39 @@ void BEulerInt::setTolerances(double reltol, size_t n, double* abstol)
}
m_reltol = reltol;
}
//================================================================================================
void BEulerInt::setTolerances(double reltol, double abstol)
{
m_itol = 0;
m_reltol = reltol;
m_abstols = abstol;
}
//================================================================================================
void BEulerInt::setProblemType(int jacFormMethod)
{
m_jacFormMethod = jacFormMethod;
}
//================================================================================================
void BEulerInt::setMethodBEMT(BEulerMethodType t)
{
m_method = t;
}
//================================================================================================
void BEulerInt::setMaxStep(doublereal hmax)
{
m_hmax = hmax;
}
//================================================================================================
void BEulerInt::setMaxNumTimeSteps(int maxNumTimeSteps)
{
m_max_time_step_attempts = maxNumTimeSteps;
}
//================================================================================================
void BEulerInt::setNumInitialConstantDeltaTSteps(int num)
{
m_numInitialConstantDeltaTSteps = num;
}
//================================================================================================
/*
*
* setPrintSolnOptins():
*
* This routine controls when the solution is printed
*
* @param printStepInterval If greater than 0, then the
* soln is printed every printStepInterval
* steps.
*
* @param printNumberToTout The solution is printed at
* regular invervals a total of
* "printNumberToTout" times.
*
* @param printSolnFirstSteps The solution is printed out
* the first "printSolnFirstSteps"
* steps. After these steps the other
* parameters determine the printing.
* default = 0
*
* @param dumpJacobians Dump jacobians to disk.
*
* default = false
*
*/
void BEulerInt::setPrintSolnOptions(int printSolnStepInterval,
int printSolnNumberToTout,
int printSolnFirstSteps,
@ -182,24 +143,12 @@ void BEulerInt::setPrintSolnOptions(int printSolnStepInterval,
m_printSolnFirstSteps = printSolnFirstSteps;
m_dumpJacobians = dumpJacobians;
}
//================================================================================================
void BEulerInt::setIterator(IterType t)
{
m_iter = t;
}
//================================================================================================
/*
*
* setNonLinOptions()
*
* Set the options for the nonlinear method
*
* Defaults are set in the .h file. These are the defaults:
* min_newt_its = 0
* matrixConditioning = false
* colScaling = false
* rowScaling = true
*/
void BEulerInt::setNonLinOptions(int min_newt_its, bool matrixConditioning,
bool colScaling, bool rowScaling)
{
@ -214,34 +163,17 @@ void BEulerInt::setNonLinOptions(int min_newt_its, bool matrixConditioning,
m_rowScales.assign(m_neq, 1.0);
}
}
//================================================================================================
/*
*
* setInitialTimeStep():
*
* Set the initial time step. Right now, we set the
* time step by setting delta_t_np1.
*/
void BEulerInt::setInitialTimeStep(double deltaT)
{
delta_t_np1 = deltaT;
}
//================================================================================================
/*
* setPrintFlag():
*
*/
void BEulerInt::setPrintFlag(int print_flag)
{
m_print_flag = print_flag;
}
//================================================================================================
/*
*
* initialize():
*
* Find the initial conditions for y and ydot.
*/
void BEulerInt::initializeRJE(double t0, ResidJacEval& func)
{
m_neq = func.nEquations();
@ -265,12 +197,7 @@ void BEulerInt::initializeRJE(double t0, ResidJacEval& func)
delta_t_n = 0.0;
delta_t_nm1 = 0.0;
}
//================================================================================================
/*
*
* reinitialize():
*
*/
void BEulerInt::reinitializeRJE(double t0, ResidJacEval& func)
{
m_neq = func.nEquations();
@ -291,12 +218,7 @@ void BEulerInt::reinitializeRJE(double t0, ResidJacEval& func)
m_func = &func;
}
//================================================================================================
/*
*
* getPrintTime():
*
*/
double BEulerInt::getPrintTime(double time_current)
{
double tnext;
@ -311,24 +233,12 @@ double BEulerInt::getPrintTime(double time_current)
}
return 1.0E300;
}
//================================================================================================
/*
* nEvals():
*
* Return the total number of function evaluations
*/
int BEulerInt::nEvals() const
{
return m_nfe;
}
//================================================================================================
/*
*
* internalMalloc():
*
* Internal routine that sets up the fixed length storage based on
* the size of the problem to solve.
*/
void BEulerInt::internalMalloc()
{
m_ewt.assign(m_neq, 0.0);
@ -348,15 +258,7 @@ void BEulerInt::internalMalloc()
}
tdjac_ptr = new SquareMatrix(m_neq);
}
//================================================================================================
/*
* setSolnWeights():
*
* Set the solution weights
* This is a very important routine as it affects quite a few
* operations involving convergence.
*
*/
void BEulerInt::setSolnWeights()
{
int i;
@ -378,32 +280,24 @@ void BEulerInt::setSolnWeights()
}
}
}
//================================================================================================
/*
*
* setColumnScales():
*
* Set the column scaling vector at the current time
*/
void BEulerInt::setColumnScales()
{
m_func->calcSolnScales(time_n, &m_y_n[0], &m_y_nm1[0], &m_colScales[0]);
}
//================================================================================================
/*
* computeResidWts():
*
* We compute residual weights here, which we define as the L_0 norm
* of the Jacobian Matrix, weighted by the solution weights.
* This is the proper way to guage the magnitude of residuals. However,
* it does need the evaluation of the jacobian, and the implementation
* below is slow, but doesn't take up much memory.
*
* Here a small weighting indicates that the change in solution is
* very sensitive to that equation.
*/
void BEulerInt::computeResidWts(GeneralMatrix& jac)
{
/*
* We compute residual weights here, which we define as the L_0 norm
* of the Jacobian Matrix, weighted by the solution weights.
* This is the proper way to guage the magnitude of residuals. However,
* it does need the evaluation of the jacobian, and the implementation
* below is slow, but doesn't take up much memory.
*
* Here a small weighting indicates that the change in solution is
* very sensitive to that equation.
*/
int i, j;
double* data = &(*(jac.begin()));
double value;
@ -415,18 +309,12 @@ void BEulerInt::computeResidWts(GeneralMatrix& jac)
}
}
}
//================================================================================================
/*
* filterNewStep():
*
* void BEulerInt::
*
*/
double BEulerInt::filterNewStep(double timeCurrent, double* y_current, double* ydot_current)
{
return 0.0;
}
//==================================================================================================
static void print_line(const char* str, int n)
{
for (int i = 0; i < n; i++) {
@ -434,7 +322,7 @@ static void print_line(const char* str, int n)
}
printf("\n");
}
//==================================================================================================
/*
* Print out for relevant time step information
*/
@ -471,7 +359,7 @@ static void print_time_step1(int order, int n_time_step, double time,
}
printf("\n\tdelta_t_nm1 = %8.5e\n", delta_t_nm1);
}
//================================================================================================
/*
* Print out for relevant time step information
*/
@ -490,7 +378,7 @@ static void print_time_step2(int time_step_num, int order,
print_line("=", 80);
printf("\n");
}
//================================================================================================
/*
* Print Out descriptive information on why the current step failed
*/
@ -523,7 +411,7 @@ static void print_time_fail(bool convFailure, int time_step_num,
printf("\n");
print_line("=", 80);
}
//================================================================================================
/*
* Print out the final results and counters
*/
@ -551,7 +439,7 @@ static void print_final(double time, int step_failed,
printf("\n");
print_line("=", 80);
}
//================================================================================================
/*
* Header info for one line comment about a time step
*/
@ -577,7 +465,7 @@ static void print_lvl1_Header(int nTimes)
printf("\n");
print_line("-", 80);
}
//================================================================================================
/*
* One line entry about time step
* rslt -> 4 letter code
@ -595,9 +483,8 @@ static void print_lvl1_summary(
}
printf("\n");
}
//================================================================================================
/*
* subtractRD():
* This routine subtracts 2 numbers. If the difference is less
* than 1.0E-14 times the magnitude of the smallest number,
* then diff returns an exact zero.
@ -626,18 +513,7 @@ double subtractRD(double a, double b)
}
return diff;
}
//================================================================================================
/*
*
* Function called by BEuler to evaluate the Jacobian matrix and the
* current residual at the current time step.
* @param N = The size of the equation system
* @param J = Jacobian matrix to be filled in
* @param f = Right hand side. This routine returns the current
* value of the rhs (output), so that it does
* not have to be computed again.
*
*/
void BEulerInt::beuler_jac(GeneralMatrix& J, double* const f,
double time_curr, double CJ,
double* const y,
@ -751,36 +627,6 @@ void BEulerInt::beuler_jac(GeneralMatrix& J, double* const f,
}
/*
* Function to calculate the predicted solution vector, m_y_pred_n for the
* (n+1)th time step. This routine can be used by a first order - forward
* Euler / backward Euler predictor / corrector method or for a second order
* Adams-Bashforth / Trapezoidal Rule predictor / corrector method. See Nachos
* documentation Sand86-1816 and Gresho, Lee, Sani LLNL report UCRL - 83282 for
* more information.
*
* variables:
*
* on input:
*
* N - number of unknowns
* order - indicates order of method
* = 1 -> first order forward Euler/backward Euler
* predictor/corrector
* = 2 -> second order Adams-Bashforth/Trapezoidal Rule
* predictor/corrector
*
* delta_t_n - magnitude of time step at time n (i.e., = t_n+1 - t_n)
* delta_t_nm1 - magnitude of time step at time n - 1 (i.e., = t_n - t_n-1)
* y_n[] - solution vector at time n
* y_dot_n[] - acceleration vector from the predictor at time n
* y_dot_nm1[] - acceleration vector from the predictor at time n - 1
*
* on output:
*
* m_y_pred_n[] - predicted solution vector at time n + 1
*/
void BEulerInt::calc_y_pred(int order)
{
int i;
@ -807,42 +653,8 @@ void BEulerInt::calc_y_pred(int order)
*/
m_func->filterSolnPrediction(time_n, &m_y_pred_n[0]);
} /* calc_y_pred */
}
/* Function to calculate the acceleration vector ydot for the first or
* second order predictor/corrector time integrator. This routine can be
* called by a first order - forward Euler / backward Euler predictor /
* corrector or for a second order Adams - Bashforth / Trapezoidal Rule
* predictor / corrector. See Nachos documentation Sand86-1816 and Gresho,
* Lee, Sani LLNL report UCRL - 83282 for more information.
*
* variables:
*
* on input:
*
* N - number of local unknowns on the processor
* This is equal to internal plus border unknowns.
* order - indicates order of method
* = 1 -> first order forward Euler/backward Euler
* predictor/corrector
* = 2 -> second order Adams-Bashforth/Trapezoidal Rule
* predictor/corrector
*
* delta_t_n - Magnitude of the current time step at time n
* (i.e., = t_n - t_n-1)
* y_curr[] - Current Solution vector at time n
* y_nm1[] - Solution vector at time n-1
* ydot_nm1[] - Acceleration vector at time n-1
*
* on output:
*
* ydot_curr[] - Current acceleration vector at time n
*
* Note we use the current attribute to denote the possibility that
* y_curr[] may not be equal to m_y_n[] during the nonlinear solve
* because we may be using a look-ahead scheme.
*/
void BEulerInt::
calc_ydot(int order, double* y_curr, double* ydot_curr)
{
@ -863,32 +675,8 @@ calc_ydot(int order, double* y_curr, double* ydot_curr)
}
return;
}
} /************* END calc_ydot () ****************************************/
}
/* This function calculates the time step truncation error estimate
* from a very simple formula based on Gresho et al. This routine can be
* called for a
* first order - forward Euler/backward Euler predictor/ corrector and
* for a
* second order Adams- Bashforth/Trapezoidal Rule predictor/corrector. See
* Nachos documentation Sand86-1816 and Gresho, Lee, LLNL report
* UCRL - 83282
* for more information.
*
* variables:
*
* on input:
*
* abs_error - Generic absolute error tolerance
* rel_error - Generic realtive error tolerance
* x_coor[] - Solution vector from the implicit corrector
* x_pred_n[] - Solution vector from the explicit predictor
*
* on output:
*
* delta_t_n - Magnitude of next time step at time t_n+1
* delta_t_nm1 - Magnitude of previous time step at time t_n
*/
double BEulerInt::time_error_norm()
{
int i;
@ -943,47 +731,6 @@ double BEulerInt::time_error_norm()
return sqrt(rel_norm / m_neq);
}
/*************************************************************************
* Time step control function for the selection of the time step size based on
* a desired accuracy of time integration and on an estimate of the relative
* error of the time integration process. This routine can be called for a
* first order - forward Euler/backward Euler predictor/ corrector and for a
* second order Adams- Bashforth/Trapezoidal Rule predictor/corrector. See
* Nachos documentation Sand86-1816 and Gresho, Lee, Sani LLNL report UCRL -
* 83282 for more information.
*
* variables:
*
* on input:
*
* order - indicates order of method
* = 1 -> first order forward Euler/backward Euler
* predictor/corrector
* = 2 -> second order forward Adams-Bashforth/Trapezoidal
* rule predictor/corrector
*
* delta_t_n - Magnitude of time step at time t_n
* delta_t_nm1 - Magnitude of time step at time t_n-1
* rel_error - Generic realtive error tolerance
* time_error_factor - Estimated value of the time step truncation error
* factor. This value is a ratio of the computed
* error norms. The premultiplying constants
* and the power are not yet applied to normalize the
* predictor/corrector ratio. (see output value)
*
* on output:
*
* return - delta_t for the next time step
* If delta_t is negative, then the current time step is
* rejected because the time-step truncation error is
* too large. The return value will contain the negative
* of the recommended next time step.
*
* time_error_factor - This output value is normalized so that
* values greater than one indicate the current time
* integration error is greater than the user
* specified magnitude.
*/
double BEulerInt::time_step_control(int order, double time_error_factor)
{
double factor = 0.0, power = 0.0, delta_t;
@ -1019,15 +766,8 @@ double BEulerInt::time_step_control(int order, double time_error_factor)
delta_t = factor * delta_t_n;
}
return delta_t;
} /************ END of time_step_control()********************************/
//================================================================================================
/**************************************************************************
*
* integrate():
*
* defaults are located in the .h file. They are as follows:
* time_init = 0.0
*/
}
double BEulerInt::integrateRJE(double tout, double time_init)
{
double time_current;
@ -1186,14 +926,6 @@ double BEulerInt::integrateRJE(double tout, double time_init)
return time_current;
}
/**************************************************************************
*
* step():
*
* This routine advances the calculations one step using a predictor
* corrector approach. We use an implicit algorithm here.
*
*/
double BEulerInt::step(double t_max)
{
double CJ;
@ -1590,8 +1322,6 @@ double BEulerInt::step(double t_max)
return time_n;
}
//-----------------------------------------------------------
// Constants
//-----------------------------------------------------------
@ -1599,17 +1329,10 @@ double BEulerInt::step(double t_max)
const double DampFactor = 4;
const int NDAMP = 10;
//-----------------------------------------------------------
// MultiNewton methods
//-----------------------------------------------------------
/**
* L2 Norm of a delta in the solution
*
* The second argument has a default of false. However,
* if true, then a table of the largest values is printed
* out to standard output.
*/
double BEulerInt::soln_error_norm(const double* const delta_y,
bool printLargest)
{
@ -1663,15 +1386,7 @@ double BEulerInt::soln_error_norm(const double* const delta_y,
#ifdef DEBUG_HKM_JAC
SquareMatrix jacBack();
#endif
/**************************************************************************
*
* doNewtonSolve():
*
* Compute the undamped Newton step. The residual function is
* evaluated at the current time, t_n, at the current values of the
* solution vector, m_y_n, and the solution time derivative, m_ydot_n,
* but the Jacobian is not recomputed.
*/
void BEulerInt::doNewtonSolve(double time_curr, double* y_curr,
double* ydot_curr, double* delta_y,
GeneralMatrix& jac, int loglevel)
@ -1838,39 +1553,6 @@ void BEulerInt::doNewtonSolve(double time_curr, double* y_curr,
m_numTotalLinearSolves++;
}
//================================================================================================
// Bound the Newton step while relaxing the solution
/*
* Return the factor by which the undamped Newton step 'step0'
* must be multiplied in order to keep all solution components in
* all domains between their specified lower and upper bounds.
* Other bounds may be applied here as well.
*
* Currently the bounds are hard coded into this routine:
*
* Minimum value for all variables: - 0.01 * m_ewt[i]
* Maximum value = none.
*
* Thus, this means that all solution components are expected
* to be numerical greater than zero in the limit of time step
* truncation errors going to zero.
*
* Delta bounds: The idea behind these is that the Jacobian
* couldn't possibly be representative if the
* variable is changed by a lot. (true for
* nonlinear systems, false for linear systems)
* Maximum increase in variable in any one newton iteration:
* factor of 2
* Maximum decrease in variable in any one newton iteration:
* factor of 5
*
* @param y Current value of the solution
* @param step0 Current raw step change in y[]
* @param loglevel Log level. This routine produces output if loglevel
* is greater than one
*
* @return Returns the damping coefficient
*/
double BEulerInt::boundStep(const double* const y,
const double* const step0, int loglevel)
{
@ -1937,18 +1619,7 @@ double BEulerInt::boundStep(const double* const y,
}
return fbound;
}
//================================================================================================
/**************************************************************************
*
* dampStep():
*
* On entry, step0 must contain an undamped Newton step for the
* solution x0. This method attempts to find a damping coefficient
* such that the next undamped step would have a norm smaller than
* that of step0. If successful, the new solution after taking the
* damped step is returned in y1, and the undamped step at y1 is
* returned in step1.
*/
int BEulerInt::dampStep(double time_curr, const double* y0,
const double* ydot0, const double* step0,
double* y1, double* ydot1, double* step1,
@ -1956,8 +1627,6 @@ int BEulerInt::dampStep(double time_curr, const double* y0,
int& loglevel, bool writetitle,
int& num_backtracks)
{
// Compute the weighted norm of the undamped step size step0
double s0 = soln_error_norm(step0);
@ -2104,28 +1773,7 @@ int BEulerInt::dampStep(double time_curr, const double* y0,
return -2;
}
}
//================================================================================================
// Solve a nonlinear system
/*
* Find the solution to F(X, xprime) = 0 by damped Newton iteration. On
* entry, y_comm[] contains an initial estimate of the solution and
* ydot_comm[] contains an estimate of the derivative.
* On successful return, y_comm[] contains the converged solution
* and ydot_comm[] contains the derivative
*
*
* @param y_comm[] Contains the input solution. On output y_comm[] contains
* the converged solution
* @param ydot_comm Contains the input derivative solution. On output y_comm[] contains
* the converged derivative solution
* @param CJ Inverse of the time step
* @param time_curr Current value of the time
* @param jac Jacobian
* @param num_newt_its number of newton iterations
* @param num_linear_solves number of linear solves
* @param num_backtracks number of backtracs
* @param loglevel Log level
*/
int BEulerInt::solve_nonlinear_problem(double* const y_comm,
double* const ydot_comm, double CJ,
double time_curr,
@ -2281,11 +1929,7 @@ done:
}
return m;
}
//================================================================================================
/*
*
*
*/
void BEulerInt::
print_solnDelta_norm_contrib(const double* const solnDelta0,
const char* const s0,
@ -2341,7 +1985,5 @@ print_solnDelta_norm_contrib(const double* const solnDelta0,
printf("\t\t ");
print_line("-", 90);
}
//===============================================================================================
} // End of namespace Cantera