Self defined PetscPreconditioner #1134
Comments
|
@stefanozampini or @jandrej, do you think you can help? |
|
What exactly is not working? Can you provide the error message or provide a minimal failing example so that I can run and test the code? |
|
I think what you need is a shell preconditioner. See https://github.com/mfem/mfem/blob/tds-mhd-dev/miniapps/mhd/imResistiveMHDOperatorp.hpp Take a look at how MyBlockSolver is defined and used. You can take a look at imMHDp and play with it. Here is the instruction https://github.com/mfem/mfem/blob/tds-mhd-dev/miniapps/mhd/README.md |
What kind of errors? compilation errors? or issues at runtime? Or what else? |
|
Thank you for your reply @jandrej , @tangqi , @stefanozampini . Yes, @tangqi you are right I need a shell preconditioner. I will dig into Sorry about MWC @jandrej and @stefanozampini. Below is an example code which is obtained from by simplifying #include "mfem.hpp"
#include <iostream>
using namespace std;
using namespace mfem;
class SelfDefined_Preconditioner: public PetscPreconditioner
{
private:
Operator &op1, &op2;
PetscLinearSolver *ksp1, *ksp2;
Array<int>& block_offsets;
public:
SelfDefined_Preconditioner(MPI_Comm comm, Operator& op1_, Operator& op2_,
const string& prefix, Array<int>& offsets)
: PetscPreconditioner(comm, prefix), op1(op1_), op2(op2_), block_offsets(offsets)
{
height = op2.Height() + op2.Height();
width = op1.Width() + op2.Width();
ksp1 = new PetscLinearSolver(static_cast<PetscParMatrix&> (op1), "prec1_");
ksp2 = new PetscLinearSolver(static_cast<PetscParMatrix&> (op2), "prec2_");
}
virtual ~SelfDefined_Preconditioner() {}
virtual void SetOperator(const Operator& op) {} // avoid calling base class SetOperator()
virtual void Mult(const Vector& b, Vector& x) const
{
// here divide Vector b into 2 sub Vector b1, b2, but how to do it?
// Vector b1, b2;
// Vector x1, x2;
// ksp1->Mult(b1, x1);
// ksp2->Mult(b2, x2);
// here merge 2 sub Vector x1, x2 into x, but how to do it?
x = b; // just for avoiding compiling failure
}
};
int main(int argc, char *argv[])
{
int num_procs, myid;
MPI_Init(&argc, &argv);
MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
MPI_Comm_rank(MPI_COMM_WORLD, &myid);
const char *petscrc_file = "./petsc_opts.txt";
MFEMInitializePetsc(NULL, NULL, petscrc_file, NULL);
Mesh *mesh = new Mesh(20, 20, Element::TRIANGLE, true, 1.0, 1.0);
int dim = mesh->Dimension();
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
delete mesh;
FiniteElementCollection *hdiv_coll(new RT_FECollection(1, dim));
FiniteElementCollection *l2_coll(new L2_FECollection(1, dim));
ParFiniteElementSpace *R_space = new ParFiniteElementSpace(pmesh, hdiv_coll);
ParFiniteElementSpace *W_space = new ParFiniteElementSpace(pmesh, l2_coll);
Array<int> block_offsets(3);
block_offsets[0] = 0;
block_offsets[1] = R_space->GetVSize();
block_offsets[2] = W_space->GetVSize();
block_offsets.PartialSum();
Array<int> block_trueOffsets(3);
block_trueOffsets[0] = 0;
block_trueOffsets[1] = R_space->TrueVSize();
block_trueOffsets[2] = W_space->TrueVSize();
block_trueOffsets.PartialSum();
ConstantCoefficient one(1.0);
Vector two_(dim);
two_ = 2.0;
VectorConstantCoefficient two(two_);
// This example code solves a simple 2D mixed Darcy problem
// k*u + grad p = f
// -div u = g
// Assemble the finite element matrices for the Darcy operator
// D = [ M B^T ]
// [ B 0 ]
// where:
// M = \int_\Omega k u_h \cdot v_h d\Omega u_h, v_h \in R_h
// B = -\int_\Omega \div u_h q_h d\Omega u_h \in R_h, q_h \in W_h
//
// Construct the operators for preconditioner
// P = [ M 0 ]
// [ 0 S ],
// where:
// S = \int_\Omega p_h q_h d\Omega p_h, q_h \in W_h
//
ParLinearForm* fForm(new ParLinearForm);
ParLinearForm* gForm(new ParLinearForm);
ParBilinearForm* mVarf(new ParBilinearForm(R_space));
ParMixedBilinearForm* bVarf(new ParMixedBilinearForm(R_space, W_space));
ParBilinearForm* sVarf(new ParBilinearForm(W_space));
BlockVector x(block_offsets), rhs(block_offsets);
BlockVector trueX(block_trueOffsets), trueRhs(block_trueOffsets);
PetscParMatrix *pM = NULL, *pB = NULL, *pBT = NULL, *pS = NULL;
Operator::Type tid = Operator::PETSC_MATAIJ;
OperatorHandle Mh(tid), Bh(tid), Sh(tid);
fForm->Update(R_space, rhs.GetBlock(0), 0);
fForm->AddDomainIntegrator(new VectorFEDomainLFIntegrator(two));
fForm->Assemble();
fForm->ParallelAssemble(trueRhs.GetBlock(0));
gForm->Update(W_space, rhs.GetBlock(1), 0);
gForm->AddDomainIntegrator(new DomainLFIntegrator(one));
gForm->Assemble();
gForm->ParallelAssemble(trueRhs.GetBlock(1));
mVarf->AddDomainIntegrator(new VectorFEMassIntegrator);
mVarf->Assemble();
mVarf->Finalize();
mVarf->ParallelAssemble(Mh);
Mh.Get(pM);
Mh.SetOperatorOwner(false);
bVarf->AddDomainIntegrator(new VectorFEDivergenceIntegrator);
bVarf->Assemble();
bVarf->Finalize();
bVarf->ParallelAssemble(Bh);
Bh.Get(pB);
Bh.SetOperatorOwner(false);
sVarf->AddDomainIntegrator(new MassIntegrator);
sVarf->Assemble();
sVarf->Finalize();
sVarf->ParallelAssemble(Sh);
Sh.Get(pS);
Sh.SetOperatorOwner(false);
(*pB) *= -1;
pBT = pB->Transpose();
Operator *darcyOperator = NULL;
{
BlockOperator *tdarcyOp = new BlockOperator(block_trueOffsets);
tdarcyOp->SetBlock(0,0,pM);
tdarcyOp->SetBlock(0,1,pBT);
tdarcyOp->SetBlock(1,0,pB);
darcyOperator = new PetscParMatrix(pM->GetComm(),tdarcyOp, Operator::PETSC_MATAIJ);
delete tdarcyOp;
}
SelfDefined_Preconditioner darcyPreconditioner(MPI_COMM_WORLD, *pM, *pS, "prec_", block_offsets);
trueX = 0.0;
{
PetscLinearSolver *solver = new PetscLinearSolver(MPI_COMM_WORLD, "main_");
solver->SetOperator(*darcyOperator);
solver->SetPreconditioner(darcyPreconditioner);
solver->SetPrintLevel(2);
solver->Mult(trueRhs, trueX);
delete solver;
}
// 17. Free the used memory.
delete fForm;
delete gForm;
delete darcyOperator;
delete pBT;
delete pB;
delete pM;
delete mVarf;
delete bVarf;
delete sVarf;
delete W_space;
delete R_space;
delete l2_coll;
delete hdiv_coll;
delete pmesh;
// We finalize PETSc
MFEMFinalizePetsc();
MPI_Finalize();
return 0;
}
|
|
I compiled your code but I don't know what you think is the issue . |
|
I run your example and got the error below, This is because PetscPreconditioner has a
|
|
Thanks @stefanozampini . It works. I should be careful about protected members when defining myself class next time. |
Hi all,
I want to define a very general petsc user-defined preconditioner, which is just like this one. And I know I just need to solve several separate linear systems in sequence like this explanation. I also read about using Operator, but it's different if I use PetscParMatrix.
The below
class SelfDefined_Preconditionerhas some errors but I don't where and how to fix it.Thanks for any help in advance!
The text was updated successfully, but these errors were encountered: