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

Iteration on getting this class to work.

Browse source 
There was a member function that wasn't even defined.
This commit is contained in:
Harry Moffat 2010-07-27 01:48:00 +00:00
parent 74f33b8829
commit 36add9bc08
3 changed files with 210 additions and 6 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -79,6 +79,8 @@ namespace Cantera {
m_numTotalNewtIts(0),
m_min_newt_its(0),
filterNewstep(0),
m_jacFormMethod(NSOLN_JAC_NUM),
m_nJacEval(0),
time_n(0.0),
m_matrixConditioning(0),
m_order(1),
@ -104,7 +106,25 @@ namespace Cantera {
}
}
NonlinearSolver::NonlinearSolver(const NonlinearSolver &right) {
NonlinearSolver::NonlinearSolver(const NonlinearSolver &right) :
m_func(right.m_func),
neq_(0),
delta_t_n(-1.0),
m_nfe(0),
m_colScaling(0),
m_rowScaling(0),
m_numTotalLinearSolves(0),
m_numTotalNewtIts(0),
m_min_newt_its(0),
filterNewstep(0),
m_jacFormMethod(NSOLN_JAC_NUM),
m_nJacEval(0),
time_n(0.0),
m_matrixConditioning(0),
m_order(1),
rtol_(1.0E-3),
atolBase_(1.0E-10)
{
*this =operator=(right);
}
@ -137,6 +157,8 @@ namespace Cantera {
m_numTotalNewtIts = right.m_numTotalNewtIts;
m_min_newt_its = right.m_min_newt_its;
filterNewstep = right.filterNewstep;
m_jacFormMethod = right.m_jacFormMethod;
m_nJacEval = right.m_nJacEval;
time_n = right.time_n;
m_matrixConditioning = right.m_matrixConditioning;
m_order = right.m_order;
@ -146,7 +168,7 @@ namespace Cantera {
return *this;
}
//====================================================================================================================
// Create solution weights for convergence criteria
/*
* We create soln weights from the following formula
@ -163,7 +185,7 @@ namespace Cantera {
m_ewt[i] = rtol_ * fabs(y[i]) + atolk_[i];
}
}
//====================================================================================================================
// set bounds constraints for all variables in the problem
/*
*
@ -966,6 +988,164 @@ namespace Cantera {
mdp::mdp_safe_free((void **) &imax);
}
//================================================================================================
/*
* 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.
* It also returns an exact zero if the difference is less than
* 1.0E-300.
*
* returns: a - b
*
* This routine is used in numerical differencing schemes in order
* to avoid roundoff errors resulting in creating Jacobian terms.
* Note: This is a slow routine. However, jacobian errors may cause
* loss of convergence. Therefore, in practice this routine
* has proved cost-effective.
*/
static inline double subtractRD(double a, double b) {
double diff = a - b;
double d = MIN(fabs(a), fabs(b));
d *= 1.0E-14;
double ad = fabs(diff);
if (ad < 1.0E-300) {
diff = 0.0;
}
if (ad < d) {
diff = 0.0;
}
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 NonlinearSolver::beuler_jac(SquareMatrix &J, double * const f,
double time_curr, double CJ,
double * const y,
double * const ydot,
int num_newt_its)
{
int i, j;
double* col_j;
double ysave, ydotsave, dy;
/*
* Clear the factor flag
*/
J.clearFactorFlag();
if (m_jacFormMethod == NSOLN_JAC_ANAL) {
/********************************************************************
* Call the function to get a jacobian.
*/
m_func->evalJacobian(time_curr, delta_t_n, y, ydot, J, f);
#ifdef DEBUG_HKM
//double dddd = J(89, 89);
//checkFinite(dddd);
#endif
m_nJacEval++;
m_nfe++;
} else {
/*******************************************************************
* Generic algorithm to calculate a numerical Jacobian
*/
/*
* Calculate the current value of the rhs given the
* current conditions.
*/
m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, f);
m_nfe++;
m_nJacEval++;
/*
* Malloc a vector and call the function object to return a set of
* deltaY's that are appropriate for calculating the numerical
* derivative.
*/
double *dyVector = mdp::mdp_alloc_dbl_1(neq_, MDP_DBL_NOINIT);
m_func->calcDeltaSolnVariables(time_curr, y, ydot, dyVector, DATA_PTR(m_ewt));
#ifdef DEBUG_HKM
bool print_NumJac = false;
if (print_NumJac) {
FILE *idy = fopen("NumJac.csv", "w");
fprintf(idy, "Unk m_ewt y "
"dyVector ResN\n");
for (int iii = 0; iii < neq_; iii++){
fprintf(idy, " %4d %16.8e %16.8e %16.8e %16.8e \n",
iii, m_ewt[iii], y[iii], dyVector[iii], f[iii]);
}
fclose(idy);
}
#endif
/*
* Loop over the variables, formulating a numerical derivative
* of the dense matrix.
* For the delta in the variable, we will use a variety of approaches
* The original approach was to use the error tolerance amount.
* This may not be the best approach, as it could be overly large in
* some instances and overly small in others.
* We will first protect from being overly small, by using the usual
* sqrt of machine precision approach, i.e., 1.0E-7,
* to bound the lower limit of the delta.
*/
for (j = 0; j < neq_; j++) {
/*
* Get a pointer into the column of the matrix
*/
col_j = (double *) J.ptrColumn(j);
ysave = y[j];
dy = dyVector[j];
//dy = fmaxx(1.0E-6 * m_ewt[j], fabs(ysave)*1.0E-7);
y[j] = ysave + dy;
dy = y[j] - ysave;
ydotsave = ydot[j];
ydot[j] += dy * CJ;
/*
* Call the functon
*/
m_func->evalResidNJ(time_curr, delta_t_n, y, ydot, DATA_PTR(m_y_nm1),
true, j, dy);
m_nfe++;
double diff;
for (i = 0; i < neq_; i++) {
diff = subtractRD(m_y_nm1[i], f[i]);
col_j[i] = diff / dy;
//col_j[i] = (m_wksp[i] - f[i])/dy;
}
y[j] = ysave;
ydot[j] = ydotsave;
}
/*
* Release memory
*/
mdp::mdp_safe_free((void **) &dyVector);
}
}
}

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

@ -23,12 +23,19 @@
namespace Cantera {
// I think steady state is the only option I'm gunning for
#define NSOLN_TYPE_PSEUDO_TIME_DEPENDENT 2
#define NSOLN_TYPE_TIME_DEPENDENT 1
#define NSOLN_TYPE_STEADY_STATE 0
#define NSOLN_JAC_NUM 1
#define NSOLN_JAC_ANAL 2
//! Class that calculates the solution to a nonlinear system
/*!
* UNDER CONSTRUCTION - do not use!!!!!!!!!!!!!!!!!!!!!!!!!
*
* @ingroup numerics
*/
@ -192,8 +199,16 @@ namespace Cantera {
*/
void calc_ydot(int order, double * const y_curr, double * const ydot_curr);
void beuler_jac(SquareMatrix &, double * const,
double, double, double * const, double * const, int);
//! Function called to evaluate the jacobian matrix and the curent
//! residual vector.
/*!
*
*
*/
void beuler_jac(SquareMatrix &J, double * const f,
double time_curr, double CJ, double * const y,
double * const ydot, int num_newt_its);
double filterNewStep(double, double *, double *);
@ -311,6 +326,16 @@ namespace Cantera {
int filterNewstep;
//! Jacobian formation method
/*!
* 1 = numerical (default)
* 2 = analytical
*/
int m_jacFormMethod;
int m_nJacEval;
//! Current system time
/*!
* Note, we assume even for steady state problems that the residual

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

@ -158,7 +158,6 @@ namespace Cantera {
* Backwards Euler is assumed.
*/
virtual void evalJacobian(const doublereal t, const doublereal deltaT,
const double* const y,
const double* const ydot,
SquareMatrix &J,