gokhalen
commented
1 year ago @chldkdtn @dylan-copeland Any ideas?
Closed gokhalen closed 1 year ago
gokhalen
commented
1 year ago @chldkdtn @dylan-copeland Any ideas?
dylan-copeland
commented
1 year ago @gokhalen if you would like to make a branch to share your code, we can help debug. One thing I am curious about is the dimensions you reported here:
??MR2D:Bad submatrix:i=0,j=0,m=1361,n=2000,M=1361,N=1361
Is n=2000 the full-order dimension? It seems this is coming from some linear algebra routine (lapack or scalapack?). Usually the dense linear algebra is only at the reduced level. Seeing a 1361^2 submatrix leads me to guess that you have 1361 as the number of basis vectors or hyperreduction samples? If so, maybe that is too large, resulting in numerical instability. Since numerical linear algebra is not exact, problems like this may arise if there are too many basis vectors. For example, if you include basis vectors corresponding to singular values extremely close to 0, there may be singularity. For diffusion problems, a small number of basis vectors (at most 10 to 20) and hyperreduction samples should be sufficient. That is just my guess. If you think something else is happening, I would probably need more information, ideally the code.
gokhalen
commented
1 year ago Hi @dylan-copeland,
Thanks for your response. Since git is blocked at work, I can't quickly create a new branch. My code is at the end of the message, since Github does not allow me to upload cpp files.
The code is deliberately kept simple in 4 parts - common part, offline, merge and online, even though this duplicates the time-loop. This keeps the online and offline parts completely separate. You can ignore the online phase, it is still work in progress.
I've placed it into the examples/prom directory of libROM and modified CMakelists.txt appropriately in order to get it to build.
My goal is to create a simple ROM for a linear problem, without hyper-reduction, and later modify it to include non-linear terms and use hyper-reduction.
Answers to your questions are:
1) 2000 is not the full order model dimension. 1361 is the FOM dimension.
2000 seems to be roughly equal to nsets*max_num_snapshots (for nsets=4 and max_num_snapshots=502)
When nsets=3 and max_num_snapshots=502, n=2000 changes to n=1500, again approximately equal to nsets*max_num_snapshots
2) 1361 is the number of dofs in the problem. -On the same machine I've merged two basis files 2 nsets for the non_linear_elasticity_global_rom problem, which has approx 30000 dofs (15000 for vel, and 1500 for def). -So, I doubt that the number of dofs is the problem.
3) How do I set the number of basis vectors to generate? Is this set in the second argument to options?
CAROM::Options options(fespace.GetTrueVSize(), max_num_snapshots, 1, update_right_SV);
max_num_snapshots is 502 for the current choice of parameters
Setting max_num_snapshots to 10 does not help, I get the error
Creating file: basis0_u.000000
P=0000000:Program abort called in file ``/home/test/libROM/lib/linalg/svd/SVD.h'' at line 169
P=0000000:ERROR MESSAGE:
P=0000000:Failed verify: d_max_time_intervals == -1 || num_time_intervals < d_max_time_intervals
--------------------------------------------------------------------------4) I'm running inside Kevin Seung Whan Chung's docker container
5) My command-lines are for basis generation are:
sudo mpirun -np 1 --allow-run-as-root ../linear_heat_conduction_nachiket_prom --offline -dt 0.0001 -tf 0.05 -s 14 -k 0.96 -id 0
sudo mpirun -np 1 --allow-run-as-root ../linear_heat_conduction_nachiket_prom --offline -dt 0.0001 -tf 0.05 -s 14 -k 0.90 -id 1
sudo mpirun -np 1 --allow-run-as-root ../linear_heat_conduction_nachiket_prom --offline -dt 0.0001 -tf 0.05 -s 14 -k 0.85 -id 2
sudo mpirun -np 1 --allow-run-as-root ../linear_heat_conduction_nachiket_prom --offline -dt 0.0001 -tf 0.05 -s 14 -k 0.80 -id 36) My command-lines for merge are
for -ns 2 (which works):
sudo mpirun -np 1 --allow-run-as-root ../linear_heat_conduction_nachiket_prom --merge -ns 2 -dt 0.0001 -tf 0.05 -s 14for -ns 3 (which fiails)
sudo mpirun -np 1 --allow-run-as-root ../linear_heat_conduction_nachiket_prom --merge -ns 3 -dt 0.0001 -tf 0.05 -s 14Error message
Creating file: basisu.000000
??MR2D:Bad submatrix:i=0,j=0,m=1361,n=1500,M=1361,N=1361for -ns 4 (which fails)
sudo mpirun -np 1 --allow-run-as-root ../linear_heat_conduction_nachiket_prom --merge -ns 4 -dt 0.0001 -tf 0.05 -s 14Error Message
Creating file: basisu.000000
??MR2D:Bad submatrix:i=0,j=0,m=1361,n=2000,M=1361,N=13617) My code follows
// This does the FOM, and ROM (without hyper-reduction)
// for the linear heat equation. Based of mfem non-linear heat conduction
// example, ex16p.cpp. The ROM implementation is based on the
// non_linear_elasticity_global_rom example from libROM
// We will first do the explicit implementation then implicit
// This code only includes the exlicit implementation
// This example should be run with one process only
// Disclaimer: This example does not work
#include "mfem.hpp"
#include "linalg/Vector.h"
#include "linalg/BasisGenerator.h"
#include "linalg/BasisReader.h"
#include "hyperreduction/DEIM.h"
#include "hyperreduction/GNAT.h"
#include "hyperreduction/S_OPT.h"
#include "mfem/SampleMesh.hpp"
#include "mfem/Utilities.hpp"
#include <fstream>
#include <iostream>
#include <memory>
#include <cmath>
#include <limits>
using namespace std;
using namespace mfem;
/** After spatial discretization, the conduction model can be written as:
*
* du/dt = M^{-1}(-Ku)
*
* where u is the vector representing the temperature, M is the mass matrix,
* and K is the diffusion operator with diffusivity depending on u:
* (\kappa + \alpha u). We do not use alpha here. We are doing linear heat cond. alpha = 0.0
*
* Class ConductionOperator represents the right-hand side of the above ODE.
*/
class ConductionOperator : public TimeDependentOperator
{
protected:
ParFiniteElementSpace &fespace;
Array<int> ess_tdof_list; // this list remains empty for pure Neumann b.c.
ParBilinearForm *M;
ParBilinearForm *K;
double current_dt;
HypreSmoother M_prec; // Preconditioner for the mass matrix M
CGSolver T_solver; // Implicit solver for T = M + dt K
HypreSmoother T_prec; // Preconditioner for the implicit solver
double alpha, kappa;
mutable Vector z; // auxiliary vector
public:
ConductionOperator(ParFiniteElementSpace &f, double alpha, double kappa,
const Vector &u);
virtual void Mult(const Vector &u, Vector &du_dt) const;
/** Solve the Backward-Euler equation: k = f(u + dt*k, t), for the unknown k.
This is the only requirement for high-order SDIRK implicit integration.*/
// virtual void ImplicitSolve(const double dt, const Vector &u, Vector &k);
/// Update the diffusion BilinearForm K using the given true-dof vector `u`.
// void SetParameters(const Vector &u);
mutable Vector du_dt_sp;
HypreParMatrix Mmat;
HypreParMatrix Kmat;
HypreParMatrix *T; // T = M + dt K
CGSolver M_solver; // Krylov solver for inverting the mass matrix M
virtual ~ConductionOperator();
};
class RomOperator : public TimeDependentOperator
{
private:
int rtdim;
protected:
ConductionOperator *fomSp;
CAROM::Matrix *K_hat;
CAROM::Matrix *M_hat;
CAROM::Matrix *M_hat_inv;
public:
ConductionOperator *fom;
RomOperator(ConductionOperator* fom_, ConductionOperator* fomSp_, const int rtdim,
const Vector* u0_, const Vector u0_fom_, const CAROM::Matrix* V_u_);
virtual void Mult(const Vector& u, Vector& du_dt) const;
void Mult_Hyperreduced(const Vector& u, Vector& du_dt) const;
void Mult_FullOrder(const Vector& u, Vector& du_dt) const;
CAROM::Matrix V_u;
const Vector* u0;
Vector u0_fom;
mutable Vector z_mut_rom; // auxiliary vector
bool hyperreduce = false;
virtual ~RomOperator();
};
// TODO: move this to the library?
CAROM::Matrix* GetFirstColumns(const int N, const CAROM::Matrix* A)
{
CAROM::Matrix* S = new CAROM::Matrix(A->numRows(), std::min(N, A->numColumns()),
A->distributed());
for (int i = 0; i < S->numRows(); ++i)
{
for (int j = 0; j < S->numColumns(); ++j)
(*S)(i, j) = (*A)(i, j);
}
// delete A; // TODO: find a good solution for this.
return S;
}
// TODO: move this to the library?
void BasisGeneratorFinalSummary(CAROM::BasisGenerator* bg,
const double energyFraction, int& cutoff, const std::string cutoffOutputPath)
{
const int rom_dim = bg->getSpatialBasis()->numColumns();
const CAROM::Vector* sing_vals = bg->getSingularValues();
MFEM_VERIFY(rom_dim <= sing_vals->dim(), "");
double sum = 0.0;
for (int sv = 0; sv < sing_vals->dim(); ++sv) {
sum += (*sing_vals)(sv);
}
vector<double> energy_fractions = { 0.9999999, 0.999999, 0.99999, 0.9999, 0.999, 0.99, 0.9 };
bool reached_cutoff = false;
ofstream outfile(cutoffOutputPath);
double partialSum = 0.0;
for (int sv = 0; sv < sing_vals->dim(); ++sv) {
partialSum += (*sing_vals)(sv);
for (int i = energy_fractions.size() - 1; i >= 0; i--)
{
if (partialSum / sum > energy_fractions[i])
{
outfile << "For energy fraction: " << energy_fractions[i] << ", take first "
<< sv + 1 << " of " << sing_vals->dim() << " basis vectors" << endl;
energy_fractions.pop_back();
}
else
{
break;
}
}
if (!reached_cutoff && partialSum / sum > energyFraction)
{
cutoff = sv + 1;
reached_cutoff = true;
}
}
if (!reached_cutoff) cutoff = sing_vals->dim();
outfile << "Take first " << cutoff << " of " << sing_vals->dim() <<
" basis vectors" << endl;
outfile.close();
}
void MergeBasis(const int dimFOM, const int nparam, const int max_num_snapshots,
std::string name)
{
MFEM_VERIFY(nparam > 0, "Must specify a positive number of parameter sets");
bool update_right_SV = false;
bool isIncremental = false;
CAROM::Options options(dimFOM, nparam * max_num_snapshots, 1, update_right_SV);
CAROM::BasisGenerator generator(options, isIncremental, "basis" + name);
for (int paramID = 0; paramID < nparam; ++paramID)
{
std::string snapshot_filename = "basis" + std::to_string(
paramID) + "_" + name + "_snapshot";
generator.loadSamples(snapshot_filename, "snapshot");
}
generator.endSamples(); // save the merged basis file
int cutoff = 0;
BasisGeneratorFinalSummary(&generator, 0.9999999, cutoff,
"mergedSV_" + name + ".txt");
}
double InitialTemperature(const Vector &x);
int main(int argc, char *argv[])
{
// 1. Initialize MPI and HYPRE.
Mpi::Init(argc, argv);
int num_procs = Mpi::WorldSize();
int myid = Mpi::WorldRank();
Hypre::Init();
// 2. Parse command-line options.
const char *mesh_file = "../../data/star.mesh";
int ser_ref_levels = 2;
// We are not doing any parallel refinement, hence set to zero
int par_ref_levels = 0;
int order = 2;
int ode_solver_type = 14;
double t_final = 0.5;
double dt = 1.0e-2;
// non-linearity in heat conduction is set to zero
// because we are doing linear heat conduction
const double alpha = 0.0;
double kappa = 0.5;
bool visualization = true;
bool visit = false;
int vis_steps = 5;
int precision = 8;
cout.precision(precision);
// ROM parameters
bool offline = false;
bool merge = false;
bool online = false;
bool hyperreduce = false;
// number of basis sets.
// If we want to use 3 FOM simulations nsets = 3
int nsets = 0;
int id_param = 0;
bool update_right_SV = false;
bool isIncremental = false;
// number of basis vectors to use for the temperature
int rtdim = -1;
OptionsParser args(argc, argv);
args.AddOption(&mesh_file, "-m", "--mesh",
"Mesh file to use.");
args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
"Number of times to refine the mesh uniformly in serial.");
args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
"Number of times to refine the mesh uniformly in parallel.");
args.AddOption(&order, "-o", "--order",
"Order (degree) of the finite elements.");
args.AddOption(&ode_solver_type, "-s", "--ode-solver",
"ODE solver: 1 - Backward Euler, 2 - SDIRK2, 3 - SDIRK3,\n\t"
"\t 11 - Forward Euler, 12 - RK2, 13 - RK3 SSP, 14 - RK4.");
args.AddOption(&t_final, "-tf", "--t-final",
"Final time; start time is 0.");
args.AddOption(&dt, "-dt", "--time-step",
"Time step.");
args.AddOption(&kappa, "-k", "--kappa",
"Kappa coefficient offset.");
args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
"--no-visualization",
"Enable or disable GLVis visualization.");
args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
"--no-visit-datafiles",
"Save data files for VisIt (visit.llnl.gov) visualization.");
args.AddOption(&vis_steps, "-vs", "--visualization-steps",
"Visualize every n-th timestep.");
args.AddOption(&nsets, "-ns", "--nset", "Number of parametric snapshot sets");
args.AddOption(&offline, "-offline", "--offline", "-no-offline", "--no-offline",
"Enable or disable the offline phase.");
args.AddOption(&online, "-online", "--online", "-no-online", "--no-online",
"Enable or disable the online phase.");
args.AddOption(&merge, "-merge", "--merge", "-no-merge", "--no-merge",
"Enable or disable the merge phase.");
args.AddOption(&rtdim, "-rtdim", "--rtdim",
"Basis dimension for temperature solution space.");
args.AddOption(&id_param, "-id", "--id", "Parametric index");
args.AddOption(&hyperreduce, "-hyp", "--hyperreduce", "-no-hyp",
"--no-hyperreduce", "Hyper reduce nonlinear term");
args.Parse();
if (!args.Good())
{
args.PrintUsage(cout);
return 1;
}
if (myid == 0)
{
args.PrintOptions(cout);
}
// dependent quantities
int max_num_snapshots = int(t_final / dt) + 2;
if (myid == 0)
cout<<"max_num_snapshots = "<< max_num_snapshots << endl;
const std::string basisFileName = "basis" + std::to_string(id_param);
const bool check = (offline && !merge && !online) || (!offline && merge
&& !online) || (!offline && !merge && online);
MFEM_VERIFY(check, "only one of offline, merge, or online must be true!");
StopWatch solveTimer, offlineTimer, mergeTimer, onlineTimer, totalTimer;
totalTimer.Start();
// 3. Read the serial mesh from the given mesh file on all processors. We can
// handle triangular, quadrilateral, tetrahedral and hexahedral meshes
// with the same code.
Mesh *mesh = new Mesh(mesh_file, 1, 1);
int dim = mesh->Dimension();
// 4. Define the ODE solver used for time integration. Several implicit
// singly diagonal implicit Runge-Kutta (SDIRK) methods, as well as
// explicit Runge-Kutta methods are available.
// We are simply looking at explicit methods now, because of simplicity
// We will look at implicit methods later
// All runs have been carried out with integrator 14 - RK4Solver
ODESolver *ode_solver;
switch (ode_solver_type)
{
// Explicit methods
case 11: ode_solver = new ForwardEulerSolver; break;
case 12: ode_solver = new RK2Solver(0.5); break; // midpoint method
case 13: ode_solver = new RK3SSPSolver; break;
case 14: ode_solver = new RK4Solver; break;
case 15: ode_solver = new GeneralizedAlphaSolver(0.5); break;
default:
if (myid == 0)
cout << "Unknown ODE solver type: " << ode_solver_type << '\n';
delete mesh;
MPI_Finalize();
return 3;
}
// 5. Refine the mesh in serial to increase the resolution. In this example
// we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
// a command-line parameter.
for (int lev = 0; lev < ser_ref_levels; lev++)
{
mesh->UniformRefinement();
}
// 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
// this mesh further in parallel to increase the resolution. Once the
// parallel mesh is defined, the serial mesh can be deleted.
ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
delete mesh;
for (int lev = 0; lev < par_ref_levels; lev++)
{
pmesh->UniformRefinement();
}
// 7. Define the vector finite element space representing the current and the
// initial temperature, u_ref.
H1_FECollection fe_coll(order, dim);
ParFiniteElementSpace fespace(pmesh, &fe_coll);
HYPRE_BigInt fe_size = fespace.GlobalTrueVSize();
if (myid == 0)
{
cout << "Number of temperature unknowns: " << fe_size << endl;
}
ParGridFunction u_gf(&fespace);
// 8. Set the initial conditions for u. All boundaries are considered
// natural.
FunctionCoefficient u_0(InitialTemperature);
u_gf.ProjectCoefficient(u_0);
Vector u;
u_gf.GetTrueDofs(u);
// Initialize and copy the initial solution
Vector u0 = Vector(u);
Vector u_diff = Vector(u);
// 9. Offline Phase
if (offline)
{
offlineTimer.Start();
// 9.1 Initialize various things
cout<<"Entering offline phase..."<<endl;
ConductionOperator oper(fespace,alpha,kappa,u);
// Declare du_dt - it is not filled as of now
Vector du_dt(fe_size);
// Declare basis generators
CAROM::BasisGenerator* basis_generator_u = 0;
// Set options
cout<<"fespace.GetTrueVSize()= "<<fespace.GetTrueVSize()<<endl;
cout<<"fespace.GlobalTrueVSize()= "<<fespace.GlobalTrueVSize()<<endl;
CAROM::Options options(fespace.GetTrueVSize(), max_num_snapshots, 1, update_right_SV);
// Create basis generator
basis_generator_u = new CAROM::BasisGenerator(options, isIncremental, basisFileName + "_u");
// Write initial solution
u_gf.SetFromTrueDofs(u);
{
ostringstream mesh_name, sol_name;
mesh_name << "LHC_NHG_prom_offline-mesh." << setfill('0') << setw(6) << myid;
sol_name << "LHC_NHG_prom_offline-init." << setfill('0') << setw(6) << myid;
ofstream omesh(mesh_name.str().c_str());
omesh.precision(precision);
pmesh->Print(omesh);
ofstream osol(sol_name.str().c_str());
osol.precision(precision);
u_gf.Save(osol);
}
// 9.2 Perform time-integration (looping over the time iterations, ti, with a
// time-step dt).
ode_solver->Init(oper);
double t = 0.0;
bool last_step = false;
for (int ti = 1; !last_step; ti++)
{
if ( t + dt >= t_final - dt/2)
{
last_step = true;
}
ode_solver->Step(u, t, dt);
if ( basis_generator_u->isNextSample(t) || last_step )
{
du_dt = oper.du_dt_sp.GetData();
u_diff = Vector(u);
u_diff -= u0;
// See mixed_nonlinear_diffusion.cpp in libRROM for this usage of GetData
basis_generator_u->takeSample(u_diff.GetData(), t, dt);
basis_generator_u->computeNextSampleTime(u_diff.GetData(), du_dt.GetData(),t);
}
if ( last_step || (ti % vis_steps) == 0 )
{
if (myid == 0)
{
cout << " step " << ti << ", t = "<< t << endl;
}
u_gf.SetFromTrueDofs(u);
}
if (last_step)
{
basis_generator_u->writeSnapshot();
delete basis_generator_u;
}
} // time loop ti
// 9.3 Save the final solution in parallel.
// s output can be viewed later using GLVis: "glvis -np <np> -m ex16-mesh -g ex16-final".
{
ostringstream sol_name;
sol_name << "LHC_NHG_prom_offline-final." << setfill('0') << setw(6) << myid;
ofstream osol(sol_name.str().c_str());
osol.precision(precision);
u_gf.Save(osol);
}
offlineTimer.Stop();
cout<<"Elapsed time for offline phase " << offlineTimer.RealTime() << "s" << endl << flush;
}// offline
// 10. Merge Phase
if (merge)
{
mergeTimer.Start();
// 10.1 Initialize various things
cout<<"Entering merge phase...nsets = "<<nsets<<" ..fe_size = "<<fe_size<<endl;
MergeBasis(fe_size, nsets, max_num_snapshots, "u");
mergeTimer.Stop();
cout<<" Elapsed time for merge phase " << mergeTimer.RealTime() << "s" << endl << flush;
MPI_Finalize();
return 0;
}
// 11. Online Phase
if (online)
{
onlineTimer.Start();
// 11.1 Declare and Initialize various things
const CAROM::Matrix* BU_librom = 0;
CAROM::BasisReader* readerU = 0;
RomOperator* romop = 0;
ConductionOperator oper(fespace,alpha,kappa,u);
// 11.2 Read bases
readerU = new CAROM:: BasisReader("basisu");
BU_librom = readerU->getSpatialBasis(0.0);
if (rtdim == -1) // Change rtdim
{
rtdim = BU_librom->numColumns();
cout<< "rtdim set from BU_librom->numColumns() = "<<rtdim << endl;
}
else
{
BU_librom = GetFirstColumns(rtdim, BU_librom);
cout<< "rtdim after GetFirstColumns = "<<rtdim << endl;
}
MFEM_VERIFY(BU_librom->numRows() == fe_size, "Number of rows of BU_librom not equal to fe_size");
// initialize rom operator
// fomSp_ = fom_ (we're not using fomSp_ here)
// fomSp_ is just a dummy
romop = new RomOperator(&oper, &oper, rtdim, &u0, u0, BU_librom);
// write the initial solution
u_gf.SetFromTrueDofs(u);
{
ostringstream mesh_name, sol_name;
mesh_name << "LHC_NHG_prom_online-mesh." << setfill('0') << setw(6) << myid;
sol_name << "LHC_NHG_prom_online-init." << setfill('0') << setw(6) << myid;
ofstream omesh(mesh_name.str().c_str());
omesh.precision(precision);
pmesh->Print(omesh);
ofstream osol(sol_name.str().c_str());
osol.precision(precision);
u_gf.Save(osol);
}
// 11.3 time integration
ode_solver->Init(*romop);
double t = 0.0;
bool last_step = false;
for (int ti = 1; !last_step; ti++)
{
if ( t + dt >= t_final - dt/2)
{
last_step = true;
}
ode_solver->Step(u,t,dt);
// set u_gf
if ( last_step || (ti % vis_steps) == 0 )
{
if (myid == 0)
{
cout << " step " << ti << ", t = "<< t << endl;
}
u_gf.SetFromTrueDofs(u);
}
} // time integration - online phase
// write final online solution
{
ostringstream sol_name;
sol_name << "LHC_NHG_prom_online-final." << setfill('0') << setw(6) << myid;
ofstream osol(sol_name.str().c_str());
osol.precision(precision);
u_gf.Save(osol);
}
onlineTimer.Stop();
cout<<" Elapsed time for online phase " << onlineTimer.RealTime() << "s" << endl << flush;
}//online phase
totalTimer.Stop();
cout << "Elapsed time for entire simulation " << totalTimer.RealTime() << "s"<< endl << flush;
return 0;
}
ConductionOperator::ConductionOperator(ParFiniteElementSpace &f, double al,
double kap, const Vector &u)
: TimeDependentOperator(f.GetTrueVSize(), 0.0), fespace(f), M(NULL), K(NULL),
T(NULL), current_dt(0.0), M_solver(f.GetComm()), T_solver(f.GetComm()), z(height)
{
const double rel_tol = 1e-8;
M = new ParBilinearForm(&fespace);
M->AddDomainIntegrator(new MassIntegrator());
M->Assemble(0); // keep sparsity pattern of M and K the same
M->FormSystemMatrix(ess_tdof_list, Mmat);
M_solver.iterative_mode = false;
M_solver.SetRelTol(rel_tol);
M_solver.SetAbsTol(0.0);
M_solver.SetMaxIter(100);
M_solver.SetPrintLevel(0);
M_prec.SetType(HypreSmoother::Jacobi);
M_solver.SetPreconditioner(M_prec);
M_solver.SetOperator(Mmat);
alpha = al;
kappa = kap;
K = new ParBilinearForm(&fespace);
ConstantCoefficient u_coeff(kappa);
K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
K->Assemble(0);
K->FormSystemMatrix(ess_tdof_list,Kmat);
T_solver.iterative_mode = false;
T_solver.SetRelTol(rel_tol);
T_solver.SetAbsTol(0.0);
T_solver.SetMaxIter(100);
T_solver.SetPrintLevel(0);
T_solver.SetPreconditioner(T_prec);
// SetParameters(u); // We don't have it anymore
}
//------------------------------------------------------------------------------
void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
{
// Compute:
// du_dt = M^{-1}*-Ku
// for du_du, where K is lineaized by using u from the previous timestep
Kmat.Mult(u,z);
z.Neg(); // z = -z
M_solver.Mult(z, du_dt);
du_dt_sp = du_dt;
}
//------------------------------------------------------------------------------
/*
void ConductionOperator::SetParameters(const Vector &u)
{
ParGridFunction u_alpha_gf(&fespace);
u_alpha_gf.SetFromTrueDofs(u);
for ( int i = 0; i < u_alpha_gf.Size(); i++)
{
u_alpha_gf(i) = kappa + alpha*u_alpha_gf(i);
}
delete K;
K = new ParBilinearForm(&fespace);
GridFunctionCoefficient u_coeff(&u_alpha_gf);
K->AddDomainIntegrator(new DiffusionIntegrator(u_coeff));
K->Assemble(0); // keep sparsity pattern of M and K the same
K->FormSystemMatrix(ess_tdof_list, Kmat);
delete T;
T = NULL; // re-compute T on the next ImplicitSolve
}*/
//------------------------------------------------------------------------------
ConductionOperator::~ConductionOperator()
{
delete T;
delete M;
delete K;
}
//------------------------------------------------------------------------------
double InitialTemperature(const Vector &x)
{
if (x.Norml2() < 0.5)
{
return 2.0;
}
else
{
return 1.0;
}
}
//------------------------------------------------------------------------------
RomOperator::RomOperator(ConductionOperator* fom_, ConductionOperator* fomSp_,
const int rtdim,
const Vector* u0_, const Vector u0_fom_, const CAROM::Matrix* V_u_)
: TimeDependentOperator(1361, 0.0), fom(fom_), fomSp(fomSp_), rtdim(rtdim),
u0(u0_), u0_fom(u0_fom_), V_u(*V_u_), z_mut_rom(u0->Size())
{
// fom_ : full order model
// fomSp_ : full order model on sample mesh not used now
// rtdim : number of basis vectors for temperature
// u0_ : u0 in the full dimension through a pointer
// u0_fom_ : u0 in the full dimension through a variable
cout << " in RomOperator Constructor rtdim == " << rtdim;
K_hat = new CAROM::Matrix(rtdim,rtdim,false);
M_hat = new CAROM::Matrix(rtdim,rtdim,false);
M_hat_inv = new CAROM::Matrix(rtdim,rtdim,false);
// if we look at ex10.cpp in the mfem examples: the hyperelastic system is
// dv/dt = -M^{-1}*(H(x) + S*v)
// dx/dt = v,
// if we look at ex16.cpp the heat conduction system is
// du/dt = M^{-1}(-Ku)
// Since there is a negative on both sides, the S matrix in the ROM
// translates to K for conduction
// Create K_hat
ComputeCtAB(fom->Kmat, V_u, V_u, *K_hat);
// Create M_hat
ComputeCtAB(fom->Mmat, V_u, V_u, *M_hat);
// Invert M_hat and store
M_hat->inverse(*M_hat_inv);
}
//------------------------------------------------------------------------------
void RomOperator::Mult(const Vector& u, Vector& du_dt) const
{
if (hyperreduce)
Mult_Hyperreduced(u, du_dt);
else
Mult_FullOrder(u, du_dt);
}
//-----------------------------------------------------------------------------
void RomOperator::Mult_FullOrder(const Vector& u, Vector& du_dt) const
{
cout<<"Calling Mult_FullOrder u0->Size()= "<<u0->Size()<<" u.Size()="<<u.Size()<<" du_dt.Size()="<<du_dt.Size()<<endl;
// this works
fom->Mult(u,du_dt);
// this blows up
//fom->Kmat.Mult(u,z_mut_rom);
//z_mut_rom.Neg();
//fom->M_solver.Mult(z_mut_rom, du_dt);
}
//-----------------------------------------------------------------------------
void RomOperator::Mult_Hyperreduced(const Vector& u, Vector& du_dt) const
{
}
//------------------------------------------------------------------------------
RomOperator::~RomOperator()
{
delete K_hat;
delete M_hat;
delete M_hat_inv;
}
//-----------------------------------------------------------------------------
Hi @gokhalen, the problem occurs when the number of snapshots exceeds the FOM dimension. In this example, the FOM dimension is 1361, and each parameter has 500 snapshots. For -ns 2, 1000 snapshots works, but for -ns 3, 1500 > 1361 fails. Subsampling the snapshots for each parameter would enable keeping the total number of snapshots below the FOM dimension.
gokhalen
commented
1 year ago Hi @dylan-copeland
Thanks for your response.
1) Thanks. How do I set the number of snapshots taken? Do I have to do this manually by controlling when takeSample is called?
It seems that simply setting max_num_snapshots=100 or changing the max_num_snapshots to 100 ONLY in the arguments to
CAROM::Options options(fespace.GetTrueVSize(), max_num_snapshots, 1, update_right_SV);Leads to an error in the offline phase itself.
Creating file: basis0_u.000000
P=0000000:Program abort called in file ``/home/test/libROM/lib/linalg/svd/SVD.h'' at line 169
P=0000000:ERROR MESSAGE:
P=0000000:Failed verify: d_max_time_intervals == -1 || num_time_intervals < d_max_time_intervals2) Is there a mathematical reason why total number of snapshots cannot be greater than FOM dimension? It seems to me that we should be able to take as many snapshots as we want.
dylan-copeland
commented
1 year ago Options, along with isNextSample to determine when takeSample should be called. This basically checks for how much the snapshot variable is changing over time, to give an indicator of when a new sample should be taken with takeSample. I am not aware of any examples where this works well. In my experience, it is best just to subsample manually, by calling takeSample every 1 out of N times, where N is a parameter that you set in your simulation.gokhalen
commented
1 year ago Thank you, this helped.
gokhalen
commented
1 year ago Thank you, it seems I have a simple ROM for the linear heat conduction equation.
I'm trying to create a simple ROM for a transient heat conduction problem based on ex16.cpp in MFEM. I'm trying to base my ROM on nonlinear_elasticity_global_rom.cpp.
I seem to be able to generate basis snapshots but I have a problem merging the basis. With 2 basis snapshots generated by 2 offline runs, the merge is fine (at least there are no errors reported by libROM), but with 3 or 4 basis snapshots, I get the error message
??MR2D:Bad submatrix:i=0,j=0,m=1361,n=2000,M=1361,N=1361Can someone tell me what is happening here?
Edit: FWIW, the basis files seem to have non-zero size.
When the merge succeeds, with 2 basis snapshots, I get
When the merge fails, with 3/4 basis snapshots, I get
Thanks,
Nachiket