Commit 29c68ffa authored by Andreas Lauser's avatar Andreas Lauser
Browse files

*vangenuchten, box jacobian assembler, box local jacobian, box local residual: various changes

- box jacobian assembler: Do no longer use the PDELab grid operator space to assemble the jacobian.
- vangenuchten: fix typos to make it compile again
- box local jacobian: evaluate local residual inside local jacobian to increase efficiency slightly
- box local residual: correct comment about fluxes

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4693 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 7c6132c3
......@@ -592,11 +592,10 @@ public:
int numEdges; //!< number of edges
int numFaces; //!< number of faces (0 in < 3D)
const LocalFEMSpace feMap_;
static const LocalFEMSpace feMap_;
bool computeGradientAtScvCenters;
BoxFVElementGeometry()
: feMap_()
{
computeGradientAtScvCenters = false;
}
......@@ -843,6 +842,9 @@ public:
}
};
template<class TypeTag>
const typename BoxFVElementGeometry<TypeTag>::LocalFEMSpace BoxFVElementGeometry<TypeTag>::feMap_;
}
#endif
......
......@@ -126,6 +126,7 @@ public:
// set the current grid element and update the element's
// finite volume geometry
elemPtr_ = &element;
fvElemGeom_.update(gridView_(), element);
reset_();
int elemColor = jacAsm_().elementColor(element);
......@@ -134,7 +135,6 @@ public:
return;
}
fvElemGeom_.update(gridView_(), elem_());
bcTypes_.update(problem_(), elem_(), fvElemGeom_);
// this is pretty much a HACK because the internal state of
......@@ -156,17 +156,13 @@ public:
elem_(),
fvElemGeom_,
false /* isOldSol? */);
if (numDiffMethod != 0) {
// for forward and backward differences we need the local
// residual at the evaluation point. For efficiency we
// only calculate this once.
// calculate the local residual
localResidual().eval(elem_(),
fvElemGeom_,
prevVolVars_,
curVolVars_,
bcTypes_);
evalPointResid_ = localResidual().residual();
};
residual_ = localResidual().residual();
// calculate the local jacobian matrix
ElementSolutionVector partialDeriv(numVertices);
......@@ -186,13 +182,15 @@ public:
}
/*!
* \brief Returns a reference to the local residual.
* \brief Returns a reference to the object which calculates the
* local residual.
*/
const LocalResidual &localResidual() const
{ return localResidual_; }
/*!
* \brief Returns a reference to the local residual.
* \brief Returns a reference to the object which calculates the
* local residual.
*/
LocalResidual &localResidual()
{ return localResidual_; }
......@@ -209,6 +207,15 @@ public:
const MatrixBlock &mat(int i, int j) const
{ return A_[i][j]; }
/*!
* \brief Returns the residual of the equations at vertex i.
*
* \param i The local vertex (or sub-contol volume) index on which
* the equations are defined
*/
const PrimaryVariables &residual(int i) const
{ return residual_[i]; }
protected:
Implementation &asImp_()
{ return *static_cast<Implementation*>(this); }
......@@ -312,7 +319,7 @@ protected:
}
else {
// backward differences
dest = evalPointResid_;
dest = residual_;
}
......@@ -337,7 +344,7 @@ protected:
}
else {
// forward differences
dest -= evalPointResid_;
dest -= residual_;
}
dest /= delta;
......@@ -395,8 +402,8 @@ protected:
}
const Element *elemPtr_;
FVElementGeometry fvElemGeom_;
ElementBoundaryTypes bcTypes_;
// The problem we would like to solve
......@@ -406,11 +413,12 @@ protected:
// levels
ElementVolumeVariables prevVolVars_;
ElementVolumeVariables curVolVars_;
ElementSolutionVector evalPointResid_; // residual at evaluation point
LocalResidual localResidual_;
LocalBlockMatrix A_;
};
ElementSolutionVector residual_;
};
}
#endif
......@@ -119,19 +119,22 @@ public:
*/
void eval(const Element &element)
{
FVElementGeometry fvGeom;
fvGeom.update(gridView_(), element);
FVElementGeometry fvElemGeom;
fvElemGeom.update(gridView_(), element);
fvElemGeomPtr_ = &fvElemGeom;
ElementVolumeVariables volVarsPrev, volVarsCur;
volVarsPrev.update(problem_(),
element,
fvGeom,
fvElemGeom_(),
true /* oldSol? */);
volVarsCur.update(problem_(),
element,
fvGeom,
fvElemGeom_(),
false /* oldSol? */);
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeom);
bcTypes.update(problem_(), element, fvElemGeom_());
// this is pretty much a HACK because the internal state of
// the problem is not supposed to be changed during the
......@@ -140,7 +143,7 @@ public:
// not thread save.) The real solution are context objects!
problem_().updateCouplingParams(element);
asImp_().eval(element, fvGeom, volVarsPrev, volVarsCur, bcTypes);
asImp_().eval(element, fvElemGeom_(), volVarsPrev, volVarsCur, bcTypes);
}
/*!
......@@ -154,21 +157,26 @@ public:
*/
void evalStorage(const Element &element)
{
FVElementGeometry fvGeom;
fvGeom.update(gridView_(), element);
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeom);
ElementVolumeVariables volVars;
volVars.update(problem_(), element, fvGeom, false);
elemPtr_ = &element;
residual_.resize(fvGeom.numVertices);
residual_ = 0;
FVElementGeometry fvElemGeom;
fvElemGeom.update(gridView_(), element);
fvElemGeomPtr_ = &fvElemGeom;
elemPtr_ = &element;
fvElemGeomPtr_ = &fvGeom;
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvElemGeom_());
bcTypesPtr_ = &bcTypes;
// no previous volume variables!
prevVolVarsPtr_ = 0;
ElementVolumeVariables volVars;
volVars.update(problem_(), element, fvElemGeom_(), false);
curVolVarsPtr_ = &volVars;
residual_.resize(fvElemGeom_().numVertices);
residual_ = 0;
asImp_().evalStorage_();
}
......@@ -183,23 +191,21 @@ public:
void evalFluxes(const Element &element,
const ElementVolumeVariables &curVolVars)
{
FVElementGeometry fvGeom;
fvGeom.update(gridView_(), element);
elemPtr_ = &element;
fvElemGeomPtr_ = &model_().fvElemGeom(element);
ElementBoundaryTypes bcTypes;
bcTypes.update(problem_(), element, fvGeom);
bcTypes.update(problem_(), element, fvElemGeom_());
residual_.resize(fvGeom.numVertices);
residual_.resize(fvElemGeom_().numVertices);
residual_ = 0;
elemPtr_ = &element;
fvElemGeomPtr_ = &fvGeom;
bcTypesPtr_ = &bcTypes;
prevVolVarsPtr_ = 0;
curVolVarsPtr_ = &curVolVars;
asImp_().evalFluxes_();
}
/*!
* \brief Compute the local residual, i.e. the deviation of the
* equations from zero.
......@@ -418,15 +424,16 @@ protected:
//
// dStorage/dt = Flux + Source
//
// Re-arranging this, we get
// where the 'Flux' and the 'Source' terms represent the
// mass per second which _ENTER_ the finite
// volume. Re-arranging this, we get
//
// dStorage/dt - Source - Flux = 0
//
// The residual already contains the "dStorage/dt -
// Source" term, so we have to subtract the flux
// term. Since the calculated flux goes from sub-control
// volume i to sub-control volume j, we need to add the
// flux to j and subtract it from i
// Since the flux calculated by computeFlux() goes _OUT_
// of sub-control volume i and _INTO_ sub-control volume
// j, we need to add the flux to finite volume i and
// subtract it from finite volume j
residual_[i] += flux;
residual_[j] -= flux;
}
......@@ -613,6 +620,7 @@ protected:
const ElementBoundaryTypes *bcTypesPtr_;
};
}
#endif
......@@ -25,7 +25,6 @@
#include "pdelabboxlocaloperator.hh"
namespace Dumux {
namespace PDELab {
/*!
......@@ -42,6 +41,7 @@ class BoxAssembler
typedef typename GET_PROP_TYPE(TypeTag, PTAG(VertexMapper)) VertexMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ElementMapper)) ElementMapper;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(FVElementGeometry)) FVElementGeometry;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Constraints)) Constraints;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ScalarGridFunctionSpace)) ScalarGridFunctionSpace;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridFunctionSpace)) GridFunctionSpace;
......@@ -60,6 +60,7 @@ class BoxAssembler
typedef typename GridView::IntersectionIterator IntersectionIterator;
typedef typename GridView::template Codim<dim>::Entity Vertex;
typedef typename GridView::template Codim<dim>::EntityPointer VertexPointer;
typedef typename GridView::template Codim<dim>::Iterator VertexIterator;
typedef SolutionVector Vector;
......@@ -123,7 +124,6 @@ public:
*/
void init(Problem& problem)
{
inJacobianAssemble_ = false;
problemPtr_ = &problem;
fem_ = new FEM();
//cn_ = new Constraints(*problemPtr_);
......@@ -149,16 +149,6 @@ public:
*matrix_ = 0;
reuseMatrix_ = false;
// calculate the ghost vertices
VertexIterator vIt = gridView_().template begin<dim>();
VertexIterator vEndIt = gridView_().template end<dim>();
for (; vIt != vEndIt; ++vIt) {
if (vIt->partitionType() == Dune::GhostEntity) {
int vIdx = vertexMapper_().map(*vIt);
ghostIndices_.push_back(vIdx);
}
};
int numVerts = gridView_().size(dim);
int numElems = gridView_().size(0);
residual_.resize(numVerts);
......@@ -196,8 +186,26 @@ public:
*/
void assemble()
{
(*matrix_) = 0.0;
residual_ = 0.0;
ElementIterator elemIt = gridView_().template begin<0>();
ElementIterator elemEndIt = gridView_().template end<0>();
for (; elemIt != elemEndIt; ++elemIt) {
const Element &elem = *elemIt;
if (elem.partitionType() == Dune::GhostEntity)
assembleGhostElement_(elem);
else
assembleElement_(elem);
};
return;
#if 0
// assemble the global jacobian matrix
if (!reuseMatrix_) {
// we actually need to reassemle!
if (enablePartialReassemble) {
// move the evaluation points of red vertices
for (int i = 0; i < vertexColor_.size(); ++i) {
......@@ -208,18 +216,16 @@ public:
reassembleTolerance_ = nextReassembleTolerance_;
// we actually need to reassemle!
resetMatrix_();
inJacobianAssemble_ = true;
gridOperatorSpace_->jacobian(model_().curSol(), *matrix_);
inJacobianAssemble_ = false;
assemble_();
}
reuseMatrix_ = false;
// calculate the global residual
residual_ = 0;
gridOperatorSpace_->residual(model_().curSol(), residual_);
//gridOperatorSpace_->residual(model_().curSol(), residual_);
#if 0
typedef typename Matrix::block_type BlockType;
// set the entries for the ghost nodes
BlockType Id(0.0);
......@@ -234,6 +240,8 @@ public:
residual_[globI] = 0;
model_().curSol()[globI] = 0;
}
#endif
#endif
}
/*!
......@@ -484,9 +492,53 @@ public:
private:
// assemble a non-ghost element
void assembleElement_(const Element &elem)
{
model_().localJacobian().assemble(elem);
int numVertices = elem.template count<dim>();
for (int i=0; i < numVertices; ++ i) {
int globI = vertexMapper_().map(elem, i, dim);
// update the global residual
residual_[globI] += model_().localJacobian().residual(i);
// update the global jacobian
for (int j=0; j < numVertices; ++ j) {
int globJ = vertexMapper_().map(elem, j, dim);
(*matrix_)[globI][globJ] +=
model_().localJacobian().mat(i,j);
}
}
}
// "assemble" a ghost element
void assembleGhostElement_(const Element &elem)
{
int n = elem.template count<dim>();
for (int i=0; i < n; ++i) {
const VertexPointer vp = elem.template subEntity<dim>(i);
if (vp->partitionType() != Dune::GhostEntity)
continue; // ignore a ghost cell's non-ghost vertices
// set main diagonal entries for the vertex
int vIdx = vertexMapper_().map(*vp);
typedef typename Matrix::block_type BlockType;
BlockType &J = (*matrix_)[vIdx][vIdx];
for (int j = 0; j < BlockType::rows; ++j)
J[j][j] = 1.0;
// set residual for the vertex
residual_[vIdx] = 0;
}
}
void resetMatrix_()
{
if (!enablePartialReassemble) {
// If partial reassembly of the jacobian is not enabled,
// we can just reset everything!
(*matrix_) = 0;
return;
}
......@@ -534,7 +586,6 @@ private:
Matrix *matrix_;
bool reuseMatrix_;
bool inJacobianAssemble_;
SolutionVector evalPoint_;
std::vector<EntityColor> vertexColor_;
std::vector<EntityColor> elementColor_;
......
......@@ -111,8 +111,13 @@ int startFromDGF(int argc, char **argv)
// create grid
// -> load the grid from file
GridPointer gridPtr(dgfFileName);
if (mpiHelper.size() > 1)
if (mpiHelper.size() > 1) {
if (!Dune::Capabilities::isParallel<Grid>::v) {
std::cerr << "WARNING: THE PROGRAM IS STARTED USING MPI, BUT THE GRID IMPLEMENTATION\n"
<< " YOU HAVE CHOSEN IS NOT PARALLEL!\n";
}
gridPtr->loadBalance();
}
// instantiate and run the concrete problem
TimeManager timeManager;
......
......@@ -249,7 +249,7 @@ public:
return sp.eval(Swe);
}
return VanGenuchten::krw(params, Sw);
return VanGenuchten::krw(params, Swe);
};
/*!
......@@ -285,7 +285,7 @@ public:
return sp.eval(Swe);
}
return VanGenuchten::krn(params, Sw);
return VanGenuchten::krn(params, Swe);
}
private:
......
......@@ -166,7 +166,7 @@ public:
const Scalar x = 1 - std::pow(Swe, 1.0/params.vgM());
const Scalar xToM = std::pow(x, params.vgM());
return (1 - xToM)/std::sqrt(Sw) * ( (1 - xToM)/2 + 2*xToM*(1-x)/x );
return (1 - xToM)/std::sqrt(Swe) * ( (1 - xToM)/2 + 2*xToM*(1-x)/x );
};
......@@ -207,7 +207,7 @@ public:
const Scalar x = std::pow(Swe, 1.0/params.vgM());
return
-std::pow(1 - x, 2*params.vgM())
*std::pow(1 - Sw, -2/3)
*std::pow(1 - Swe, -2/3)
*(1.0/3 + 2*x/Swe);
}
......
......@@ -26,6 +26,7 @@
#include <dumux/common/exceptions.hh>
#include <dumux/common/propertysystem.hh>
#include <dune/common/timer.hh>
#include <dune/istl/istlexception.hh>
#include <iostream>
......@@ -118,6 +119,10 @@ protected:
SolutionVector &u = model().curSol();
JacobianAssembler &jacobianAsm = model().jacobianAssembler();
Dune::Timer assembleTimer(false);
Dune::Timer solveTimer(false);
Dune::Timer updateTimer(false);
// tell the controller that we begin solving
ctl.newtonBegin(*this, u);
......@@ -136,8 +141,11 @@ protected:
std::cout << "Assembling global jacobian";
std::cout.flush();
}
assembleTimer.start();
// linearize the problem at the current solution
jacobianAsm.assemble();
assembleTimer.stop();
// solve the resultuing linear equation system
if (ctl.verbose()) {
......@@ -150,17 +158,20 @@ protected:
std::cout.flush();
}
solveTimer.start();
// set the delta vector to zero before solving the linear system!
u = 0;
// ask the controller to solve the linearized system
ctl.newtonSolveLinear(jacobianAsm.matrix(),
u,
jacobianAsm.residual());
solveTimer.stop();
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_);
updateTimer.stop();
// tell the controller that we're done with this iteration
ctl.newtonEndStep(u, uOld_);
......@@ -169,6 +180,17 @@ protected:
// 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: "
<< 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";
}
if (!ctl.newtonConverged()) {
ctl.newtonFail();
return false;
......
......@@ -51,8 +51,7 @@ public:
static Dune::UGGrid<2> *create(const Dune::FieldVector<Scalar, 2> &upperRight,
const Dune::FieldVector<int, 2> &cellRes)
{
Dune::UGGrid<2> *grid = new Dune::UGGrid<2>;
Dune::GridFactory<Dune::UGGrid<2> > factory(grid);
Dune::GridFactory<Dune::UGGrid<2> > factory;
for (int i=0; i<=cellRes[0]; i++) {
for (int j=0; j<=cellRes[1]; j++) {
Dune::FieldVector<double,2> pos;
......@@ -95,7 +94,7 @@ public:
}
}
factory.createGrid();
Dune::UGGrid<2> *grid = factory.createGrid();
grid->loadBalance();
return grid;
}
......@@ -199,7 +198,6 @@ int main(int argc, char** argv)
//res[1] = 16;
std::auto_ptr<Grid> grid(CreateGrid<Grid, Scalar>::create(upperRight, res));
////////////////////////////////////////////////////////////
// instantiate and run the concrete problem
////////////////////////////////////////////////////////////
......
......@@ -51,7 +51,7 @@ NEW_TYPE_TAG(WaterAirProblem, INHERITS_FROM(BoxTwoPTwoCNI));
// Set the grid type
SET_PROP(WaterAirProblem, Grid)
{
typedef Dune::SGrid<2,2> type;
typedef Dune::YaspGrid<2> type;
};
SET_PROP(WaterAirProblem, LocalFEMSpace)
......@@ -82,6 +82,9 @@ SET_TYPE_PROP(WaterAirProblem,
// Enable gravity
SET_BOOL_PROP(WaterAirProblem, EnableGravity, true);
// Use forward differences instead of central differences
SET_INT_PROP(WaterAirProblem, NumericDifferenceMethod, +1);
// Write newton convergence
SET_BOOL_PROP(WaterAirProblem, NewtonWriteConvergence, false);
}
......
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