Skip to content

changed: cache L2 projection matrix #437

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
akva2 wants to merge 1 commit into
base: master
Choose a base branch
from
Open

changed: cache L2 projection matrix #437

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 9 additions & 0 deletions src/ASM/ASMbase.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
#include "ASM3D.h"
#include "IFEM.h"
#include "MPC.h"
#include "SparseMatrix.h"
#include "Tensor.h"
#include "Vec3.h"
#include "Vec3Oper.h"
Expand Down Expand Up @@ -63,6 +64,7 @@ ASMbase::ASMbase (unsigned char n_p, unsigned char n_s, unsigned char n_f)
nel = nnod = 0;
idx = 0;
firstIp = 0;
glbL2_A = nullptr;
}


Expand All @@ -80,6 +82,7 @@ ASMbase::ASMbase (const ASMbase& patch, unsigned char n_f)
nnod = patch.nnod;
idx = patch.idx;
firstIp = patch.firstIp;
glbL2_A = nullptr;
// Note: Properties are _not_ copied
}

Expand All @@ -96,6 +99,7 @@ ASMbase::ASMbase (const ASMbase& patch)
nnod = patch.nnod;
idx = patch.idx;
firstIp = patch.firstIp;
glbL2_A = nullptr;

// Only copy the regular part of the FE data, leave out any extraordinaries

Expand Down Expand Up @@ -127,6 +131,8 @@ ASMbase::~ASMbase ()
{
for (MPC* mpc : mpcs)
delete mpc;

delete glbL2_A;
}


Expand Down Expand Up @@ -174,6 +180,9 @@ void ASMbase::clear (bool retainGeometry)
BCode.clear();
dCode.clear();
mpcs.clear();

delete glbL2_A;
glbL2_A = nullptr;
}


Expand Down
2 changes: 2 additions & 0 deletions src/ASM/ASMbase.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -918,6 +918,8 @@ class ASMbase
static int gEl; //!< Global element counter
static int gNod; //!< Global node counter

mutable SparseMatrix* glbL2_A = nullptr; //!< Global L2-projection (mass) matrix

private:
std::vector<char> myLMTypes; //!< Type of Lagrange multiplier ('L' or 'G')
std::set<size_t> myLMs; //!< Nodal indices of the Lagrange multipliers
Expand Down
110 changes: 60 additions & 50 deletions src/ASM/GlbL2projector.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -166,7 +166,7 @@ public:
L2Mats(GlbL2& p, size_t nen, size_t nf, LocalIntegral* q = nullptr)
: gl2Int(p), elmData(q)
{
this->resize(1,nf);
this->resize(p.pA ? 1 : 0,nf);
this->redim(nen);
}

Expand All @@ -193,7 +193,8 @@ class L2GlobalInt : public GlobalIntegral
{
public:
//! \brief The constructor initializes the system matrix references.
L2GlobalInt(SparseMatrix& A_, StdVector& B_) : A(A_), B(B_) {}
L2GlobalInt(int dim_, SparseMatrix* A_, StdVector& B_)
: dim(dim_), A(A_), B(B_) {}

//! \brief Empty destructor.
virtual ~L2GlobalInt() {}
Expand All @@ -205,41 +206,47 @@ public:
for (size_t i = 0; i < elm->mnpc.size(); i++)
{
int inod = elm->mnpc[i]+1;
for (size_t j = 0; j < elm->mnpc.size(); j++)
{
int jnod = elm->mnpc[j]+1;
A(inod,jnod) += elm->A.front()(i+1,j+1);
if (A) {
for (size_t j = 0; j < elm->mnpc.size(); j++)
{
int jnod = elm->mnpc[j]+1;
(*A)(inod,jnod) += elm->A.front()(i+1,j+1);
}
}
for (const Vector& b : elm->b)
{
B(inod) += b[i];
inod += A.dim();
inod += dim;
}
}
return true;
}

private:
SparseMatrix& A; //!< Reference to left-hand-side matrix
int dim; //!< Dimension of matrix
SparseMatrix* A; //!< Reference to left-hand-side matrix
StdVector& B; //!< Reference to right-hand-side vector
};


GlbL2::GlbL2 (IntegrandBase* p, size_t n) : problem(p)
GlbL2::GlbL2 (IntegrandBase* p, size_t n, SparseMatrix* A) :
pA(A), problem(p)
{
nrhs = p->getNoFields(2);
this->allocate(n);
}


GlbL2::GlbL2 (FunctionBase* f, size_t n) : problem(nullptr), functions({f})
GlbL2::GlbL2 (FunctionBase* f, size_t n, SparseMatrix* A) :
pA(A), problem(nullptr), functions({f})
{
nrhs = f->dim();
this->allocate(n);
}


GlbL2::GlbL2 (const FunctionVec& f, size_t n) : problem(nullptr), functions(f)
GlbL2::GlbL2 (const FunctionVec& f, size_t n, SparseMatrix* A) :
pA(A), problem(nullptr), functions(f)
{
nrhs = 0;
for (FunctionBase* func : f)
Expand All @@ -250,32 +257,32 @@ GlbL2::GlbL2 (const FunctionVec& f, size_t n) : problem(nullptr), functions(f)

GlbL2::~GlbL2()
{
delete pA;
delete pB;
#ifdef HAS_PETSC
delete adm;
#endif
}


void GlbL2::allocate (size_t n)
{
#ifdef HAS_PETSC
adm = nullptr;
if (GlbL2::MatrixType == LinAlg::PETSC && GlbL2::SolverParams)
{
adm = new ProcessAdm();
pA = new PETScMatrix(*adm,*GlbL2::SolverParams);
pB = new PETScVector(*adm,n*nrhs);
static ProcessAdm adm;
if (!pA) {
pA = new PETScMatrix(adm,*GlbL2::SolverParams);
pA->redim(n,n);
}

pB = new PETScVector(adm,n*nrhs);
}
else
#endif
{
pA = new SparseMatrix(SparseMatrix::SUPERLU);
if (!pA) {
pA = new SparseMatrix(SparseMatrix::SUPERLU);
pA->redim(n,n);
}
pB = new StdVector(n*nrhs);
}

pA->redim(n,n);
}


Expand Down Expand Up @@ -353,7 +360,8 @@ bool GlbL2::evalInt (LocalIntegral& elmInt,
solPt.insert(solPt.end(),funcPt.begin(),funcPt.end());
}

gl2.A.front().outer_product(fe.N,fe.N,true,fe.detJxW);
if (!gl2.A.empty())
gl2.A.front().outer_product(fe.N,fe.N,true,fe.detJxW);
for (size_t j = 0; j < solPt.size(); j++)
gl2.b[j].add(fe.N,solPt[j]*fe.detJxW);

Expand All @@ -376,7 +384,8 @@ bool GlbL2::evalIntMx (LocalIntegral& elmInt,
if (!problem->diverged(fe.iGP+1))
return false;

gl2.A.front().outer_product(fe.N,fe.N,true,fe.detJxW);
if (!gl2.A.empty())
gl2.A.front().outer_product(fe.N,fe.N,true,fe.detJxW);
for (size_t j = 0; j < solPt.size(); j++)
gl2.b[j].add(fe.N,solPt[j]*fe.detJxW);

Expand Down Expand Up @@ -418,6 +427,7 @@ bool GlbL2::solve (Matrix& sField)
#if SP_DEBUG > 1
std::cout <<"\nSolution:"<< sField;
#endif

return true;
}

Expand Down Expand Up @@ -473,24 +483,24 @@ bool ASMbase::L2projection (Matrix& sField,
const TimeDomain& time)
{
PROFILE2("ASMbase::L2projection");

GlbL2 gl2(integrand,this->getNoNodes(1));
L2GlobalInt dummy(*gl2.pA,*gl2.pB);

gl2.preAssemble(MNPC,this->getNoElms(true));
GlbL2 gl2(integrand,this->getNoNodes(1),glbL2_A);
L2GlobalInt dummy(this->getNoNodes(1),glbL2_A ? nullptr : gl2.pA,*gl2.pB);
if (!glbL2_A)
gl2.preAssemble(MNPC,this->getNoElms(true));
glbL2_A = gl2.pA;
return this->integrate(gl2,dummy,time) && gl2.solve(sField);
}


bool ASMbase::L2projection (Matrix& sField, FunctionBase* function, double t)
{
PROFILE2("ASMbase::L2projection");

GlbL2 gl2(function,this->getNoNodes(1));
L2GlobalInt dummy(*gl2.pA,*gl2.pB);
GlbL2 gl2(function,this->getNoNodes(1),glbL2_A);
L2GlobalInt dummy(this->getNoNodes(1),glbL2_A ? nullptr : gl2.pA,*gl2.pB);
if (!glbL2_A)
gl2.preAssemble(MNPC,this->getNoElms(true));
glbL2_A = gl2.pA;
TimeDomain time; time.t = t;

gl2.preAssemble(MNPC,this->getNoElms(true));
return this->integrate(gl2,dummy,time) && gl2.solve(sField);
}

Expand All @@ -499,12 +509,12 @@ bool ASMbase::L2projection (const std::vector<Matrix*>& sField,
const FunctionVec& function, double t)
{
PROFILE2("ASMbase::L2projection");

GlbL2 gl2(function,this->getNoNodes(1));
L2GlobalInt dummy(*gl2.pA,*gl2.pB);
GlbL2 gl2(function,this->getNoNodes(1),glbL2_A);
L2GlobalInt dummy(this->getNoNodes(1),glbL2_A ? nullptr : gl2.pA,*gl2.pB);
if (!glbL2_A)
gl2.preAssemble(MNPC,this->getNoElms(true));
glbL2_A = gl2.pA;
TimeDomain time; time.t = t;

gl2.preAssemble(MNPC,this->getNoElms(true));
return this->integrate(gl2,dummy,time) && gl2.solve(sField);
}

Expand All @@ -520,44 +530,45 @@ bool ASMbase::globalL2projection (Matrix& sField,
// Assemble the projection matrices
size_t i, nnod = this->getNoProjectionNodes();
size_t j, ncomp = integrand.dim();
SparseMatrix* A;
StdVector* B;
switch (GlbL2::MatrixType) {
case LinAlg::UMFPACK:
A = new SparseMatrix(SparseMatrix::UMFPACK);
if (!glbL2_A)
glbL2_A = new SparseMatrix(SparseMatrix::UMFPACK);
B = new StdVector(nnod*ncomp);
break;
#ifdef HAS_PETSC
case LinAlg::PETSC:
if (GlbL2::SolverParams)
{
A = new PETScMatrix(ProcessAdm(), *GlbL2::SolverParams);
B = new PETScVector(ProcessAdm(), nnod*ncomp);
static ProcessAdm adm;
if (!glbL2_A)
glbL2_A = new PETScMatrix(adm, *GlbL2::SolverParams);
B = new PETScVector(adm, nnod*ncomp);
break;
}
#endif
default:
A = new SparseMatrix(SparseMatrix::SUPERLU);
glbL2_A = new SparseMatrix(SparseMatrix::SUPERLU);
B = new StdVector(nnod*ncomp);
}
A->redim(nnod,nnod);
glbL2_A->redim(nnod,nnod);

if (!this->assembleL2matrices(*A,*B,integrand,continuous))
if (!this->assembleL2matrices(*glbL2_A,*B,integrand,continuous))
{
delete A;
delete B;
return false;
}

#if SP_DEBUG > 1
std::cout <<"---- Matrix A -----\n"<< *A
std::cout <<"---- Matrix A -----\n"<< *glbL2_A
<<"-------------------"<< std::endl;
std::cout <<"---- Vector B -----"<< *B
<<"-------------------"<< std::endl;
#endif

// Solve the patch-global equation system
if (!A->solve(*B)) return false;
if (!glbL2_A->solve(*B)) return false;

// Store the control-point values of the projected field
sField.resize(ncomp,nnod);
Expand All @@ -569,7 +580,6 @@ bool ASMbase::globalL2projection (Matrix& sField,
std::cout <<"- Solution Vector -"<< sField
<<"-------------------"<< std::endl;
#endif
delete A;
delete B;
return true;
}
13 changes: 6 additions & 7 deletions src/ASM/GlbL2projector.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,7 +24,6 @@ class FunctionBase;
class LinSolParams;
class SparseMatrix;
class StdVector;
class ProcessAdm;

typedef std::vector<int> IntVec; //!< Vector of integers
typedef std::vector<size_t> uIntVec; //!< Vector of unsigned integers
Expand Down Expand Up @@ -113,15 +112,18 @@ class GlbL2 : public Integrand
//! \brief The constructor initializes the projection matrices.
//! \param[in] p The main problem integrand
//! \param[in] n Dimension of the L2-projection matrices (number of nodes)
GlbL2(IntegrandBase* p, size_t n);
//! \param[in] A Sparse matrix to use
GlbL2(IntegrandBase* p, size_t n, SparseMatrix* A);
//! \brief Alternative constructor for projection of an explicit function.
//! \param[in] f The function to do L2-projection on
//! \param[in] n Dimension of the L2-projection matrices (number of nodes)
GlbL2(FunctionBase* f, size_t n);
//! \param[in] A Sparse matrix to use
GlbL2(FunctionBase* f, size_t n, SparseMatrix* A);
//! \brief Alternative constructor for projection of explicit functions.
//! \param[in] f The functions to do L2-projection on
//! \param[in] n Dimension of the L2-projection matrices (number of nodes)
GlbL2(const FunctionVec& f, size_t n);
//! \param[in] A Sparse matrix to use
GlbL2(const FunctionVec& f, size_t n, SparseMatrix* A);
//! \brief The destructor frees the system matrix and system vector.
virtual ~GlbL2();

Expand Down Expand Up @@ -203,9 +205,6 @@ class GlbL2 : public Integrand
IntegrandBase* problem; //!< The main problem integrand
FunctionVec functions; //!< Explicit functions to L2-project
size_t nrhs; //!< Number of right-hand-size vectors
#ifdef HAS_PETSC
ProcessAdm* adm; //!< Process administrator for PETSc
#endif
};

#endif
Loading