Skip to content

Commit

Permalink
ROMInterfaceForm for interface EQP (#38)
Browse files Browse the repository at this point in the history 
* moved DGLaxFriedrichsFluxIntegrator to interface integrator.

* DGLaxFriedrichsFluxIntegrator interface gradient verified.

* AddInterfaceIntegrator typo fix.

* SampleInfo added itf.

* InterfaceIntegrator::AssembleQuadratureVector for interface

* ROMInterfaceForm::InterfaceAddMult tested.

* ROMInterfaceForm::InterfaceGetGradient. needs debugging.

* ROMInterfaceForm::InterfaceGetGradient verified.

* ROMInterfaceForm- eqp setup/training verified for a port.

* ROMInterfaceForm EQP training takes column indices of two snapshot matrices.
  • Loading branch information
@dreamer2368
dreamer2368 authored May 7, 2024
1 parent 4586b27 commit e30dfab
Show file tree
Hide file tree
Showing 22 changed files with 2,819 additions and 891 deletions.
8 changes: 8 additions & 0 deletions .github/workflows/ci.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,10 +53,18 @@ jobs:
run: |
cd ${GITHUB_WORKSPACE}/build/test
./nonlinear_integ_grad
- name: Test nonlinear interface integrator gradient
run: |
cd ${GITHUB_WORKSPACE}/build/test
./interfaceinteg_grad
- name: Test ROM NonlinearForm
run: |
cd ${GITHUB_WORKSPACE}/build/test
./test_rom_nonlinearform
- name: Test ROM InterfaceForm
run: |
cd ${GITHUB_WORKSPACE}/build/test
./test_rom_interfaceform
- name: Test Poisson DD solver
run: |
cd ${GITHUB_WORKSPACE}/build/test
Expand Down
4 changes: 2 additions & 2 deletions CMakeLists.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -191,8 +191,8 @@ set(scaleupROMObj_SOURCES
include/interface_form.hpp
src/interface_form.cpp

# include/multiblock_bilinearform.hpp
# src/multiblock_bilinearform.cpp
include/rom_interfaceform.hpp
src/rom_interfaceform.cpp

include/multiblock_solver.hpp
src/multiblock_solver.cpp
Expand Down
77 changes: 4 additions & 73 deletions include/hyperreduction_integ.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,12 +11,15 @@ namespace mfem
{

struct SampleInfo {
// TODO(kevin): all integrators could use the same integer variable, reducing the memory size of this struct.
int el; // element index (used for DomainIntegrator)
int face; // face index (used for InteriorFaceIntegrator)
int be; // boundary element index (used for BdrFaceIntegrator)
int itf; // interface info index (used for InterfaceIntegrator)
// can add dofs for other hyper reductions.

int qp; // quadrature point
double qw; // quadrature weight
// can add dofs for other hyper reductions.
};

class HyperReductionIntegrator : virtual public NonlinearFormIntegrator
Expand Down Expand Up @@ -224,78 +227,6 @@ class IncompressibleInviscidFluxNLFIntegrator :
const Vector &x, DenseMatrix &jac) override;
};

/*
Interior face, interface integrator
< [v], {uu \dot n} + \Lambda * [u] >
\Lambda = max( | u- \dot n |, | u+ \dot n | )
*/
class DGLaxFriedrichsFluxIntegrator : public HyperReductionIntegrator
{
private:
int dim, ndofs1, ndofs2, nvdofs;
double w, un1, un2, un;
Coefficient *Q{};

Vector nor, flux, shape1, shape2, u1, u2, tmp_vec;
DenseMatrix udof1, udof2, elv1, elv2;
DenseMatrix elmat_comp11, elmat_comp12, elmat_comp21, elmat_comp22, tmp;

// precomputed basis value at the sample point.
Array<DenseMatrix *> shapes1, shapes2;
public:
DGLaxFriedrichsFluxIntegrator(Coefficient &q)
: HyperReductionIntegrator(true), Q(&q) {}

void AssembleFaceVector(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &Tr,
const Vector &elfun, Vector &elvect) override;

/// @brief Assemble the local action of the gradient of the
/// NonlinearFormIntegrator resulting from a face integral term.
void AssembleFaceGrad(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &Tr,
const Vector &elfun, DenseMatrix &elmat) override;

void AssembleQuadratureVector(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &T,
const IntegrationPoint &ip,
const double &iw,
const Vector &eltest,
Vector &elquad) override;

void AssembleQuadratureGrad(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &T,
const IntegrationPoint &ip,
const double &iw,
const Vector &eltest,
DenseMatrix &quadmat) override;

void AppendPrecomputeInteriorFaceCoeffs(const FiniteElementSpace *fes,
DenseMatrix &basis,
const SampleInfo &sample) override;
void AppendPrecomputeBdrFaceCoeffs(const FiniteElementSpace *fes,
DenseMatrix &basis,
const SampleInfo &sample) override;

void AddAssembleVector_Fast(const int s, const double qw,
FaceElementTransformations &T, const IntegrationPoint &ip,
const Vector &x, Vector &y) override;
void AddAssembleGrad_Fast(const int s, const double qw,
FaceElementTransformations &T, const IntegrationPoint &ip,
const Vector &x, DenseMatrix &jac) override;

private:
void AppendPrecomputeFaceCoeffs(const FiniteElementSpace *fes,
FaceElementTransformations *T,
DenseMatrix &basis,
const SampleInfo &sample);

};

} // namespace mfem

