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

Initial commits of BEulerInt. This is a backwards euler stepper

Browse source 
with lots of hooks. Using this to develop a subgrid electrode object.
This commit is contained in:
Harry Moffat 2010-05-23 15:11:55 +00:00
parent 4856da5083
commit 3167a96541
4 changed files with 2853 additions and 2 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

2359
Cantera/src/numerics/BEulerInt.cpp Normal file
View file

File diff suppressed because it is too large Load diff

465
Cantera/src/numerics/BEulerInt.h Normal file
View file

@ -0,0 +1,465 @@
/**
* @file BEulerInt.h
*/
/* $Author: hkmoffa $
* $Date: 2009/01/27 16:50:53 $
* $Revision: 1.19 $
*/
/*
* Copywrite 2004 Sandia Corporation. Under the terms of Contract
* DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
* retains certain rights in this software.
* See file License.txt for licensing information.
*/
#ifndef CT_BEULERINT_H
#define CT_BEULERINT_H
#include "ct_defs.h"
#include "ctlapack.h"
#include "utilities.h"
#include "ctexceptions.h"
#include "Integrator.h"
#include "ResidJacEval.h"
#include "SquareMatrix.h"
#include "NonlinearSolver.h"
#include "mdp_allo.h"
#ifndef MAX
# define MAX(x,y) (( (x) > (y) ) ? (x) : (y))
#endif
#ifndef MIN
# define MIN(x,y) (( (x) < (y) ) ? (x) : (y))
#endif
#define OPT_SIZE 10
#define SUCCESS 0
#define FAILURE 1
#define STEADY 0
#define TRANSIENT 1
namespace Cantera {
enum BEulerMethodType {
BEulerFixedStep,
BEulerVarStep
};
/**
* Exception class thrown when a BEuler error is encountered.
*/
class BEulerErr : public CanteraError {
public:
BEulerErr(std::string msg);
};
#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.
*/
BEulerInt();
virtual ~BEulerInt();
virtual void setTolerances(double reltol, int n, double* abstol);
virtual void setTolerances(double reltol, double abstol);
virtual void setProblemType(int probtype);
virtual void initializeRJE(double t0, ResidJacEval& func);
virtual void reinitializeRJE(double t0, ResidJacEval& func);
virtual double integrateRJE(double tout, double tinit = 0.0);
virtual doublereal step(double tout);
virtual void setSolnWeights();
virtual double& solution(int k){ return m_y_n[k]; }
double* solution(){ return m_y_n; }
int nEquations() const { return m_neq;}
virtual int nEvals() const;
virtual void setMethodBEMT(BEulerMethodType t);
virtual void setIterator(IterType t);
virtual void setMaxStep(double hmax);
virtual void setMaxNumTimeSteps(int);
virtual void setNumInitialConstantDeltaTSteps(int);
void print_solnDelta_norm_contrib(const double * const soln0,
const char * const s0,
const double * const soln1,
const char * const s1,
const char * const title,
const double * const y0,
const double * const y1,
double damp,
int num_entries);
virtual void setPrintSolnOptions(int printSolnStepInterval,
int printSolnNumberToTout,
int printSolnFirstSteps = 0,
bool dumpJacobians = false);
void setNonLinOptions(int min_newt_its = 0,
bool matrixConditioning = false,
bool colScaling = false,
bool rowScaling = true);
virtual void setPrintFlag(int print_flag);
virtual void setColumnScales();
/**
* calculate the solution error norm
*/
virtual double soln_error_norm(const double * const,
bool printLargest = false);
virtual void setInitialTimeStep(double delta_t);
void beuler_jac(SquareMatrix &, double * const,
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.
*/
void calc_y_pred(int);
/**
* Internal function to calculate the time derivative at the
* new step
*/
void calc_ydot(int, double *, double *);
/**
* Internal function to calculate the time step truncation
* error for a predictor corrector time step
*/
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
*/
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
*
*
*/
int solve_nonlinear_problem(double * const y_comm,
double * const ydot_comm, double CJ,
double time_curr,
SquareMatrix& jac,
int &num_newt_its,
int &num_linear_solves,
int &num_backtracks,
int loglevelInput);
/**
* Compute the undamped Newton step. The residual function is
* evaluated at x, but the Jacobian is not recomputed.
*/
void doNewtonSolve(double, double *, double*, double *,
SquareMatrix&, int);
//! 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 boundStep(const double * const y, const double * const step0, int loglevel);
/*
* Damp step
*/
int dampStep(double, const double*, const double*,
const double *, double*, double*,
double*, double&, SquareMatrix&, int&, bool, int&);
/*
* Compute Residual Weights
*/
void computeResidWts(SquareMatrix &jac);
/*
* Filter a new step
*/
double filterNewStep(double, double *, double *);
/*
* get the next time to print out
*/
double getPrintTime(double time_current);
/********************** Member data ***************************/
/*********************
* METHOD FLAGS
*********************/
//! IterType is used to specify how the nonlinear equations are
//! to be relaxed at each time step.
IterType m_iter;
/**
* MethodType is used to specify how the time step is to be
* chosen. Currently, there are two choices, one is a fixed
* step method while the other is based on a predictor-corrector
* algorithm and a time-step truncation error tolerance.
*/
BEulerMethodType m_method;
/**
* m_jacFormMethod determines how a matrix is formed.
*/
int m_jacFormMethod;
/**
* m_rowScaling is a boolean. If true then row sum scaling
* of the Jacobian matrix is carried out when solving the
* linear systems.
*/
bool m_rowScaling;
/**
* m_colScaling is a boolean. If true, then column scaling
* is performed on each solution of the linear system.
*/
bool m_colScaling;
/**
* m_matrixConditioning is a boolean. If true, then the
* Jacobian and every rhs is multiplied by the inverse
* of a matrix that is suppose to reduce the condition
* number of the matrix. This is done before row scaling.
*/
bool m_matrixConditioning;
/**
* If m_itol =1 then each component has an individual
* value of atol. If m_itol = 0, the all atols are equal.
*/
int m_itol;
/**
* Relative time truncation error tolerances
*/
double m_reltol;
/**
* Absolute time truncation error tolerances, when uniform
* for all variables.
*/
double m_abstols;
/**
* Vector of absolute time truncation error tolerance
* when not uniform for all variables.
*/
double *m_abstol;
/**
* Error Weights. This is a surprisingly important quantity.
*/
double *m_ewt;
//! Maximum step size
double m_hmax;
/**
* Maximum integration order
*/
int m_maxord;
/**
* Current integration order
*/
int m_order;
/**
* Time step number
*/
int m_time_step_num;
int m_time_step_attempts;
/**
* 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
*/
int m_numInitialConstantDeltaTSteps;
/**
* Failure Counter -> keeps track of the number
* of consequetive failures
*/
int m_failure_counter;
/**
* Minimum Number of Newton Iterations per nonlinear step
* default = 0
*/
int m_min_newt_its;
/************************
* PRINTING OPTIONS
************************/
/**
* Step Interval at which to print out the solution
* default = 1;
* If set to zero, there is no printout
*/
int m_printSolnStepInterval;
/**
* Number of evenly spaced printouts of the solution
* If zero, there is no printout from this option
* default 1
* If set to zero there is no printout.
*/
int m_printSolnNumberToTout;
/**
* Number of initial steps that the solution is
* printed out.
* default = 0
*/
int m_printSolnFirstSteps;
/**
* Dump Jacobians to disk
* default false
*/
bool m_dumpJacobians;
/*********************
* INTERNAL SOLUTION VALUES
*********************/
/**
* Number of equations in the ode integrator
*/
int m_neq;
double *m_y_n;
double *m_y_nm1;
double *m_y_pred_n;
double *m_ydot_n;
double *m_ydot_nm1;
/************************
* TIME VARIABLES
************************/
/**
* Initial time at the start of the integration
*/
double m_t0;
/**
* Final time
*/
double m_time_final;
/**
*
*/
double time_n;
double time_nm1;
double time_nm2;
double delta_t_n;
double delta_t_nm1;
double delta_t_nm2;
double delta_t_np1;
/**
* Maximum permissible time step
*/
double delta_t_max;
double *m_resid;
double *m_residWts;
double *m_wksp;
ResidJacEval *m_func;
double *m_rowScales;
double *m_colScales;
/**
* Pointer to the jacobian representing the
* time dependent problem
*/
SquareMatrix *tdjac_ptr;
/**
* Determines the level of printing for each time
* step.
* 0 -> absolutely nothing is printed for
* a single time step.
* 1 -> One line summary per time step
* 2 -> short description, points of interest
* 3 -> Lots printed per time step (default)
*/
int m_print_flag;
/***************************************************************************
* COUNTERS OF VARIOUS KINDS
***************************************************************************/
/**
* 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
*/
int m_numTotalNewtIts;
/**
* Total number of linear iterations
*/
int m_numTotalLinearSolves;
/**
* Total number of convergence failures.
*/
int m_numTotalConvFails;
/**
* Total Number of time truncation error failures
*/
int m_numTotalTruncFails;
/*
*
*/
int num_failures;
};
} // namespace
#endif // CT_BEULER

