Commit e667ab13 authored by Andreas Lauser's avatar Andreas Lauser
Browse files

box models: make them work without PDELab

since parallel the linear solver is hard-coded into PDElab, this only
works for sequential runs. to try it out, change the value of
HAVE_DUNE_PDELAB to 0 in the config.h file and add

...

around the definition of the LocalFEMSpace property.

git-svn-id: svn://svn.iws.uni-stuttgart.de/DUMUX/dumux/trunk@4826 2fb0f335-1f38-0410-981e-8018bf24f1b0
parent 95bab6dd
......@@ -76,6 +76,7 @@ NEW_PROP_TAG(TimeManager); //!< Manages the simulation time
NEW_PROP_TAG(NewtonMethod); //!< The type of the newton method
NEW_PROP_TAG(NewtonController); //!< The type of the newton controller
#if HAVE_DUNE_PDELAB
// properties for the PDELab wrapper
NEW_PROP_TAG(LocalFEMSpace); //!< The local finite element space used for the finite element interpolation
NEW_PROP_TAG(ScalarGridFunctionSpace); //!< The used grid function space for a single finite element function
......@@ -84,6 +85,7 @@ NEW_PROP_TAG(Constraints); //!< The constraints on the grid function space
NEW_PROP_TAG(ConstraintsTrafo); //!< The type of PDELab's constraints transformation
NEW_PROP_TAG(LocalOperator); //!< The type of the local operator used by PDELab
NEW_PROP_TAG(GridOperatorSpace); //!< The used grid operator space
#endif
//! Specify whether the jacobian matrix of the last iteration of a
//! time step should be re-used as the jacobian of the first iteration
......
......@@ -126,6 +126,7 @@ SET_TYPE_PROP(BoxModel, LocalJacobian, Dumux::BoxLocalJacobian<TypeTag>);
/*!
* \brief The type of a solution for the whole grid at a fixed time.
*/
#if HAVE_DUNE_PDELAB
SET_PROP(BoxModel, SolutionVector)
{ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
......@@ -133,6 +134,15 @@ SET_PROP(BoxModel, SolutionVector)
public:
typedef typename GridFunctionSpace::template VectorContainer<Scalar>::Type type;
};
#else
SET_PROP(BoxModel, SolutionVector)
{ private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
public:
typedef Dune::BlockVector<Dune::FieldVector<Scalar, numEq> > type;
};
#endif
/*!
* \brief The type of a solution for a whole element.
......@@ -177,6 +187,33 @@ public:
*/
SET_TYPE_PROP(BoxModel, JacobianAssembler, Dumux::PDELab::BoxAssembler<TypeTag>);
//! use central differences to calculate the jacobian by default
SET_INT_PROP(BoxModel, NumericDifferenceMethod, 0);
// disable jacobian matrix recycling by default
SET_BOOL_PROP(BoxModel, EnableJacobianRecycling, false);
// disable partial reassembling by default
SET_BOOL_PROP(BoxModel, EnablePartialReassemble, false);
// disable time-step ramp up by default
SET_BOOL_PROP(BoxModel, EnableTimeStepRampUp, false);
#if ! HAVE_DUNE_PDELAB
//! Set the type of a global jacobian matrix from the solution types
SET_PROP(BoxModel, JacobianMatrix)
{
private:
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
enum { numEq = GET_PROP_VALUE(TypeTag, PTAG(NumEq)) };
typedef typename Dune::FieldMatrix<Scalar, numEq, numEq> MatrixBlock;
public:
typedef typename Dune::BCRSMatrix<MatrixBlock> type;
};
#endif
#if HAVE_DUNE_PDELAB
//! Extract the type of a global jacobian matrix from the solution types
SET_PROP(BoxModel, JacobianMatrix)
{
......@@ -248,15 +285,7 @@ public:
SET_PROP(BoxModel, LocalOperator)
{ typedef typename GET_PROP(TypeTag, PTAG(GridOperatorSpace))::LocalOperator type; };
//! use central differences to calculate the jacobian by default
SET_INT_PROP(BoxModel, NumericDifferenceMethod, 0);
// disable jacobian matrix recycling by default
SET_BOOL_PROP(BoxModel, EnableJacobianRecycling, false);
// disable partial reassembling by default
SET_BOOL_PROP(BoxModel, EnablePartialReassemble, false);
// disable time-step ramp up by default
SET_BOOL_PROP(BoxModel, EnableTimeStepRampUp, false);
#endif // HAVE_DUNE_PDELAB
} // namespace Properties
} // namespace Dumux
......
......@@ -37,17 +37,20 @@ class BoxAssembler
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Problem)) Problem;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridView)) GridView;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
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;
#if HAVE_DUNE_PDELAB
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalFEMSpace)) FEM;
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;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(ConstraintsTrafo)) ConstraintsTrafo;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(LocalOperator)) LocalOperator;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(GridOperatorSpace)) GridOperatorSpace;
#endif
typedef typename GET_PROP_TYPE(TypeTag, PTAG(SolutionVector)) SolutionVector;
typedef typename GET_PROP_TYPE(TypeTag, PTAG(JacobianMatrix)) JacobianMatrix;
......@@ -113,7 +116,7 @@ public:
BoxAssembler()
{
problemPtr_ = 0;
#if HAVE_DUNE_PDELAB
fem_ = 0;
cn_ = 0;
scalarGridFunctionSpace_ = 0;
......@@ -121,6 +124,9 @@ public:
constraintsTrafo_ = 0;
localOperator_ = 0;
gridOperatorSpace_ = 0;
#endif // HAVE_DUNE_PDELAB
problemPtr_ = 0;
matrix_ = 0;
// set reassemble tolerance to 0, so that if partial
......@@ -133,6 +139,8 @@ public:
~BoxAssembler()
{
delete matrix_;
#if HAVE_DUNE_PDELAB
delete gridOperatorSpace_;
delete localOperator_;
delete constraintsTrafo_;
......@@ -140,6 +148,7 @@ public:
delete scalarGridFunctionSpace_;
delete cn_;
delete fem_;
#endif
}
/*!
......@@ -154,6 +163,11 @@ public:
void init(Problem& problem)
{
problemPtr_ = &problem;
#if !HAVE_DUNE_PDELAB
// initialize the BCRS matrix
createMatrix_();
#else
fem_ = new FEM();
//cn_ = new Constraints(*problemPtr_);
cn_ = new Constraints();
......@@ -171,10 +185,11 @@ public:
gridOperatorSpace_ =
new GridOperatorSpace(*gridFunctionSpace_, *constraintsTrafo_,
*gridFunctionSpace_, *constraintsTrafo_, *localOperator_);
matrix_ = new Matrix(*gridOperatorSpace_);
#endif
// initialize the jacobian matrix and the right hand side
// vector
matrix_ = new Matrix(*gridOperatorSpace_);
*matrix_ = 0;
reuseMatrix_ = false;
......@@ -496,6 +511,7 @@ public:
return elementColor_[globalElementIdx];
}
#if HAVE_DUNE_PDELAB
/*!
* \brief Returns a pointer to the PDELab's grid function space.
*/
......@@ -512,6 +528,7 @@ public:
{
return *constraintsTrafo_;
}
#endif // HAVE_DUNE_PDELAB
/*!
* \brief Return constant reference to global Jacobian matrix.
......@@ -527,6 +544,54 @@ public:
private:
#if !HAVE_DUNE_PDELAB
// Construct the BCRS matrix for the global jacobian
void createMatrix_()
{
int nVerts = gridView_().size(dim);
// allocate raw matrix
matrix_ = new Matrix(nVerts, nVerts, Matrix::random);
// find out how many neighbors each vertex has
typedef std::set<int> NeighborSet;
std::vector<NeighborSet > neighbors(nVerts);
ElementIterator eIt = gridView_().template begin<0>();
const ElementIterator eEndIt = gridView_().template end<0>();
for (; eIt != eEndIt; ++eIt) {
const Element &elem = *eIt;
// loop over all element vertices
int n = elem.template count<dim>();
for (int i = 0; i < n; ++i) {
int globalI = vertexMapper_().map(*eIt, i, dim);
for (int j = 0; j < n; ++j) {
int globalJ = vertexMapper_().map(*eIt, j, dim);
// insert into BCRS matrix
neighbors[globalI].insert(globalJ);
}
}
};
// allocate space for the rows
for (int i = 0; i < nVerts; ++i) {
matrix_->setrowsize(i, neighbors[i].size());
}
matrix_->endrowsizes();
// allocate space for the rows
for (int i = 0; i < nVerts; ++i) {
// off-diagonal entries
typename NeighborSet::iterator nIt = neighbors[i].begin();
typename NeighborSet::iterator nEndIt = neighbors[i].end();
for (; nIt != nEndIt; ++nIt) {
matrix_->addindex(i, *nIt);
}
}
matrix_->endindices();
};
#endif
// reset the global linear system of equations. if partial
// reassemble is enabled, this means that the jacobian matrix must
// only be erased partially!
......@@ -647,6 +712,7 @@ private:
Scalar nextReassembleTolerance_;
Scalar reassembleTolerance_;
#if HAVE_DUNE_PDELAB
// PDELab stuff
Constraints *cn_;
FEM *fem_;
......@@ -655,6 +721,7 @@ private:
ConstraintsTrafo *constraintsTrafo_;
LocalOperator *localOperator_;
GridOperatorSpace *gridOperatorSpace_;
#endif
};
} // namespace PDELab
......
......@@ -685,6 +685,28 @@ protected:
if (gridView_().comm().rank() != 0)
verbosity = 0;
#if ! HAVE_DUNE_PDELAB
typedef Dune::SeqILU0<JacobianMatrix, Vector, Vector> Preconditioner;
Preconditioner precond(A, 1.0);
typedef Dune::MatrixAdapter<JacobianMatrix,Vector,Vector> MatrixAdapter;
MatrixAdapter operatorA(A);
typedef Dune::BiCGSTABSolver<Vector> Solver;
Solver solver(operatorA, precond, residReduction, 500, verbosity);
// typedef Dune::RestartedGMResSolver<Vector> Solver;
// Solver solver(operatorA, precond, residReduction, 50, 500, verbosity);
Dune::InverseOperatorResult result;
Vector bTmp(b);
solver.apply(x, bTmp, result);
if (!result.converged)
DUNE_THROW(Dumux::NumericalProblem,
"Solving the linear system of equations did not converge.");
#else // HAVE_DUNE_PDELAB
#if HAVE_PARDISO
typedef Dumux::PDELab::ISTLBackend_NoOverlap_Loop_Pardiso<TypeTag> Solver;
Solver solver(problem_(), 500, verbosity);
......@@ -707,6 +729,7 @@ protected:
if (!solver.result().converged)
DUNE_THROW(Dumux::NumericalProblem,
"Solving the linear system of equations did not converge.");
#endif // HAVE_DUNE_PDELAB
// make sure the solver didn't produce a nan or an inf
// somewhere. this should never happen but for some strange
......
......@@ -64,6 +64,7 @@ SET_PROP(LensProblem, Grid)
#endif
};
#if HAVE_DUNE_PDELAB
SET_PROP(LensProblem, LocalFEMSpace)
{
typedef typename GET_PROP_TYPE(TypeTag, PTAG(Scalar)) Scalar;
......@@ -77,6 +78,7 @@ public:
typedef Dune::PDELab::P1LocalFiniteElementMap<Scalar,Scalar,dim> type; // for simplices
#endif
};
#endif
// Set the problem property
SET_PROP(LensProblem, 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