Commit 523db2d1 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

newton: rename: u->uLastIter, uOld->uCurrentIter, use separate memory for the...

newton: rename: u->uLastIter, uOld->uCurrentIter, use separate memory for the next iterative solution and the delta vector

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4861 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 6c33aa50
......@@ -524,9 +524,8 @@ public:
* \param oldGlobalSol The previous global solution
*/
void updateStaticData(SolutionVector &curGlobalSol,
SolutionVector &oldGlobalSol)
const SolutionVector &oldGlobalSol)
{
bool wasSwitched = false;
for (unsigned i = 0; i < staticVertexDat_.size(); ++i)
......@@ -554,8 +553,11 @@ public:
i,
false);
const GlobalPosition &global = it->geometry().corner(i);
wasSwitched = primaryVarSwitch_(curGlobalSol, volVars,
globalIdx, global) || wasSwitched;
if (primaryVarSwitch_(curGlobalSol,
volVars,
globalIdx,
global))
wasSwitched = true;
}
}
......
......@@ -74,14 +74,15 @@ public:
* iterations required or on the variable switch
*
* \param u The current global solution vector
* \param uOld The previous global solution vector
* \param uLastIter The previous global solution vector
*
*/
void newtonEndStep(SolutionVector &u, SolutionVector &uOld)
void newtonEndStep(SolutionVector &uCurrentIter,
const SolutionVector &uLastIter)
{
// call the method of the base class
this->method().model().updateStaticData(u, uOld);
ParentType::newtonEndStep(u, uOld);
this->method().model().updateStaticData(uCurrentIter, uLastIter);
ParentType::newtonEndStep(uCurrentIter, uLastIter);
}
/*!
......@@ -92,33 +93,36 @@ public:
*
* Different update strategies, such as line search and chopped
* updates can be implemented. The default behaviour is just to
* subtract deltaU from uOld.
* subtract deltaU from uLastIter.
*
* \param uCurrentIter The solution after the current Newton iteration
* \param uLastIter The solution after the last Newton iteration
* \param deltaU The delta as calculated from solving the linear
* system of equations. This parameter also stores
* the updated solution.
* \param uOld The solution of the last iteration
*/
void newtonUpdate(SolutionVector &deltaU, const SolutionVector &uOld)
void newtonUpdate(SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
this->writeConvergence_(uOld, deltaU);
this->writeConvergence_(uLastIter, deltaU);
this->newtonUpdateRelError(uOld, deltaU);
this->newtonUpdateRelError(uLastIter, deltaU);
// compute the vertex and element colors for partial
// reassembly
if (enablePartialReassemble) {
Scalar reassembleTol = Dumux::geometricMean(this->error_, 0.1*this->tolerance_);
reassembleTol = std::max(reassembleTol, 0.1*this->tolerance_);
this->model_().jacobianAssembler().updateDiscrepancy(uOld, deltaU);
this->model_().jacobianAssembler().updateDiscrepancy(uLastIter, deltaU);
this->model_().jacobianAssembler().computeColors(reassembleTol);
}
if (GET_PROP_VALUE(TypeTag, PTAG(NewtonUseLineSearch)))
lineSearchUpdate_(deltaU, uOld);
lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
else {
deltaU *= - 1.0;
deltaU += uOld;
uCurrentIter = uLastIter;
uCurrentIter -= deltaU;
}
}
......@@ -137,32 +141,32 @@ public:
};
private:
void lineSearchUpdate_(SolutionVector &u, const SolutionVector &uOld)
void lineSearchUpdate_(SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
Scalar lambda = 1.0;
Scalar globDef;
SolutionVector tmp(u);
Scalar oldGlobDef = this->method().model().globalResidual(tmp);
int n = 0;
// calculate the residual of the current solution
SolutionVector tmp(uLastIter);
Scalar oldGlobDef = this->method().model().globalResidual(tmp, uLastIter);
while (true) {
u *= -lambda;
u += uOld;
uCurrentIter = deltaU;
uCurrentIter *= -lambda;
uCurrentIter += uLastIter;
globDef = this->method().model().globalResidual(tmp);
// calculate the residual of the current solution
globDef = this->method().model().globalResidual(tmp, uCurrentIter);
if (globDef < oldGlobDef || lambda <= 1.0/8) {
this->endIterMsg() << ", defect " << oldGlobDef << "->" << globDef << "@lambda=2^-" << n;
this->endIterMsg() << ", defect " << oldGlobDef << "->" << globDef << "@lambda=" << lambda;
return;
}
// undo the last iteration
u -= uOld;
u /= - lambda;
// try with a smaller update
lambda /= 2;
++n;
}
};
......
......@@ -72,33 +72,33 @@ public:
* This is basically the step
* \f[ u^{k+1} = u^k - \Delta u^k \f]
*
* \param deltaU When the method is called, this contains the
* vector of differences between the current
* iterative solution and the next \f$\Delta
* u\f$. After the method is finished it should
* contain the next iterative solution \f$ u^{k+1} \f$.
* \param uOld The current iterative solution \f$ u^k \f$
* \param uCurrentIter The solution after the current Newton iteration \f$ u^{k+1} \f$
* \param uLastIter The solution after the last Newton iteration \f$ u^k \f$
* \param deltaU The vector of differences between the last
* iterative solution and the next one \f$ \Delta u^k \f$
*/
void newtonUpdate(SolutionVector &deltaU, const SolutionVector &uOld)
void newtonUpdate(SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
this->writeConvergence_(uOld, deltaU);
this->newtonUpdateRelError(uOld, deltaU);
this->writeConvergence_(uLastIter, deltaU);
this->newtonUpdateRelError(uLastIter, deltaU);
// compute the vertex and element colors for partial
// reassembly
if (enablePartialReassemble) {
Scalar reassembleTol = Dumux::geometricMean(this->error_, 0.1*this->tolerance_);
reassembleTol = std::max(reassembleTol, 0.1*this->tolerance_);
this->model_().jacobianAssembler().updateDiscrepancy(uOld, deltaU);
this->model_().jacobianAssembler().updateDiscrepancy(uLastIter, deltaU);
this->model_().jacobianAssembler().computeColors(reassembleTol);
}
if (GET_PROP_VALUE(TypeTag, PTAG(NewtonUseLineSearch)))
lineSearchUpdate_(deltaU, uOld);
lineSearchUpdate_(uCurrentIter, uLastIter, deltaU);
else {
// update the solution vector
deltaU *= -1;
deltaU += uOld;
uCurrentIter = uLastIter;
uCurrentIter -= deltaU;
// clamp saturation change to at most 20% per iteration
FVElementGeometry fvElemGeom;
......@@ -114,7 +114,7 @@ public:
const SpatialParameters &sp = this->problem_().spatialParameters();
const MaterialLawParams &mp = sp.materialLawParams(*eIt, fvElemGeom, i);
Scalar pcMin = MaterialLaw::pC(mp, 1.0);
Scalar pW = uOld[globI][pwIdx];
Scalar pW = uLastIter[globI][pwIdx];
Scalar pN = std::max(this->problem_().referencePressure(*eIt, fvElemGeom, i),
pW + pcMin);
Scalar pcOld = pN - pW;
......@@ -126,38 +126,40 @@ public:
Scalar pwMax = pN - MaterialLaw::pC(mp, SwOld + 0.2);
// clamp the result
deltaU[globI][pwIdx] = std::max(pwMin, std::min(deltaU[globI][pwIdx], pwMax));
pW = uCurrentIter[globI][pwIdx];
pW = std::max(pwMin, std::min(pW, pwMax));
uCurrentIter[globI][pwIdx] = pW;
}
}
}
}
private:
void lineSearchUpdate_(SolutionVector &u, const SolutionVector &uOld)
void lineSearchUpdate_(SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
Scalar lambda = 1.0;
Scalar globDef;
SolutionVector tmp(u);
Scalar oldGlobDef = this->model_().globalResidual(tmp, uOld);
SolutionVector tmp(uLastIter);
Scalar oldGlobDef = this->model_().globalResidual(tmp, uLastIter);
int n = 0;
while (true) {
u *= -lambda;
u += uOld;
globDef = this->model_().globalResidual(tmp);
uCurrentIter = deltaU;
uCurrentIter *= -lambda;
uCurrentIter += uLastIter;
// calculate the residual of the current solution
globDef = this->model_().globalResidual(tmp, uCurrentIter);
if (globDef < oldGlobDef || lambda <= 1.0/64) {
this->endIterMsg() << ", defect " << oldGlobDef << "->" << globDef << "@lambda=2^-" << n;
this->endIterMsg() << ", defect " << oldGlobDef << "->" << globDef << "@lambda=" << lambda;
return;
}
// undo the last iteration
u -= uOld;
u /= - lambda;
// try with a smaller update
lambda /= 2;
++n;
}
};
};
......
......@@ -157,10 +157,10 @@ struct NewtonConvergenceWriter
gv);
};
void writeFields(const SolutionVector &uOld,
void writeFields(const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
ctl_.method().model().addConvergenceVtkFields(*vtkMultiWriter_, uOld, deltaU);
ctl_.method().model().addConvergenceVtkFields(*vtkMultiWriter_, uLastIter, deltaU);
};
void endIteration()
......@@ -204,7 +204,7 @@ struct NewtonConvergenceWriter<TypeTag, false>
void beginIteration(const GridView &gv)
{ };
void writeFields(const SolutionVector &uOld,
void writeFields(const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{ };
......@@ -323,7 +323,7 @@ public:
*
* \param u The current solution
*/
bool newtonProceed(const SolutionVector &u)
bool newtonProceed(const SolutionVector &uCurrentIter)
{
if (numSteps_ < rampUpSteps() + 2)
return true; // we always do at least two iterations
......@@ -356,7 +356,7 @@ public:
* \param method The object where the NewtonMethod is executed
* \param u The initial solution
*/
void newtonBegin(NewtonMethod &method, SolutionVector &u)
void newtonBegin(NewtonMethod &method, const SolutionVector &u)
{
method_ = &method;
numSteps_ = 0;
......@@ -399,10 +399,10 @@ public:
* The relative error can be seen as a norm of the difference
* between the current and the next iteration.
*
* \param uOld The current iterative solution
* \param uLastIter The current iterative solution
* \param deltaU The difference between the current and the next solution
*/
void newtonUpdateRelError(const SolutionVector &uOld,
void newtonUpdateRelError(const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
// calculate the relative error as the maximum relative
......@@ -412,12 +412,12 @@ public:
int idxI = -1;
int aboveTol = 0;
for (int i = 0; i < int(uOld.size()); ++i) {
PrimaryVariables uNewI = uOld[i];
for (int i = 0; i < int(uLastIter.size()); ++i) {
PrimaryVariables uNewI = uLastIter[i];
uNewI -= deltaU[i];
Scalar vertErr =
model_().relativeErrorVertex(i,
uOld[i],
uLastIter[i],
uNewI);
if (vertErr > tolerance_)
......@@ -432,18 +432,18 @@ public:
}
/*!
* \brief Solve the linear system of equations \f$\mathbf{A}u - b = 0\f$.
* \brief Solve the linear system of equations \f$\mathbf{A}x - b = 0\f$.
*
* Throws Dumux::NumericalProblem if the linear solver didn't
* converge.
*
* \param A The matrix of the linear system of equations
* \param u The vector which solves the linear system
* \param x The vector which solves the linear system
* \param b The right hand side of the linear system
*/
template <class Vector>
void newtonSolveLinear(const JacobianMatrix &A,
Vector &u,
Vector &x,
const Vector &b)
{
// if the deflection of the newton method is large, we do not
......@@ -454,7 +454,7 @@ public:
Scalar residReduction = 1e-6;
try {
solveLinear_(A, u, b, residReduction);
solveLinear_(A, x, b, residReduction);
// make sure all processes converged
int converged = 1;
......@@ -487,40 +487,44 @@ public:
*
* Different update strategies, such as line search and chopped
* updates can be implemented. The default behaviour is just to
* subtract deltaU from uOld, i.e.
* subtract deltaU from uLastIter, i.e.
* \f[ u^{k+1} = u^k - \Delta u^k \f]
*
* \param uCurrentIter The solution vector after the current iteration
* \param uLastIter The solution vector after the last iteration
* \param deltaU The delta as calculated from solving the linear
* system of equations. This parameter also stores
* the updated solution.
* \param uOld The solution vector of the last iteration
*/
void newtonUpdate(SolutionVector &deltaU, const SolutionVector &uOld)
void newtonUpdate(SolutionVector &uCurrentIter,
const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
writeConvergence_(uOld, deltaU);
writeConvergence_(uLastIter, deltaU);
newtonUpdateRelError(uOld, deltaU);
newtonUpdateRelError(uLastIter, deltaU);
// compute the vertex and element colors for partial
// reassembly
if (enablePartialReassemble) {
Scalar reassembleTol = Dumux::geometricMean(error_, 0.1*tolerance_);
reassembleTol = std::max(reassembleTol, 0.1*tolerance_);
this->model_().jacobianAssembler().updateDiscrepancy(uOld, deltaU);
this->model_().jacobianAssembler().updateDiscrepancy(uLastIter, deltaU);
this->model_().jacobianAssembler().computeColors(reassembleTol);
}
deltaU *= -1;
deltaU += uOld;
uCurrentIter = uLastIter;
uCurrentIter -= deltaU;
}
/*!
* \brief Indicates that one newton iteration was finished.
*
* \param u The solution after the current iteration
* \param uOld The solution at the beginning of the current iteration
* \param uLastIter The solution at the beginning of the current iteration
*/
void newtonEndStep(SolutionVector &u, SolutionVector &uOld)
void newtonEndStep(const SolutionVector &uCurrentIter,
const SolutionVector &uLastIter)
{
++numSteps_;
......@@ -677,11 +681,11 @@ protected:
const Implementation &asImp_() const
{ return *static_cast<const Implementation*>(this); }
void writeConvergence_(const SolutionVector &uOld,
void writeConvergence_(const SolutionVector &uLastIter,
const SolutionVector &deltaU)
{
convergenceWriter_.beginIteration(this->gridView_());
convergenceWriter_.writeFields(uOld, deltaU);
convergenceWriter_.writeFields(uLastIter, deltaU);
convergenceWriter_.endIteration();
};
......
......@@ -115,8 +115,10 @@ public:
protected:
bool execute_(NewtonController &ctl)
{
// TODO (?): u shouldn't be hard coded to the model
SolutionVector &u = model().curSol();
SolutionVector &uCurrentIter = model().curSol();
SolutionVector uLastIter(uCurrentIter);
SolutionVector deltaU(uCurrentIter);
JacobianAssembler &jacobianAsm = model().jacobianAssembler();
Dune::Timer assembleTimer(false);
......@@ -124,30 +126,39 @@ protected:
Dune::Timer updateTimer(false);
// tell the controller that we begin solving
ctl.newtonBegin(*this, u);
ctl.newtonBegin(*this, uCurrentIter);
// execute the method as long as the controller thinks
// that we should do another iteration
while (ctl.newtonProceed(u))
while (ctl.newtonProceed(uCurrentIter))
{
// notify the controller that we're about to start
// a new timestep
ctl.newtonBeginStep();
// make the current solution to the old one
uOld_ = u;
uLastIter = uCurrentIter;
if (ctl.verbose()) {
std::cout << "Assembling global jacobian";
std::cout.flush();
}
///////////////
// assemble
///////////////
assembleTimer.start();
// linearize the problem at the current solution
assembleTimer.start();
jacobianAsm.assemble();
assembleTimer.stop();
// solve the resultuing linear equation system
///////////////
// linear solve
///////////////
// solve the resulting linear equation system
solveTimer.start();
if (ctl.verbose()) {
std::cout << "\rSolve Mx = r";
// Clear the current line using an ansi escape
......@@ -158,36 +169,36 @@ protected:
std::cout.flush();
}
solveTimer.start();
// set the delta vector to zero before solving the linear system!
u = 0;
deltaU = 0;
// ask the controller to solve the linearized system
ctl.newtonSolveLinear(jacobianAsm.matrix(),
u,
deltaU,
jacobianAsm.residual());
solveTimer.stop();
///////////////
// update
///////////////
updateTimer.start();
// update the current solution (i.e. uOld) with the delta
// (i.e. u). The result is stored in u
ctl.newtonUpdate(u, uOld_);
ctl.newtonUpdate(uCurrentIter, uLastIter, deltaU);
updateTimer.stop();
// tell the controller that we're done with this iteration
ctl.newtonEndStep(u, uOld_);
ctl.newtonEndStep(uCurrentIter, uLastIter);
}
// tell the controller that we're done
ctl.newtonEnd();
Scalar elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
if (ctl.verbose()) {
std::cout << "Timings assemble/solve/update: "
Scalar elapsedTot = assembleTimer.elapsed() + solveTimer.elapsed() + updateTimer.elapsed();
std::cout << "Assemble/solve/update time: "
<< assembleTimer.elapsed() << "(" << 100*assembleTimer.elapsed()/elapsedTot << "%)/"
<< solveTimer.elapsed() << "(" << 100*solveTimer.elapsed()/elapsedTot << "%)/"
<< updateTimer.elapsed() << "(" << 100*updateTimer.elapsed()/elapsedTot << "%)"
//<< " localEvals: " << problem_.model().localResidual().elapsed() << "("<<100*problem_.model().localResidual().elapsed()/assembleTimer.elapsed()<<"% of asm)"
//<< " singleEvalTime: " << problem_.model().localResidual().evalTimeSingle() << " sec"
<< "\n";
}
......@@ -201,8 +212,6 @@ protected:
}
private:
SolutionVector uOld_;
Problem &problem_;
};
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment