Plug in rom interface form into SteadyNSEQPROM. Still needs to handle…

… boundaries...

include/steady_ns_solver.hpp

@@ -97,15 +97,25 @@ class SteadyNSEQPROM : public SteadyNSROM
 {
 protected:
  Array<ROMNonlinearForm *> hs; // not owned by SteadyNSEQPROM.
+ 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_, const Array<int> &block_offsets_, const bool direct_solve_=true)
- : SteadyNSROM(linearOp_, hs_.Size(), block_offsets_, direct_solve_), hs(hs_) {}
+ SteadyNSEQPROM(SparseMatrix *linearOp_, Array<ROMNonlinearForm *> &hs_, ROMInterfaceForm *itf_,
+  const Array<int> &block_offsets_, const bool direct_solve_=true);
 
  virtual ~SteadyNSEQPROM() {}
 
  virtual void Mult(const Vector &x, Vector &y) const;
  virtual Operator &GetGradient(const Vector &x) const;
+
+private:
+ void GetVel(const Vector &x, Vector &x_u) const;
+ void AddVel(const Vector &y_u, Vector &y) const;
 };
 
 class SteadyNSSolver : public StokesSolver

src/linalg_utils.cpp

@@ -137,8 +137,8 @@ namespace mfem
 void AddToBlockMatrix(const Array<int> &ridx, const Array<int> &cidx, const MatrixBlocks &mats, BlockMatrix &bmat)
 {
  assert(bmat.owns_blocks);
- assert(bmat.NumRowBlocks() >= mats.nrows);
- assert(bmat.NumColBlocks() >= mats.ncols);
+ // assert(bmat.NumRowBlocks() >= mats.nrows);
+ // assert(bmat.NumColBlocks() >= mats.ncols);
  assert(ridx.Size() == mats.nrows);
  assert(cidx.Size() == mats.ncols);
  assert((ridx.Min() >= 0) && (ridx.Max() < bmat.NumRowBlocks()));

src/steady_ns_solver.cpp

@@ -214,6 +214,37 @@ Operator& SteadyNSTensorROM::GetGradient(const Vector &x) const
  SteadyNSTensorROM
 */
 
+SteadyNSEQPROM::SteadyNSEQPROM(
+ SparseMatrix *linearOp_, Array<ROMNonlinearForm *> &hs_, ROMInterfaceForm *itf_,
+ const Array<int> &block_offsets_, const bool direct_solve_)
+ : SteadyNSROM(linearOp_, hs_.Size(), block_offsets_, direct_solve_), hs(hs_), itf(itf_)
+{
+ if (!separate_variable)
+ {
+ u_offsets = block_offsets;
+ x_u.SetSize(block_offsets.Last());
+ y_u.SetSize(block_offsets.Last());
+ u_idxs.SetSize(numSub);
+ for (int m = 0; m < numSub; m++) u_idxs[m] = m;
+ return;
+ }
+
+ u_offsets.SetSize(numSub+1);
+ u_offsets = 0;
+ for (int k = 0; k < numSub; k++)
+ {
+ u_offsets[k+1] = block_offsets[1 + k * num_var] - block_offsets[k * num_var];
+ }
+ u_offsets.PartialSum();
+
+ x_u.SetSize(u_offsets.Last());
+ y_u.SetSize(u_offsets.Last());
+
+ u_idxs.SetSize(numSub);
+ for (int m = 0; m < numSub; m++)
+ u_idxs[m] = m * num_var;
+}
+
 void SteadyNSEQPROM::Mult(const Vector &x, Vector &y) const
 {
  y = 0.0;
@@ -227,15 +258,50 @@ void SteadyNSEQPROM::Mult(const Vector &x, Vector &y) const
 
  hs[m]->AddMult(x_comp, y_comp);
  }
+
+ if (!itf) return;
+
+ GetVel(x, x_u);
+ y_u = 0.0;
+ itf->InterfaceAddMult(x_u, y_u);
+ AddVel(y_u, y);
 }
 
 Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const
 {
  delete jac_mono;
  delete jac_hypre;
  jac_mono = new SparseMatrix(*linearOp);
- DenseMatrix *jac_comp;
 
+ if (itf)
+ {
+ BlockMatrix jac(block_offsets);
+ jac.owns_blocks = true;
+
+ Array2D<SparseMatrix *> mats(numSub, numSub);
+ for (int i = 0; i < numSub; i++)
+ for (int j = 0; j < numSub; j++)
+ {
+ int iidx = (separate_variable) ? i * num_var : i;
+ int jidx = (separate_variable) ? j * num_var : j;
+ mats(i, j) = new SparseMatrix(block_offsets[iidx+1] - block_offsets[iidx],
+ block_offsets[jidx+1] - block_offsets[jidx]);
+ }
+
+ GetVel(x, x_u);
+ itf->InterfaceGetGradient(x_u, mats);
+
+ MatrixBlocks u_mats(mats);
+
+ AddToBlockMatrix(u_idxs, u_idxs, u_mats, jac);
+ jac.Finalize();
+ SparseMatrix *tmp = jac.CreateMonolithic();
+ (*jac_mono) += *tmp;
+
+ delete tmp;
+ }
+
+ DenseMatrix *jac_comp;
  for (int m = 0; m < numSub; m++)
  {
  int midx = (separate_variable) ? m * num_var : m;
@@ -256,6 +322,43 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const
  return *jac_mono;
 }
 
+void SteadyNSEQPROM::GetVel(const Vector &x, Vector &x_u) const
+{
+ if (!separate_variable)
+ {
+ x_u = x;
+ return;
+ }
+ x_u.SetSize(u_offsets.Last());
+
+ BlockVector x_block(const_cast<Vector &>(x).GetData(), block_offsets);
+ BlockVector u_block(x_u.GetData(), u_offsets);
+ for (int m = 0; m < numSub; m++)
+ {
+ Vector tmp;
+ u_block.GetBlockView(m, tmp);
+ tmp = x_block.GetBlock(m * num_var);
+ }
+}
+
+void SteadyNSEQPROM::AddVel(const Vector &y_u, Vector &y) const
+{
+ if (!separate_variable)
+ {
+ y += y_u;
+ return;
+ }
+
+ BlockVector u_block(const_cast<Vector &>(y_u).GetData(), u_offsets);
+ BlockVector y_block(y.GetData(), block_offsets);
+ for (int m = 0; m < numSub; m++)
+ {
+ Vector tmp;
+ y_block.GetBlockView(m * num_var, tmp);
+ tmp += u_block.GetBlock(m);
+ }
+}
+
 /*
  SteadyNSSolver
 */
@@ -621,7 +724,7 @@ void SteadyNSSolver::SolveROM()
  {
  assert(subdomain_eqps.Size() == numSub);
  for (int m = 0; m < numSub; m++) assert(subdomain_eqps[m]);
- rom_oper = new SteadyNSEQPROM(rom_handler->GetOperator(), subdomain_eqps, *(rom_handler->GetBlockOffsets()));
+ rom_oper = new SteadyNSEQPROM(rom_handler->GetOperator(), subdomain_eqps, itf_eqp, *(rom_handler->GetBlockOffsets()));
  }
 
  rom_handler->NonlinearSolve(*rom_oper, U_domain);
@@ -896,6 +999,7 @@ void SteadyNSSolver::AllocateROMEQPElems()
  rom_handler->GetReferenceBasis(idx, basis);
  itf_eqp->SetBasisAtComponent(c, *basis);
  }
