Commit 5fb3baff authored by Praetorius, Simon's avatar Praetorius, Simon
Browse files

some minor correction to eliminate warnings and pedantic errors

parent af1268e7
......@@ -23,12 +23,11 @@ QuasiCrystal_RB::QuasiCrystal_RB(const std::string &name_) :
Parameters::get(name + "->density", density);
Parameters::get(name + "->alpha", alpha);
minusC = -c;
};
}
void QuasiCrystal_RB::fillOperators()
{ FUNCNAME("QuasiCrystal_RB::fillOperators()");
{
// mass matrices
Operator *opM = new Operator(prob->getFeSpace(0), prob->getFeSpace(0));
opM->addTerm(new Simple_ZOT);
......@@ -150,11 +149,10 @@ void QuasiCrystal_RB::fillOperators()
prob->addVectorOperator(opM_4, 4);
prob->addMatrixOperator(opM, 4, 4, &minus1, &minus1); // theta
};
}
void QuasiCrystal_RB::fillBoundaryConditions()
{ FUNCNAME("QuasiCrystal_RB::fillBoundaryConditions()");
};
{ }
} }
\ No newline at end of file
} }
......@@ -12,7 +12,7 @@
#include "preconditioner/MTLPreconPfc.h"
#endif
#include "GenericOperatorTerm.h"
#include "Expressions.h"
using namespace AMDiS;
......@@ -25,39 +25,39 @@ public:
PfcPC(std::string name) : PhaseFieldCrystal(name) { }
/// initialize the preconditioners, i.e. set parameters
void initData()
void finalizeData() override
{
super::initData();
super::finalizeData();
// sequential PFC preconditioner
#if (defined HAVE_SEQ_PETSC) || (defined HAVE_PETSC)
PetscPreconPfc* runner = dynamic_cast<PetscPreconPfc*>(prob->getSolver()->getRunner());
if (runner) {
dynamic_cast<PetscSolver<PetscPreconPfc>*>(prob->getSolver())->setNested(true);
runner->setData(getTau());
PetscPreconPfc* precon = dynamic_cast<PetscPreconPfc*>(prob->getSolver()->getRightPrecon());
if (precon) {
precon->setData(getTau(), M0);
}
#endif
// PetscPreconPfcDiag* precon2 = dynamic_cast<PetscPreconPfcDiag*>(prob->getSolver()->getRightPrecon());
// if (precon2) {
// precon2->setData(getTau(), M0);
}
#elif !defined(HAVE_PARALLEL_DOMAIN_AMDIS)
MTLPreconPfc* precon = dynamic_cast<MTLPreconPfc*>(prob->getSolver()->getRightPrecon());
if (precon)
precon->setData(getTau(), M0);
#else
// parallel PFC preconditioner
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
Parallel::PetscSolverPfc* solver = dynamic_cast<Parallel::PetscSolverPfc*>(prob->getSolver());
if (solver)
solver->setData(getTau());
#endif
// sequential PFC preconditioner using MTL
#if (!defined HAVE_SEQ_PETSC) && (!defined HAVE_PETSC) && (!defined HAVE_PARALLEL_DOMAIN_AMDIS)
using AMDiS::extensions::MTLPreconPfc;
MTLPreconPfc* precon = dynamic_cast<MTLPreconPfc*>(prob->getSolver()->getRightPrecon());
if (precon)
precon->setData(getTau());
solver->setData(getTau(), M0);
// Parallel::PetscSolverPfcDiag* solver2 = dynamic_cast<Parallel::PetscSolverPfcDiag*>(prob->getSolver());
// if (solver2)
// solver2->setData(getTau(), M0);
#endif
}
/// generate initial solution for evolution equation
void solveInitialProblem(AdaptInfo *adaptInfo)
{ FUNCNAME("PFC_Demo::solveInitialProblem()");
void solveInitialProblem(AdaptInfo *adaptInfo) override
{
Flag initFlag = initDataFromFile(adaptInfo);
if (initFlag.isSet(DATA_ADOPTED))
return;
......@@ -79,16 +79,14 @@ int main(int argc, char** argv)
// add preconditioner / solver to the parameter list. Must be added before problem is initialized.
#if (defined HAVE_SEQ_PETSC) || (defined HAVE_PETSC)
CreatorMap<LinearSolver>::addCreator("petsc_pfc", new PetscSolver<PetscPreconPfc>::Creator);
#endif
#ifdef HAVE_PARALLEL_DOMAIN_AMDIS
CreatorMap<LinearSolver>::addCreator("p_petsc_pfc", new Parallel::PetscSolverPfc::Creator);
#endif
#if (!defined HAVE_SEQ_PETSC) && (!defined HAVE_PETSC) && (!defined HAVE_PARALLEL_DOMAIN_AMDIS)
using AMDiS::extensions::MTLPreconPfc;
CreatorMap<BasePreconditioner>::addCreator("pfc", new MTLPreconPfc::Creator);
CreatorMap<PetscPreconditionerNested>::addCreator("pfc", new PetscPreconPfc::Creator);
// CreatorMap<PetscPreconditionerNested>::addCreator("pfc_diag", new PetscPreconPfcDiag::Creator);
#elif !defined(HAVE_PARALLEL_DOMAIN_AMDIS)
CreatorMap<typename MTLPreconPfc::precon_base>::addCreator("pfc", new MTLPreconPfc::Creator);
// CreatorMap<typename MTLPreconPfc_Diag::base_precon>::addCreator("pfc", new MTLPreconPfc::Creator);
#else
CreatorMap<LinearSolverInterface>::addCreator("p_petsc_pfc", new Parallel::PetscSolverPfc::Creator);
// CreatorMap<LinearSolverInterface>::addCreator("p_petsc_pfc_diag", new Parallel::PetscSolverPfcDiag::Creator);
#endif
// create and initialize the PFC BaseProblem
......
......@@ -30,9 +30,9 @@ using namespace AMDiS::MTLTypes;
template<typename MatrixType>
struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
{
typedef BlockPreconditioner<MatrixType> super;
typedef ITL_BasePreconditioner<MatrixType, MTLTypes::MTLVector> precon_base;
typedef MTLPreconPfcBase<MatrixType> self;
typedef BlockPreconditioner<MatrixType> super;
typedef typename super::precon_base precon_base;
typedef MTLPreconPfcBase<MatrixType> self;
class Creator : public CreatorInterfaceName<precon_base>
{
......@@ -48,6 +48,8 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
: super(),
name(name),
tau(NULL),
M0(1.0),
matM(NULL), matL(NULL), matL0(NULL),
PId(NULL), PDiagM(NULL), PDiagMpL(NULL), PDiagMpL2(NULL),
solverM(NULL), solverMpL(NULL), solverMpL2(NULL)
{ }
......@@ -61,13 +63,22 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
void setData(Data* data)
{
tau = data->tau;
M0 = data->M0;
matM = &data->getM();
matL = &data->getL();
matL0 = &data->getL0();
}
void setData(double* tau_, double M0_ = 1.0)
{
tau = tau_;
M0 = M0_;
}
virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix);
virtual void init(const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override;
virtual void exit()
virtual void exit() override
{
if (PId) { delete PId; PId = NULL; }
if (PDiagM) { delete PDiagM; PDiagM = NULL; }
......@@ -78,7 +89,7 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
if (solverMpL2) { delete solverMpL2; solverMpL2 = NULL; }
}
virtual void solve(const MTLVector& b, MTLVector& x) const;
virtual void solve(const MTLVector& b, MTLVector& x) const override;
template <typename VectorX, typename VectorB>
void solveM(VectorX& x, const VectorB& b, int nIter = -1) const
......@@ -137,14 +148,21 @@ struct MTLPreconPfcBase : BlockPreconditioner<MatrixType>
}
}
const MTLMatrix& getM() const { return matM ? *matM : super::A->getSubMatrix(2,2); }
const MTLMatrix& getL0() const { return matL0 ? *matL0 : super::A->getSubMatrix(1,0); }
const MTLMatrix& getL() const { return matL ? *matL : super::A->getSubMatrix(2,1); }
const MTLMatrix& getM() const { return *matM; }
const MTLMatrix& getL() const { return *matL; }
double getTau() const { return *tau; }
protected:
std::string name;
double* tau;
double M0;
MTLMatrix* matM;
MTLMatrix* matL;
MTLMatrix* matL0;
MTLMatrix MpL;
MTLMatrix MpL2;
......@@ -161,10 +179,6 @@ protected:
mtl::matrix::umfpack::solver<MTLMatrix> *solverM;
mtl::matrix::umfpack::solver<MTLMatrix> *solverMpL;
mtl::matrix::umfpack::solver<MTLMatrix> *solverMpL2;
std::string name;
double* tau;
};
typedef MTLPreconPfcBase<MTLTypes::MTLMatrix> MTLPreconPfc;
......
......@@ -27,11 +27,11 @@ void MTLPreconPfcBase<MatrixType>::init(const SolverMatrix<Matrix<DOFMatrix*> >&
typedef typename mtl::Collection<MTLMatrix>::value_type value_type;
size_type n = num_rows(getM());
double delta = std::sqrt(getTau());
double delta = std::sqrt(M0 * getTau());
// helper-matrix MpL = M + delta*L
MpL.change_dim(n, n);
MpL = getM() + delta*getL();
MpL = getM() + (1.0/delta)*getL0();
// helper-matrix MpL = M + sqrt(delta)*L
MpL2.change_dim(n, n);
......@@ -81,11 +81,11 @@ void MTLPreconPfcBase<MatrixType>::solve(const MTLVector& b, MTLVector& x) const
const MTLVector b1(b[super::getRowRange(1)]);
const MTLVector b2(b[super::getRowRange(2)]);
double delta = std::sqrt(getTau());
double delta = std::sqrt(M0 * getTau());
solveM(y0, b0); // M*y0 = b0
y1 = getL() * y0; // y1 = K*y0
tmp = b1 - getTau()*y1; // tmp := b1 - tau*y1
y1 = getL0() * y0; // y1 = K*y0
tmp = b1 - y1; // tmp := b1 - tau*y1
solveMpL(y1, tmp); // (M + delta*K) * y1 = tmp
x0 = y0 + (1.0/delta)*y1; // x0 = y0 + (1/delta)*y1
......
......@@ -47,7 +47,7 @@ namespace AMDiS {
KSPSolve(data->kspM, b1, y1); // M*y1 = b1
MatMult(data->matK, y1, tmp); // tmp := K*y1
VecAYPX(tmp, -(data->tau), b2); // tmp := b2 - tau*tmp
VecAYPX(tmp, -sqr(data->delta), b2); // tmp := b2 - tau*tmp
KSPSolve(data->kspMpK, tmp, y2); // (M+sqrt(tau)K)*y2 = tmp
......@@ -83,7 +83,7 @@ namespace AMDiS {
PCShellSetContext(pc, &data);
double delta = sqrt(*tau);
double delta = std::sqrt((*tau) * M0);
MatNestGetSubMat(A.matrix, 2, 2, &data.matM);
MatNestGetSubMat(A.matrix, 2, 1, &data.matK);
......@@ -92,7 +92,7 @@ namespace AMDiS {
MatAXPY(MpK, delta, data.matK, SAME_NONZERO_PATTERN);
MatDuplicate(data.matM, MAT_COPY_VALUES, &MpK2);
MatAXPY(MpK2, sqrt(delta), data.matK, SAME_NONZERO_PATTERN);
MatAXPY(MpK2, std::sqrt(delta), data.matK, SAME_NONZERO_PATTERN);
// init sub-solvers
Runner::createSubSolver(data.kspM, data.matM, "mass_");
......
......@@ -49,12 +49,14 @@ namespace AMDiS {
public:
PetscPreconPfc()
: super("pfc_", "pfc"),
tau(NULL)
tau(NULL),
M0(1.0)
{ }
void setData(double *tauPtr)
void setData(double *tau_, double M0_ = 1.0)
{
tau = tauPtr;
tau = tau_;
M0 = M0_;
}
void init(PC pc, const SolverMatrix<Matrix<DOFMatrix*> >& A, const MatrixType& fullMatrix) override;
......@@ -66,6 +68,7 @@ namespace AMDiS {
PfcData data;
double *tau;
double M0;
};
}
......
......@@ -18,6 +18,7 @@
#include "PetscSolverPfc.h"
#include "parallel/PetscHelper.h"
#include "parallel/PetscSolverGlobalMatrix.h"
#include "Expressions.h"
namespace AMDiS { namespace Parallel {
......@@ -47,24 +48,24 @@ namespace AMDiS { namespace Parallel {
VecDuplicate(b1, &y2);
VecDuplicate(b1, &tmp);
KSPSolve(data->kspM, b1, y1); // M*y1 = b1
MatMult(data->matK, y1, tmp); // tmp := K*y1
VecAYPX(tmp, -(data->tau), b2); // tmp := b2 - tau*tmp
KSPSolve(data->kspM, b1, y1); // M*y1 = b1
MatMult(data->matK, y1, tmp); // tmp := K*y1
VecAYPX(tmp, -sqr(data->delta), b2); // tmp := b2 - M0*tau*tmp
KSPSolve(data->kspMpK, tmp, y2); // (M+sqrt(tau)K)*y2 = tmp
KSPSolve(data->kspMpK, tmp, y2); // (M+sqrt(tau)K)*y2 = tmp
MatMult(data->matM, y2, tmp); // tmp := M*y2
KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp
MatMult(data->matM, x2, tmp); // tmp := M*x2
KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp
MatMult(data->matM, y2, tmp); // tmp := M*y2
KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp
MatMult(data->matM, x2, tmp); // tmp := M*x2
KSPSolve(data->kspMpK2, tmp, x2); // MpK2*x2 = tmp
VecCopy(x2, x1); // x1 := x2
VecCopy(x2, x1); // x1 := x2
VecAXPBYPCZ(x1, 1.0, 1.0/(data->delta), -1.0/(data->delta), y1, y2); // x1 = 1*y1 + factor*y2 - factor*x1
MatMult(data->matK, x2, tmp); // tmp := K*x2
VecAYPX(tmp, -1.0, b3); // tmp := b3 - tmp
MatMult(data->matK, x2, tmp); // tmp := K*x2
VecAYPX(tmp, -1.0, b3); // tmp := b3 - tmp
KSPSolve(data->kspM, tmp, x3); // M*x3 = tmp
KSPSolve(data->kspM, tmp, x3); // M*x3 = tmp
VecDestroy(&y1);
VecDestroy(&y2);
......@@ -105,16 +106,15 @@ namespace AMDiS { namespace Parallel {
PCShellSetApply(getPc(), pcPfcShell);
PCShellSetContext(getPc(), &data);
double delta = sqrt(*tau);
double delta = std::sqrt((*tau) * M0);
const FiniteElemSpace *feSpace = componentSpaces[0];
// create mass-matrix
DOFMatrix matrixM(feSpace, feSpace);
Operator massOp(feSpace, feSpace);
Simple_ZOT zot;
massOp.addTerm(&zot);
matrixM.assembleOperator(massOp);
Operator opM(feSpace, feSpace);
addZOT(opM, 1.0);
matrixM.assembleOperator(opM);
solverM = createSubSolver(0, "M_");
solverM->fillPetscMatrix(&matrixM);
......@@ -123,11 +123,10 @@ namespace AMDiS { namespace Parallel {
// create MpK-matrix
DOFMatrix matrixMpK(feSpace, feSpace);
Operator laplaceOp2(feSpace, feSpace);
Simple_SOT sot2(delta);
laplaceOp2.addTerm(&zot);
laplaceOp2.addTerm(&sot2);
matrixMpK.assembleOperator(laplaceOp2);
Operator opMpK(feSpace, feSpace);
addZOT(opMpK, 1.0);
addSOT(opMpK, delta);
matrixMpK.assembleOperator(opMpK);
solverMpK = createSubSolver(0, "MpK_");
solverMpK->fillPetscMatrix(&matrixMpK);
......@@ -137,11 +136,10 @@ namespace AMDiS { namespace Parallel {
// create MpK2-matrix
DOFMatrix matrixMpK2(feSpace, feSpace);
Operator laplaceOp3(feSpace, feSpace);
Simple_SOT sot3(sqrt(delta));
laplaceOp3.addTerm(&zot);
laplaceOp3.addTerm(&sot3);
matrixMpK2.assembleOperator(laplaceOp3);
Operator opMpK2(feSpace, feSpace);
addZOT(opMpK2, 1.0);
addSOT(opMpK2, std::sqrt(delta));
matrixMpK2.assembleOperator(opMpK2);
solverMpK2 = createSubSolver(0, "MpK2_");
solverMpK2->fillPetscMatrix(&matrixMpK2);
......@@ -150,10 +148,9 @@ namespace AMDiS { namespace Parallel {
// create laplace-matrix
DOFMatrix matrixK(feSpace, feSpace);
Operator laplaceOp(feSpace, feSpace);
Simple_SOT sot;
laplaceOp.addTerm(&sot);
matrixK.assembleOperator(laplaceOp);
Operator opK(feSpace, feSpace);
addSOT(opK, 1.0);
matrixK.assembleOperator(opK);
solverK = createSubSolver(0, "K_");
solverK->fillPetscMatrix(&matrixK);
......
......@@ -48,14 +48,16 @@ namespace AMDiS { namespace Parallel {
solverMpK(NULL),
solverMpK2(NULL),
useOldInitialGuess(false),
tau(NULL)
tau(NULL),
M0(1.0)
{
Parameters::get(initFileStr + "->use old initial guess", useOldInitialGuess);
}
void setData(double *tauPtr)
void setData(double *tau_, double M0_ = 1.0)
{
tau = tauPtr;
tau = tau_;
M0 = M0_;
}
protected:
......@@ -73,10 +75,10 @@ namespace AMDiS { namespace Parallel {
PetscSolver *solverM, *solverK, *solverMpK, *solverMpK2;
bool useOldInitialGuess;
SystemVector* solution;
bool useOldInitialGuess;
double *tau;
double M0;
};
} }
......
Markdown is supported
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