Skip to content

Commit

Permalink
EQP ROM workflow for unsteady NS (#50)
Browse files Browse the repository at this point in the history 
* ROMHandler: support block ordering by variable.

* MultiBlockSolver::AssembleROMMat kernel

* UnsteadyNSSolver: assemble ROM Mat.

* UnsteadyNSSolver::SolveROM verified.

* minor fix

* UnsteadyNSSolver::SetupInitialCondition

* minor fix

* SampleGenerator::SaveSnapshot- indexing bug hot fix

* fix modulo by zero bug..
  • Loading branch information
@dreamer2368
dreamer2368 authored Jul 9, 2024
1 parent 87185c7 commit 97b5355
Show file tree
Hide file tree
Showing 17 changed files with 657 additions and 158 deletions.
3 changes: 3 additions & 0 deletions include/multiblock_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -259,6 +259,9 @@ friend class ParameterizedProblem;
void ComputeSubdomainErrorAndNorm(GridFunction *fom_sol, GridFunction *rom_sol, double &error, double &norm);
void ComputeRelativeError(Array<GridFunction *> fom_sols, Array<GridFunction *> rom_sols, Vector &error);
void CompareSolution(BlockVector &test_U, Vector &error);

protected:
virtual void AssembleROMMat(BlockMatrix &romMat);
};

#endif
2 changes: 2 additions & 0 deletions include/rom_element_collection.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -50,6 +50,8 @@ class ROMLinearElement : public ROMElementCollection
Array<Array<MatrixBlocks *> *> bdr;
Array<MatrixBlocks *> port; // reference ports.

Array<MatrixBlocks *> mass; // Size(num_components);

public:
ROMLinearElement(TopologyHandler *topol_handler_,
const Array<FiniteElementSpace *> &fes_,
Expand Down
48 changes: 29 additions & 19 deletions include/rom_handler.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -31,6 +31,12 @@ enum NonlinearHandling
NUM_NLNHNDL
};

enum class ROMOrderBy
{
VARIABLE,
DOMAIN
};

const BasisTag GetBasisTagForComponent(const int &comp_idx, const TopologyHandler *topol_handler, const std::string var_name="");
const BasisTag GetBasisTag(const int &subdomain_index, const TopologyHandler *topol_handler, const std::string var_name="");

Expand Down Expand Up @@ -90,16 +96,11 @@ class ROMHandlerBase
/*
offset for the global domain ROM blocks.
For i-th subdomain and j-th variable,
index = i * num_var + j
ROMOrderBy::DOMAIN: index = i * num_var + j
ROMOrderBy::VARIABLE: index = j * numSub + i
*/
Array<int> num_basis;
Array<int> rom_block_offsets;
/*
the global domain ROM blocks with index order switched.
For i-th subdomain and j-th variable,
index = j * num_var + i
*/
Array<int> rom_varblock_offsets;

CAROM::Options* rom_options;
CAROM::BasisGenerator *basis_generator;
Expand All @@ -109,6 +110,11 @@ class ROMHandlerBase
bool update_right_SV = false;
bool incremental = false;

ROMOrderBy ordering = ROMOrderBy::DOMAIN;

mfem::BlockVector *reduced_rhs = NULL;
mfem::BlockVector *reduced_sol = NULL;

void ParseInputs();
public:
ROMHandlerBase(TopologyHandler *input_topol, const Array<int> &input_var_offsets,
Expand All @@ -129,7 +135,6 @@ class ROMHandlerBase
const std::string GetBasisPrefix() { return basis_prefix; }
const BasisTag GetRefBasisTag(const int ref_idx) { return basis_tags[ref_idx]; }
const Array<int>* GetBlockOffsets() { return &rom_block_offsets; }
const Array<int>* GetVarBlockOffsets() { return &rom_varblock_offsets; }
virtual SparseMatrix* GetOperator() = 0;
const bool GetNonlinearMode() { return nonlinear_mode; }
void SetNonlinearMode(const bool nl_mode)
Expand All @@ -139,6 +144,12 @@ class ROMHandlerBase
nonlinear_mode = nl_mode;
}
const NonlinearHandling GetNonlinearHandling() { return nlin_handle; }
const ROMOrderBy GetOrdering() { return ordering; }
const int GetBlockIndex(const int m, const int v=-1);
void GetDomainAndVariableIndex(const int &rom_block_index, int &m, int &v);

mfem::BlockVector* GetReducedSolution() { return reduced_sol; }
mfem::BlockVector* GetReducedRHS() { return reduced_rhs; }