+ itf_eqp->UpdateBlockOffsets();
  }
 }
 
@@ -1045,7 +1149,7 @@ void SteadyNSSolver::LoadEQPElems(const std::string &filename)
  // only one integrator exists in each nonlinear form.
  comp_eqps[c]->LoadEQPForIntegrator(IntegratorType::DOMAIN, 0, grp_id, dset_name + "_integ0");
  if ((oper_type == OperType::LF) && (full_dg))
- comp_eqps[c]->LoadEQPForIntegrator(IntegratorType::INTERIORFACE, 1, grp_id, dset_name + "_integ1");
+ comp_eqps[c]->LoadEQPForIntegrator(IntegratorType::INTERIORFACE, 0, grp_id, dset_name + "_integ1");
 
  if (comp_eqps[c]->PrecomputeMode())
  comp_eqps[c]->PrecomputeCoefficients();
@@ -1054,6 +1158,9 @@ void SteadyNSSolver::LoadEQPElems(const std::string &filename)
  errf = H5Gclose(grp_id);
  assert(errf >= 0);
 
+ if (oper_type == OperType::LF)
+ itf_eqp->LoadEQPForIntegrator(0, file_id, "interface_integ0");
+
  errf = H5Fclose(file_id);
  assert(errf >= 0);
 }

test/gmsh/test_multi_comp_workflow.cpp

@@ -235,10 +235,10 @@ TEST(SteadyNS_Workflow, InterfaceEQP)
  config.dict_["main"]["mode"] = "train_eqp";
  TrainEQP(MPI_COMM_WORLD);
 
- // printf("\nBuild ROM \n\n");
+ printf("\nBuild ROM \n\n");
 
- // config.dict_["main"]["mode"] = "build_rom";
- // BuildROM(MPI_COMM_WORLD);
+ config.dict_["main"]["mode"] = "build_rom";
+ BuildROM(MPI_COMM_WORLD);
 
  // config.dict_["main"]["mode"] = "single_run";
  // double error = SingleRun(MPI_COMM_WORLD, "test_output.h5");