#endif
31 changes: 18 additions & 13 deletions include/interface_form.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@ class InterfaceForm

/// Set of interior face Integrators to be assembled (added).
Array<InterfaceNonlinearFormIntegrator*> fnfi; // owned
// Array<Array<SampleInfo> *> fnfi_sample;

// For Mult and GetGradient.
mutable BlockVector x_tmp, y_tmp;
Expand All @@ -42,11 +41,13 @@ class InterfaceForm
NonlinearFormIntegrator%s and gradient Operator. */
virtual ~InterfaceForm();

/* access functions */
const Array<int>& GetBlockOffsets() { return block_offsets; }

/// Adds new Interior Face Integrator.
void AddIntefaceIntegrator(InterfaceNonlinearFormIntegrator *nlfi)
virtual void AddInterfaceIntegrator(InterfaceNonlinearFormIntegrator *nlfi)
{
fnfi.Append(nlfi);
// fnfi_sample.Append(NULL);
}

/** @brief Access all interface integrators added with
Expand All @@ -58,24 +59,28 @@ class InterfaceForm

void AssembleInterfaceMatrixAtPort(const int p, Array<FiniteElementSpace *> &fes_comp, Array2D<SparseMatrix *> &mats_p) const;

void InterfaceAddMult(const Vector &x, Vector &y) const;
virtual void InterfaceAddMult(const Vector &x, Vector &y) const;

void InterfaceGetGradient(const Vector &x, Array2D<SparseMatrix *> &mats) const;

protected:

// BilinearForm interface operator.
void AssembleInterfaceMatrix(Mesh *mesh1, Mesh *mesh2,
FiniteElementSpace *fes1, FiniteElementSpace *fes2,
Array<InterfaceInfo> *interface_infos, Array2D<SparseMatrix*> &mats) const;
virtual void InterfaceGetGradient(const Vector &x, Array2D<SparseMatrix *> &mats) const;

/*
this is public only for the sake of testing.
TODO(kevin): bring it back to protected.
*/
// NonlinearForm interface operator.
void AssembleInterfaceVector(Mesh *mesh1, Mesh *mesh2,
FiniteElementSpace *fes1, FiniteElementSpace *fes2,
Array<InterfaceInfo> *interface_infos,
const Vector &x1, const Vector &x2,
Vector &y1, Vector &y2) const;

protected:

// BilinearForm interface operator.
void AssembleInterfaceMatrix(Mesh *mesh1, Mesh *mesh2,
FiniteElementSpace *fes1, FiniteElementSpace *fes2,
Array<InterfaceInfo> *interface_infos, Array2D<SparseMatrix*> &mats) const;

void AssembleInterfaceGrad(Mesh *mesh1, Mesh *mesh2,
FiniteElementSpace *fes1, FiniteElementSpace *fes2,
Array<InterfaceInfo> *interface_infos,
Expand Down Expand Up @@ -116,7 +121,7 @@ class MixedInterfaceForm
virtual ~MixedInterfaceForm();

/// Adds new Interior Face Integrator.
void AddIntefaceIntegrator(InterfaceNonlinearFormIntegrator *nlfi)
void AddInterfaceIntegrator(InterfaceNonlinearFormIntegrator *nlfi)
{
fnfi.Append(nlfi);
// fnfi_sample.Append(NULL);
Expand Down
122 changes: 122 additions & 0 deletions include/interfaceinteg.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -52,6 +52,24 @@ class InterfaceNonlinearFormIntegrator : virtual public HyperReductionIntegrator
FaceElementTransformations &Trans2,
Array2D<DenseMatrix*> &elmats)
{ mfem_error("Abstract method InterfaceNonlinearFormIntegrator::AssembleInterfaceMatrix!\n"); }

virtual void AssembleQuadratureVector(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &Tr1,
FaceElementTransformations &Tr2,
const IntegrationPoint &ip,
const double &iw,
const Vector &eltest1, const Vector &eltest2,
Vector &elquad1, Vector &elquad2);

virtual void AssembleQuadratureGrad(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &Tr1,
FaceElementTransformations &Tr2,
const IntegrationPoint &ip,
const double &iw,
const Vector &eltest1, const Vector &eltest2,
Array2D<DenseMatrix*> &quadmats);
};

class InterfaceDGDiffusionIntegrator : public InterfaceNonlinearFormIntegrator
Expand Down Expand Up @@ -240,6 +258,110 @@ class InterfaceDGElasticityIntegrator : public InterfaceNonlinearFormIntegrator
using InterfaceNonlinearFormIntegrator::AssembleInterfaceMatrix;
};

/*
Interior face, interface integrator
< [v], {uu \dot n} + \Lambda * [u] >
\Lambda = max( | u- \dot n |, | u+ \dot n | )
*/
class DGLaxFriedrichsFluxIntegrator : public InterfaceNonlinearFormIntegrator
{
private:
int dim, ndofs1, ndofs2, nvdofs;
double w, un1, un2, un;
Coefficient *Q{};

Vector nor, flux, shape1, shape2, u1, u2, tmp_vec;
DenseMatrix udof1, udof2, elv1, elv2;
DenseMatrix elmat_comp11, elmat_comp12, elmat_comp21, elmat_comp22, tmp;

// precomputed basis value at the sample point.
Array<DenseMatrix *> shapes1, shapes2;
public:
DGLaxFriedrichsFluxIntegrator(Coefficient &q, const IntegrationRule *ir = NULL)
: InterfaceNonlinearFormIntegrator(true, ir), Q(&q) {}

void AssembleFaceVector(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &Tr,
const Vector &elfun, Vector &elvect) override;

/// @brief Assemble the local action of the gradient of the
/// NonlinearFormIntegrator resulting from a face integral term.
void AssembleFaceGrad(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &Tr,
const Vector &elfun, DenseMatrix &elmat) override;

void AssembleQuadratureVector(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &T,
const IntegrationPoint &ip,
const double &iw,
const Vector &eltest,
Vector &elquad) override;

void AssembleQuadratureGrad(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &T,
const IntegrationPoint &ip,
const double &iw,
const Vector &eltest,
DenseMatrix &quadmat) override;

void AppendPrecomputeInteriorFaceCoeffs(const FiniteElementSpace *fes,
DenseMatrix &basis,
const SampleInfo &sample) override;
void AppendPrecomputeBdrFaceCoeffs(const FiniteElementSpace *fes,
DenseMatrix &basis,
const SampleInfo &sample) override;

void AddAssembleVector_Fast(const int s, const double qw,
FaceElementTransformations &T, const IntegrationPoint &ip,
const Vector &x, Vector &y) override;
void AddAssembleGrad_Fast(const int s, const double qw,
FaceElementTransformations &T, const IntegrationPoint &ip,
const Vector &x, DenseMatrix &jac) override;

void AssembleInterfaceVector(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &Tr1,
FaceElementTransformations &Tr2,
const Vector &elfun1, const Vector &elfun2,
Vector &elvect1, Vector &elvect2) override;

void AssembleInterfaceGrad(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &Tr1,
FaceElementTransformations &Tr2,
const Vector &elfun1, const Vector &elfun2,
Array2D<DenseMatrix *> &elmats) override;

void AssembleQuadratureVector(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &Tr1,
FaceElementTransformations &Tr2,
const IntegrationPoint &ip,
const double &iw,
const Vector &eltest1, const Vector &eltest2,
Vector &elquad1, Vector &elquad2) override;

void AssembleQuadratureGrad(const FiniteElement &el1,
const FiniteElement &el2,
FaceElementTransformations &Tr1,
FaceElementTransformations &Tr2,
const IntegrationPoint &ip,
const double &iw,
const Vector &eltest1, const Vector &eltest2,
Array2D<DenseMatrix*> &quadmats) override;

private:
void AppendPrecomputeFaceCoeffs(const FiniteElementSpace *fes,
FaceElementTransformations *T,
DenseMatrix &basis,
const SampleInfo &sample);

};

} // namespace mfem

#endif
4 changes: 4 additions & 0 deletions include/linalg_utils.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -161,6 +161,10 @@ void AddSubMatrixRtAP(const DenseMatrix& R, const Array<int> &Rrows,
const DenseMatrix& A,
const DenseMatrix& P, const Array<int> &Prows,
DenseMatrix& RAP);
void AddSubMatrixRtAP(const DenseMatrix& R, const Array<int> &Rrows,
const DenseMatrix& A,
const DenseMatrix& P, const Array<int> &Prows,
SparseMatrix& RAP);

// DenseTensor is column major and i is the fastest index.
// y_k = T_{ijk} * x_i * x_j
Expand Down
Loading

0 comments on commit e30dfab

Please sign in to comment.