/* parse inputs for supremizer. only for Stokes/SteadyNS Solver. */
void ParseSupremizerInput(Array<int> &num_ref_supreme, Array<int> &num_supreme);
Expand All @@ -157,9 +168,9 @@ class ROMHandlerBase
virtual void ProjectToDomainBasis(const Array<int> &idx_i, const Array<int> &idx_j, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;

virtual void ProjectComponentToRefBasis(const int &c, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
virtual void ProjectComponentToDomainBasis(const int &c, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
virtual void ProjectComponentToDomainBasis(const int &m, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
virtual void ProjectInterfaceToRefBasis(const int &c1, const int &c2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
virtual void ProjectInterfaceToDomainBasis(const int &c1, const int &c2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
virtual void ProjectInterfaceToDomainBasis(const int &m1, const int &m2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
virtual void ProjectVariableToDomainBasis(const int &vi, const int &vj, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
virtual void ProjectGlobalToDomainBasis(const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats) = 0;
virtual void ProjectOperatorOnReducedBasis(const Array2D<Operator*> &mats) = 0;
Expand All @@ -173,12 +184,13 @@ class ROMHandlerBase
virtual void LiftUpFromDomainBasis(const int &i, const Vector &rom_vec, Vector &vec) = 0;
virtual void LiftUpGlobal(const BlockVector &rom_vec, BlockVector &vec) = 0;

virtual void Solve(BlockVector &rhs, BlockVector &sol) = 0;
virtual void Solve(BlockVector* U) = 0;
virtual void NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec=NULL) = 0;

virtual void SaveOperator(const std::string filename) = 0;
virtual void LoadOperatorFromFile(const std::string filename) = 0;
virtual void SetRomMat(BlockMatrix *input_mat) = 0;
virtual void SetRomMat(BlockMatrix *input_mat, const bool init_direct_solver=true) = 0;
virtual void SaveRomSystem(const std::string &input_prefix, const std::string type="mm") = 0;

virtual void SaveBasisVisualization(const Array<FiniteElementSpace *> &fes, const std::vector<std::string> &var_names) = 0;
Expand Down Expand Up @@ -217,9 +229,6 @@ class MFEMROMHandler : public ROMHandlerBase
HYPRE_BigInt sys_row_starts[2];
HypreParMatrix *romMat_hypre = NULL;
MUMPSSolver *mumps = NULL;

mfem::BlockVector *reduced_rhs = NULL;
mfem::BlockVector *reduced_sol = NULL;

public:
MFEMROMHandler(TopologyHandler *input_topol, const Array<int> &input_var_offsets,
Expand All @@ -241,9 +250,9 @@ class MFEMROMHandler : public ROMHandlerBase
virtual void ProjectToDomainBasis(const Array<int> &idx_i, const Array<int> &idx_j, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);

virtual void ProjectComponentToRefBasis(const int &c, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
virtual void ProjectComponentToDomainBasis(const int &c, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
virtual void ProjectComponentToDomainBasis(const int &m, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
virtual void ProjectInterfaceToRefBasis(const int &c1, const int &c2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
virtual void ProjectInterfaceToDomainBasis(const int &c1, const int &c2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
virtual void ProjectInterfaceToDomainBasis(const int &m1, const int &m2, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
virtual void ProjectVariableToDomainBasis(const int &vi, const int &vj, const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
virtual void ProjectGlobalToDomainBasis(const Array2D<Operator*> &mats, Array2D<SparseMatrix *> &rom_mats);
virtual void ProjectOperatorOnReducedBasis(const Array2D<Operator*> &mats);
Expand All @@ -258,12 +267,13 @@ class MFEMROMHandler : public ROMHandlerBase
virtual void LiftUpFromDomainBasis(const int &i, const Vector &rom_vec, Vector &vec);
virtual void LiftUpGlobal(const BlockVector &rom_vec, BlockVector &vec);

virtual void Solve(BlockVector* U);
virtual void NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec=NULL) override;
void Solve(BlockVector &rhs, BlockVector &sol) override;
void Solve(BlockVector* U) override;
void NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec=NULL) override;

virtual void SaveOperator(const std::string input_prefix="");
virtual void LoadOperatorFromFile(const std::string input_prefix="");
virtual void SetRomMat(BlockMatrix *input_mat);
virtual void SetRomMat(BlockMatrix *input_mat, const bool init_direct_solver=true);
virtual void SaveRomSystem(const std::string &input_prefix, const std::string type="mm");

virtual void SaveBasisVisualization(const Array<FiniteElementSpace *> &fes, const std::vector<std::string> &var_names);
Expand Down
14 changes: 8 additions & 6 deletions include/steady_ns_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -60,7 +60,10 @@ class SteadyNSROM : public Operator
const int num_var = 2;
int numSub = -1;

/* block offsets of velocity */
Array<int> block_offsets;
/* block index of velocity on FOM variable */
Array<int> midxs;
Array<Array<int> *> block_idxs;
SparseMatrix *linearOp = NULL;

Expand All @@ -71,7 +74,7 @@ class SteadyNSROM : public Operator
HYPRE_BigInt sys_glob_size;
mutable HYPRE_BigInt sys_row_starts[2];
public:
SteadyNSROM(SparseMatrix *linearOp_, const int numSub_, const Array<int> &block_offsets_, const bool direct_solve_=true);
SteadyNSROM(const int numSub_, ROMHandlerBase *rom_handler, const bool direct_solve_=true);

virtual ~SteadyNSROM();

Expand All @@ -85,8 +88,8 @@ class SteadyNSTensorROM : public SteadyNSROM
Array<DenseTensor *> hs; // not owned by SteadyNSTensorROM.

public:
SteadyNSTensorROM(SparseMatrix *linearOp_, Array<DenseTensor *> &hs_, const Array<int> &block_offsets_, const bool direct_solve_=true)
: SteadyNSROM(linearOp_, hs_.Size(), block_offsets_, direct_solve_), hs(hs_) {}
SteadyNSTensorROM(ROMHandlerBase *rom_handler, Array<DenseTensor *> &hs_, const bool direct_solve_=true)
: SteadyNSROM(hs_.Size(), rom_handler, direct_solve_), hs(hs_) {}

virtual ~SteadyNSTensorROM() {}

Expand All @@ -101,13 +104,12 @@ class SteadyNSEQPROM : public SteadyNSROM
ROMInterfaceForm *itf = NULL; // not owned by SteadyNSEQPROM.

Array<int> u_offsets;
Array<int> u_idxs;

mutable Vector x_u, y_u;

public:
SteadyNSEQPROM(SparseMatrix *linearOp_, Array<ROMNonlinearForm *> &hs_, ROMInterfaceForm *itf_,
const Array<int> &block_offsets_, const bool direct_solve_=true);
SteadyNSEQPROM(ROMHandlerBase *rom_handler, Array<ROMNonlinearForm *> &hs_,
ROMInterfaceForm *itf_, const bool direct_solve_=true);

virtual ~SteadyNSEQPROM() {}

Expand Down
28 changes: 18 additions & 10 deletions include/unsteady_ns_solver.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,14 @@ friend class SteadyNSOperator;

BlockOperator *Hop = NULL;

/* function coefficients for initial condition */
VectorCoefficient *u_ic = NULL;
VectorCoefficient *p_ic = NULL;

/* mass matrix operator for ROM */
Array<int> rom_u_offsets;
BlockMatrix *rom_mass = NULL;

public:
UnsteadyNSSolver();

Expand All @@ -81,25 +89,23 @@ friend class SteadyNSOperator;

void SetParameterizedProblem(ParameterizedProblem *problem) override;

BlockVector* PrepareSnapshots(std::vector<BasisTag> &basis_tags) override;

void ProjectOperatorOnReducedBasis() override
{ mfem_error("UnsteadyNSSolver::ProjectOperatorOnReducedBasis is not implemented yet!\n"); }

void SolveROM() override
{ mfem_error("UnsteadyNSSolver::SolveROM is not implemented yet!\n"); }
void BuildCompROMLinElems() override;

void SolveROM() override;

void InitROMHandler() override;

void BuildROMTensorElems() override
{ mfem_error("UnsteadyNSSolver::BuildROMTensorElems is not implemented yet!\n"); }
void TrainROMEQPElems(SampleGenerator *sample_generator) override
{ mfem_error("UnsteadyNSSolver::TrainROMEQPElems is not implemented yet!\n"); }
void SaveROMNlinElems(const std::string &input_prefix) override
{ mfem_error("UnsteadyNSSolver::SaveROMNlinElems is not implemented yet!\n"); }
void LoadROMNlinElems(const std::string &input_prefix) override
{ mfem_error("UnsteadyNSSolver::LoadROMNlinElems is not implemented yet!\n"); }
void AssembleROMNlinOper() override
{ mfem_error("UnsteadyNSSolver::AssembleROMNlinOper is not implemented yet!\n"); }

private:
void InitializeTimeIntegration();
void SetupInitialCondition(int &initial_step, double &time);
void Step(double &time, int step);

void SanityCheck(const int step)
Expand All @@ -113,6 +119,8 @@ friend class SteadyNSOperator;
double ComputeCFL(const double dt);
void SetTime(const double time);

void AssembleROMMat(BlockMatrix &romMat) override;

};

#endif
4 changes: 3 additions & 1 deletion src/main_workflow.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -288,12 +288,14 @@ void AuxiliaryTrainROM(MPI_Comm comm, SampleGenerator *sample_generator)
bool separate_variable_basis = config.GetOption<bool>("model_reduction/separate_variable_basis", false);

/* Supremizer enrichment */
if ((separate_variable_basis) && ((solver_type == "stokes") || (solver_type == "steady-ns")))
if ((separate_variable_basis) &&
((solver_type == "stokes") || (solver_type == "steady-ns") || (solver_type == "unsteady-ns")))
{
ParameterizedProblem *problem = InitParameterizedProblem();
StokesSolver *solver = NULL;
if (solver_type == "stokes") { solver = new StokesSolver; }
else if (solver_type == "steady-ns") { solver = new SteadyNSSolver; }
else if (solver_type == "unsteady-ns") { solver = new UnsteadyNSSolver; }

if (!solver->UseRom()) mfem_error("ROM must be enabled for supremizer enrichment!\n");

Expand Down
67 changes: 25 additions & 42 deletions src/multiblock_solver.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -272,31 +272,33 @@ void MultiBlockSolver::BuildROMLinElems()

void MultiBlockSolver::AssembleROMMat()
{
assert(topol_mode == TopologyHandlerMode::COMPONENT);
assert(rom_elems);

const Array<int> *rom_block_offsets = rom_handler->GetBlockOffsets();
BlockMatrix *romMat = new BlockMatrix(*rom_block_offsets);
romMat->owns_blocks = true;

AssembleROMMat(*romMat);

romMat->Finalize();
rom_handler->SetRomMat(romMat);
}

void MultiBlockSolver::AssembleROMMat(BlockMatrix &romMat)
{
assert(topol_mode == TopologyHandlerMode::COMPONENT);
assert(rom_elems);

// component domain matrix.
for (int m = 0; m < numSub; m++)
{
int c_type = topol_handler->GetMeshType(m);
int num_block = 1;
int midx0 = m;
if (separate_variable_basis)
{
num_block *= num_var;
midx0 *= num_var;
}
int num_block = (separate_variable_basis) ? num_var : 1;

Array<int> midx(num_block);
for (int k = 0; k < num_block; k++)
midx[k] = midx0 + k;
for (int v = 0; v < num_block; v++)
midx[v] = rom_handler->GetBlockIndex(m, v);

MatrixBlocks *comp_mat = rom_elems->comp[c_type];
AddToBlockMatrix(midx, midx, *comp_mat, *romMat);
AddToBlockMatrix(midx, midx, *comp_mat, romMat);

// boundary matrices of each component.
Array<int> *bdr_c2g = topol_handler->GetBdrAttrComponentToGlobalMap(m);
Expand All @@ -312,7 +314,7 @@ void MultiBlockSolver::AssembleROMMat()
continue;

MatrixBlocks *bdr_mat = (*(rom_elems->bdr[c_type]))[b];
AddToBlockMatrix(midx, midx, *bdr_mat, *romMat);
AddToBlockMatrix(midx, midx, *bdr_mat, romMat);
} // for (int b = 0; b < bdr_c2g->Size(); b++)
} // for (int m = 0; m < numSub; m++)

Expand All @@ -325,38 +327,19 @@ void MultiBlockSolver::AssembleROMMat()

const int m1 = pInfo->Mesh1;
const int m2 = pInfo->Mesh2;
const int c1 = topol_handler->GetMeshType(m1);
const int c2 = topol_handler->GetMeshType(m2);

int num_block = 2;
int c1idx0 = c1, c2idx0 = c2;
int m1idx0 = m1, m2idx0 = m2;
if (separate_variable_basis)
{
num_block *= num_var;
c1idx0 *= num_var;
c2idx0 *= num_var;
m1idx0 *= num_var;
m2idx0 *= num_var;
}
int num_block = (separate_variable_basis) ? num_var : 1;

Array<int> midx(num_block), num_basis(num_block);
for (int k = 0, cidx = c1idx0, m = m1idx0; k < num_block/2; k++, cidx++, m++)
{
num_basis[k] = rom_handler->GetRefNumBasis(cidx);
midx[k] = m;
}
for (int k = num_block/2, cidx = c2idx0, m = m2idx0; k < num_block; k++, cidx++, m++)
{
num_basis[k] = rom_handler->GetRefNumBasis(cidx);
midx[k] = m;
}
Array<int> midx(0);

AddToBlockMatrix(midx, midx, *port_mat, *romMat);
}
for (int v = 0; v < num_block; v++)
midx.Append(rom_handler->GetBlockIndex(m1, v));

romMat->Finalize();
rom_handler->SetRomMat(romMat);
for (int v = 0; v < num_block; v++)
midx.Append(rom_handler->GetBlockIndex(m2, v));

AddToBlockMatrix(midx, midx, *port_mat, romMat);
}
}

void MultiBlockSolver::InitVisualization(const std::string& output_path)
Expand Down
Loading

0 comments on commit 97b5355

Please sign in to comment.