4
Cantera/src/numerics/Makefile.in
View file

@ -35,14 +35,14 @@ CXX_FLAGS = @CXXFLAGS@ $(LOCAL_DEFS) $(CXX_OPT) $(PIC_FLAG) $(DEBUG_FLAG)
NUMERICS_OBJ = DenseMatrix.o funcs.o Func1.o \
ODE_integrators.o BandMatrix.o DAE_solvers.o \
funcs.o sort.o SquareMatrix.o ResidJacEval.o NonlinearSolver.o \
solveProb.o
solveProb.o BEulerInt.o
NUMERICS_H = ArrayViewer.h DenseMatrix.h \
funcs.h ctlapack.h Func1.h FuncEval.h \
polyfit.h\
BandMatrix.h Integrator.h DAE_Solver.h ResidEval.h sort.h \
SquareMatrix.h ResidJacEval.h NonlinearSolver.h \
solveProb.h
solveProb.h BEulerInt.h
ifeq ($(use_sundials), 1)
ODEPACKAGE_H = CVodesIntegrator.h

27
Cantera/src/numerics/ResidEval.h
View file

@ -110,6 +110,33 @@ namespace Cantera {
//! Return the number of equations in the equation system
virtual int nEquations() const = 0;
//! Write out to a file or to standard output the current solution
/*!
* ievent is a description of the event that caused this
* function to be called.
*/
virtual void writeSolution(int ievent, const double time,
const double deltaT,
const int time_step_num,
const double *y, const double *ydot) {
int k;
printf("ResidEval::writeSolution\n");
printf(" Time = %g, ievent = %d, deltaT = %g\n", time, ievent, deltaT);
if (ydot) {
printf(" k y[] ydot[]\n");
for (k = 0; k < nEquations(); k++) {
printf("%d %g %g\n", k, y[k], ydot[k]);
}
} else {
printf(" k y[]\n");
for (k = 0; k < nEquations(); k++) {
printf("%d %g \n", k, y[k]);
}
}
}
protected: