From 9271b9a8ba57469bdb160c23fd4a5ea1b938a972 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Mon, 15 Jul 2024 06:29:57 -0700 Subject: [PATCH 01/20] minor fix --- src/main_workflow.cpp | 1 + src/rom_handler.cpp | 2 +- src/steady_ns_solver.cpp | 4 ++-- 3 files changed, 4 insertions(+), 3 deletions(-) diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index 6a1202c5..85a340cf 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -121,6 +121,7 @@ std::vector GetGlobalBasisTagList(const TopologyHandlerMode &topol_mod else if (solver_type == "stokes") var_list = StokesSolver::GetVariableNames(); else if (solver_type == "steady-ns") var_list = SteadyNSSolver::GetVariableNames(); else if (solver_type == "linelast") var_list = LinElastSolver::GetVariableNames(); + else if (solver_type == "unsteady-ns") var_list = UnsteadyNSSolver::GetVariableNames(); else { printf("Unknown MultiBlockSolver %s!\n", solver_type.c_str()); diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index bfeed644..fb17457b 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -639,7 +639,7 @@ void MFEMROMHandler::NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec else { for (int k = 0; k < reduced_sol->Size(); k++) - (*reduced_sol)(k) = 1.0e-1 * UniformRandom(); + (*reduced_sol)(k) = 1.0e-5 * UniformRandom(); } int maxIter = config.GetOption("solver/max_iter", 100); diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index 024d6443..4418dd2f 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -625,7 +625,7 @@ bool SteadyNSSolver::Solve(SampleGenerator *sample_generator) else { for (int k = 0; k < sol_byvar.Size(); k++) - sol_byvar(k) = UniformRandom(); + sol_byvar(k) = 1.0e-5 * UniformRandom(); } SteadyNSOperator oper(systemOp, hs, nl_itf, u_offsets, direct_solve); @@ -1360,4 +1360,4 @@ DenseTensor* SteadyNSSolver::GetReducedTensor(DenseMatrix *basis, FiniteElementS } // for (int i = 0; i < num_basis_c; i++) return tensor; -} \ No newline at end of file +} From 188c9ec8563584e4eda56826ff70fc6c29672b0a Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Wed, 31 Jul 2024 13:54:58 -0700 Subject: [PATCH 02/20] auxiliary train rom as a separate mode --- bin/main.cpp | 11 ++++++----- include/main_workflow.hpp | 2 +- src/main_workflow.cpp | 4 ++-- 3 files changed, 9 insertions(+), 8 deletions(-) diff --git a/bin/main.cpp b/bin/main.cpp index 4522a05b..2977c26f 100644 --- a/bin/main.cpp +++ b/bin/main.cpp @@ -28,11 +28,12 @@ int main(int argc, char *argv[]) { if (rank == 0) RunExample(); } else { - if (mode == "sample_generation") GenerateSamples(MPI_COMM_WORLD); - else if (mode == "build_rom") BuildROM(MPI_COMM_WORLD); - else if (mode == "train_rom") TrainROM(MPI_COMM_WORLD); - else if (mode == "train_eqp") TrainEQP(MPI_COMM_WORLD); - else if (mode == "single_run") double dump = SingleRun(MPI_COMM_WORLD, output_file); + if (mode == "sample_generation") GenerateSamples(MPI_COMM_WORLD); + else if (mode == "build_rom") BuildROM(MPI_COMM_WORLD); + else if (mode == "train_rom") TrainROM(MPI_COMM_WORLD); + else if (mode == "aux_train_rom") AuxiliaryTrainROM(MPI_COMM_WORLD); + else if (mode == "train_eqp") TrainEQP(MPI_COMM_WORLD); + else if (mode == "single_run") double dump = SingleRun(MPI_COMM_WORLD, output_file); else { if (rank == 0) printf("Unknown mode %s!\n", mode.c_str()); diff --git a/include/main_workflow.hpp b/include/main_workflow.hpp index 946701a2..e903b75d 100644 --- a/include/main_workflow.hpp +++ b/include/main_workflow.hpp @@ -24,7 +24,7 @@ void CollectSamplesByBasis(SampleGenerator *sample_generator, const std::string void BuildROM(MPI_Comm comm); void TrainROM(MPI_Comm comm); // supremizer-enrichment etc.. -void AuxiliaryTrainROM(MPI_Comm comm, SampleGenerator *sample_generator); +void AuxiliaryTrainROM(MPI_Comm comm); // EQP training, could include hypre-reduction optimization. void TrainEQP(MPI_Comm comm); // Input parsing routine to list out all snapshot files for training a basis. diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index 85a340cf..57032254 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -278,12 +278,12 @@ void TrainROM(MPI_Comm comm) sample_generator->FormReducedBasis(basis_prefix); - AuxiliaryTrainROM(comm, sample_generator); + AuxiliaryTrainROM(comm); delete sample_generator; } -void AuxiliaryTrainROM(MPI_Comm comm, SampleGenerator *sample_generator) +void AuxiliaryTrainROM(MPI_Comm comm) { std::string solver_type = config.GetRequiredOption("main/solver"); bool separate_variable_basis = config.GetOption("model_reduction/separate_variable_basis", false); From 199c9745bf4ba29fa236f977bd030611612b0fdc Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Thu, 29 Aug 2024 13:38:52 -0700 Subject: [PATCH 03/20] time checking and minor change. still work in progress. --- include/interface_form.hpp | 6 ++ include/rom_interfaceform.hpp | 6 ++ include/unsteady_ns_solver.hpp | 3 + src/interface_form.cpp | 49 ++++++++++- src/rom_handler.cpp | 3 +- src/rom_interfaceform.cpp | 115 +++++++++++++++++++------- src/steady_ns_solver.cpp | 5 ++ src/unsteady_ns_solver.cpp | 145 +++++++++++++++++++++++++++++++++ 8 files changed, 302 insertions(+), 30 deletions(-) diff --git a/include/interface_form.hpp b/include/interface_form.hpp index 515960b5..ed1cfcbb 100644 --- a/include/interface_form.hpp +++ b/include/interface_form.hpp @@ -15,6 +15,12 @@ namespace mfem // TODO(kevin): inherits or cherry-pick ROMNonlinearForm for hyper-reduction. class InterfaceForm { +private: + static const int Ntimer = 6; + mutable StopWatch *timers[Ntimer]; + mutable int topol_call = 0; + mutable int assemble_call = 0; + protected: int numSub = -1; int skip_zeros = 1; diff --git a/include/rom_interfaceform.hpp b/include/rom_interfaceform.hpp index e091fbd7..7e7a033a 100644 --- a/include/rom_interfaceform.hpp +++ b/include/rom_interfaceform.hpp @@ -16,6 +16,12 @@ namespace mfem class ROMInterfaceForm : public InterfaceForm { +private: + static const int Ntimer = 9; + mutable StopWatch *timers[Ntimer]; + mutable int topol_call = 0; + mutable int assemble_call = 0; + protected: const int numPorts; const int numRefPorts; diff --git a/include/unsteady_ns_solver.hpp b/include/unsteady_ns_solver.hpp index fd3e367d..d7a30ad1 100644 --- a/include/unsteady_ns_solver.hpp +++ b/include/unsteady_ns_solver.hpp @@ -67,6 +67,9 @@ friend class SteadyNSOperator; Array rom_u_offsets; BlockMatrix *rom_mass = NULL; +private: + double times[10]; + public: UnsteadyNSSolver(); diff --git a/src/interface_form.cpp b/src/interface_form.cpp index 4839a40d..845cf19d 100644 --- a/src/interface_form.cpp +++ b/src/interface_form.cpp @@ -12,7 +12,7 @@ namespace mfem InterfaceForm::InterfaceForm( Array &meshes_, Array &fes_, TopologyHandler *topol_) - : meshes(meshes_), fes(fes_), topol_handler(topol_), numSub(meshes_.Size()) + : meshes(meshes_), fes(fes_), topol_handler(topol_), numSub(meshes_.Size()), topol_call(0), assemble_call(0) { assert(fes_.Size() == numSub); @@ -21,11 +21,31 @@ InterfaceForm::InterfaceForm( for (int i = 1; i < numSub + 1; i++) block_offsets[i] = fes[i-1]->GetTrueVSize(); block_offsets.PartialSum(); + + for (int k = 0; k < Ntimer; k++) + timers[k] = new StopWatch; } InterfaceForm::~InterfaceForm() { DeletePointers(fnfi); + + printf("InterfaceForm::InterfaceAddMult\n"); + printf("topol call: %d\n", topol_call); + printf("assemble call: %d\n", assemble_call); + printf("%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\n", + "init", "port-init", "topol", "assemble", "final", "sum", "total"); + double sum = 0.0; + for (int k = 0; k < Ntimer-1; k++) + { + printf("%.3E\t", timers[k]->RealTime()); + sum += timers[k]->RealTime(); + } + printf("%.3E\t", sum); + printf("%.3E\n", timers[Ntimer-1]->RealTime()); + + for (int k = 0; k < Ntimer; k++) + delete timers[k]; } void InterfaceForm::AssembleInterfaceMatrices(Array2D &mats) const @@ -65,6 +85,10 @@ void InterfaceForm::AssembleInterfaceMatrices(Array2D &mats) con void InterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const { +printf("InterfaceForm::InterfaceAddMult\n"); + timers[5]->Start(); + timers[0]->Start(); + x_tmp.Update(const_cast(x), block_offsets); y_tmp.Update(y, block_offsets); @@ -72,8 +96,12 @@ void InterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const Mesh *mesh1, *mesh2; FiniteElementSpace *fes1, *fes2; + timers[0]->Stop(); + for (int p = 0; p < topol_handler->GetNumPorts(); p++) { + timers[1]->Start(); + const PortInfo *pInfo = topol_handler->GetPortInfo(p); midx[0] = pInfo->Mesh1; @@ -86,13 +114,22 @@ void InterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const fes2 = fes[midx[1]]; Array* const interface_infos = topol_handler->GetInterfaceInfos(p); + + timers[1]->Stop(); + AssembleInterfaceVector(mesh1, mesh2, fes1, fes2, interface_infos, x_tmp.GetBlock(midx[0]), x_tmp.GetBlock(midx[1]), y_tmp.GetBlock(midx[0]), y_tmp.GetBlock(midx[1])); } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) + timers[4]->Start(); + for (int i=0; i < y_tmp.NumBlocks(); ++i) y_tmp.GetBlock(i).SyncAliasMemory(y); + + timers[4]->Stop(); + + timers[5]->Stop(); } void InterfaceForm::InterfaceGetGradient(const Vector &x, Array2D &mats) const @@ -237,9 +274,16 @@ void InterfaceForm::AssembleInterfaceVector(Mesh *mesh1, Mesh *mesh2, for (int bn = 0; bn < interface_infos->Size(); bn++) { + timers[2]->Start(); + InterfaceInfo *if_info = &((*interface_infos)[bn]); topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); + topol_call++; + + timers[2]->Stop(); + + timers[3]->Start(); if ((tr1 != NULL) && (tr2 != NULL)) { @@ -258,11 +302,14 @@ void InterfaceForm::AssembleInterfaceVector(Mesh *mesh1, Mesh *mesh2, assert(fnfi[itg]); fnfi[itg]->AssembleInterfaceVector(*fe1, *fe2, *tr1, *tr2, el_x1, el_x2, el_y1, el_y2); + assemble_call += fnfi[itg]->GetIntegrationRule()->GetNPoints(); y1.AddElementVector(vdofs1, el_y1); y2.AddElementVector(vdofs2, el_y2); } } // if ((tr1 != NULL) && (tr2 != NULL)) + + timers[3]->Stop(); } // for (int bn = 0; bn < interface_infos.Size(); bn++) } diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index fb17457b..a4796b4c 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -181,7 +181,8 @@ void ROMHandlerBase::ParseSupremizerInput(Array &num_ref_supreme, Array &meshes_, Array &fes_, Array &comp_fes_, TopologyHandler *topol_) : InterfaceForm(meshes_, fes_, topol_), numPorts(topol_->GetNumPorts()), numRefPorts(topol_->GetNumRefPorts()), - num_comp(topol_->GetNumComponents()), comp_fes(comp_fes_) + num_comp(topol_->GetNumComponents()), comp_fes(comp_fes_), topol_call(0), assemble_call(0) { comp_basis.SetSize(num_comp); comp_basis = NULL; @@ -32,11 +32,31 @@ ROMInterfaceForm::ROMInterfaceForm( else if (nnls_str == "linf") nnls_criterion = CAROM::NNLS_termination::LINF; else mfem_error("ROMNonlinearForm: unknown NNLS criterion!\n"); + + for (int k = 0; k < Ntimer; k++) + timers[k] = new StopWatch; } ROMInterfaceForm::~ROMInterfaceForm() { DeletePointers(fnfi_ref_sample); + + printf("ROMInterfaceForm::InterfaceAddMult\n"); + printf("topol call: %d\n", topol_call); + printf("assemble call: %d\n", assemble_call); + printf("%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\n", + "init", "itg-init", "port-init", "elem-init", "fast-eqp", "slow-eqp", "final", "sum", "total", "topol"); + double sum = 0.0; + for (int k = 0; k < Ntimer-1; k++) + { + printf("%.3E\t", timers[k]->RealTime()); + sum += timers[k]->RealTime(); + } + printf("%.3E\t", sum); + printf("%.3E\n", timers[Ntimer-1]->RealTime()); + + for (int k = 0; k < Ntimer; k++) + delete timers[k]; } void ROMInterfaceForm::SetBasisAtComponent(const int c, DenseMatrix &basis_, const int offset) @@ -73,6 +93,10 @@ void ROMInterfaceForm::UpdateBlockOffsets() void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const { + timers[7]->Start(); + + timers[0]->Start(); + assert(block_offsets.Min() >= 0); assert(x.Size() == block_offsets.Last()); assert(y.Size() == block_offsets.Last()); @@ -96,15 +120,23 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const Array vdofs1, vdofs2; Vector el_x1, el_x2, el_y1, el_y2; + timers[0]->Stop(); + for (int k = 0; k < fnfi.Size(); k++) { + timers[1]->Start(); + assert(fnfi[k]); const IntegrationRule *ir = fnfi[k]->GetIntegrationRule(); assert(ir); // we enforce all integrators to set the IntegrationRule a priori. + timers[1]->Stop(); + for (int p = 0; p < numPorts; p++) { + timers[2]->Start(); + const PortInfo *pInfo = topol_handler->GetPortInfo(p); midx[0] = pInfo->Mesh1; @@ -131,61 +163,88 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const assert(interface_infos); eqp_elem = fnfi_sample[p + k * numPorts]; + + timers[2]->Stop(); + if (!eqp_elem) continue; int prev_itf = -1; for (int i = 0; i < eqp_elem->Size(); i++) { + timers[3]->Start(); + EQPSample *sample = eqp_elem->GetSample(i); int itf = sample->info.el; - InterfaceInfo *if_info = &((*interface_infos)[itf]); - topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); + if (itf != prev_itf) + { + InterfaceInfo *if_info = &((*interface_infos)[itf]); + timers[8]->Start(); + topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); + topol_call++; + assemble_call++; + timers[8]->Stop(); + + if (!precompute) + { + if ((tr1 == NULL) || (tr2 == NULL)) + mfem_error("InterfaceTransformation of the sampled face is NULL,\n" + " indicating that an empty quadrature point is sampled.\n"); + + fes1->GetElementVDofs(tr1->Elem1No, vdofs1); + fes2->GetElementVDofs(tr2->Elem1No, vdofs2); + + if (basis1_offset > 0) + for (int v = 0; v < vdofs1.Size(); v++) + vdofs1[v] += basis1_offset; + if (basis2_offset > 0) + for (int v = 0; v < vdofs2.Size(); v++) + vdofs2[v] += basis2_offset; + + // Both domains will have the adjacent element as Elem1. + fe1 = fes1->GetFE(tr1->Elem1No); + fe2 = fes2->GetFE(tr2->Elem1No); + + MultSubMatrix(*basis1, vdofs1, x1, el_x1); + MultSubMatrix(*basis2, vdofs2, x2, el_x2); + } // if (!precompute) + + prev_itf = itf; + } // if (itf != prev_itf) + const IntegrationPoint &ip = ir->IntPoint(sample->info.qp); + timers[3]->Stop(); + if (precompute) { + timers[4]->Start(); fnfi[k]->AddAssembleVector_Fast(*sample, *tr1, *tr2, x1, x2, y1, y2); + timers[4]->Stop(); continue; } - if (itf != prev_itf) - { - if ((tr1 == NULL) || (tr2 == NULL)) - mfem_error("InterfaceTransformation of the sampled face is NULL,\n" - " indicating that an empty quadrature point is sampled.\n"); - - fes1->GetElementVDofs(tr1->Elem1No, vdofs1); - fes2->GetElementVDofs(tr2->Elem1No, vdofs2); - - if (basis1_offset > 0) - for (int v = 0; v < vdofs1.Size(); v++) - vdofs1[v] += basis1_offset; - if (basis2_offset > 0) - for (int v = 0; v < vdofs2.Size(); v++) - vdofs2[v] += basis2_offset; - - // Both domains will have the adjacent element as Elem1. - fe1 = fes1->GetFE(tr1->Elem1No); - fe2 = fes2->GetFE(tr2->Elem1No); - - MultSubMatrix(*basis1, vdofs1, x1, el_x1); - MultSubMatrix(*basis2, vdofs2, x2, el_x2); - - prev_itf = itf; - } // if (itf != prev_itf) + timers[5]->Start(); fnfi[k]->AssembleQuadratureVector( *fe1, *fe2, *tr1, *tr2, ip, sample->info.qw, el_x1, el_x2, el_y1, el_y2); AddMultTransposeSubMatrix(*basis1, vdofs1, el_y1, y1); AddMultTransposeSubMatrix(*basis2, vdofs2, el_y2, y2); + + timers[5]->Stop(); } // for (int i = 0; i < sample_info->Size(); i++, sample++) } // for (int p = 0; p < numPorts; p++) } // for (int k = 0; k < fnfi.Size(); k++) + timers[6]->Start(); + for (int i=0; i < y_tmp.NumBlocks(); ++i) y_tmp.GetBlock(i).SyncAliasMemory(y); + + timers[6]->Stop(); + + timers[7]->Stop(); } void ROMInterfaceForm::InterfaceGetGradient(const Vector &x, Array2D &mats) const diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index 4418dd2f..3fff79d7 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -1032,6 +1032,8 @@ void SteadyNSSolver::AllocateROMEQPElems() itf_eqp->UpdateBlockOffsets(); itf_eqp->SetPrecomputeMode(precompute); +printf("precompute: %d\n", precompute); +printf("itf_eqp PrecomputeMode: %d\n", itf_eqp->PrecomputeMode()); } void SteadyNSSolver::TrainROMEQPElems(SampleGenerator *sample_generator) @@ -1220,7 +1222,10 @@ void SteadyNSSolver::LoadEQPElems(const std::string &filename) itf_eqp->LoadEQPForIntegrator(0, file_id, "interface_integ0"); if (itf_eqp->PrecomputeMode()) +{ +printf("\n\nI AM PRECOMPUTING!!\n\n"); itf_eqp->PrecomputeCoefficients(); +} } errf = H5Fclose(file_id); diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp index 6cdd7ec6..d1e0982d 100644 --- a/src/unsteady_ns_solver.cpp +++ b/src/unsteady_ns_solver.cpp @@ -88,6 +88,11 @@ bool UnsteadyNSSolver::Solve(SampleGenerator *sample_generator) std::string restart_file, file_fmt; file_fmt = "%s/%s_%08d.h5"; + StopWatch timers; + for (int k = 0; k < 10; k++) + times[k] = 0.0; + timers.Start(); + SetupInitialCondition(initial_step, time); int sample_interval = config.GetOption("sample_generation/time-integration/sample_interval", 0); @@ -99,11 +104,17 @@ bool UnsteadyNSSolver::Solve(SampleGenerator *sample_generator) SaveVisualization(0, time); + timers.Stop(); + times[0] = timers.RealTime(); + timers.Clear(); + double cfl = 0.0; for (int step = initial_step; step < nt; step++) { Step(time, step); + timers.Start(); + cfl = ComputeCFL(dt); SanityCheck(step); if (report_interval && @@ -125,10 +136,29 @@ bool UnsteadyNSSolver::Solve(SampleGenerator *sample_generator) if (sample_generator && sample_interval && ((step+1) > bootstrap) && (((step+1) % sample_interval) == 0)) SaveSnapshots(sample_generator); + + timers.Stop(); + times[9] += timers.RealTime(); + timers.Clear(); } + timers.Start(); + SortBySubdomains(*U_step, *U); + timers.Stop(); + times[9] += timers.RealTime(); + timers.Clear(); + + double total = 0.0; + for (int k = 0; k < 10; k++) + total += times[k]; + + printf("FOM Solve time: %.3e\n", total); + for (int k = 0; k < 10; k++) + printf("%.3e\t", times[k]); + printf("\n"); + return converged; } @@ -208,29 +238,68 @@ void UnsteadyNSSolver::SetupInitialCondition(int &initial_step, double &time) void UnsteadyNSSolver::Step(double &time, int step) { + StopWatch timers; + timers.Start(); + /* set time for forcing/boundary. At this point, time remains at the previous timestep. */ SetTime(time); /* copy velocity */ u1 = U_stepview->GetBlock(0); + timers.Stop(); + times[1] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* evaluate nonlinear advection at previous time step */ Hop->Mult(u1, Cu1); + + timers.Stop(); + times[2] += timers.RealTime(); + timers.Clear(); + timers.Start(); + nl_itf->InterfaceAddMult(u1, Cu1); + timers.Stop(); + times[3] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* Base right-hand side for boundary conditions and forcing */ SortByVariables(*RHS, *RHS_step); + timers.Stop(); + times[4] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* Add nonlinear convection */ RHS_stepview->GetBlock(0).Add(-ab1, Cu1); + timers.Stop(); + times[5] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* Add time derivative term */ // TODO: extend for high order bdf schemes massMat->AddMult(u1, RHS_stepview->GetBlock(0), -bd1 / dt); + timers.Stop(); + times[6] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* Solve for the next step */ mumps->Mult(*RHS_step, *U_step); + timers.Stop(); + times[7] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* remove pressure scalar if all dirichlet bc */ if (!pres_dbc) { @@ -239,6 +308,10 @@ void UnsteadyNSSolver::Step(double &time, int step) U_stepview->GetBlock(1) -= p_const; } + timers.Stop(); + times[8] += timers.RealTime(); + timers.Clear(); + time += dt; } @@ -431,6 +504,12 @@ void UnsteadyNSSolver::SolveROM() { assert(rom_handler->GetOrdering() == ROMOrderBy::VARIABLE); + StopWatch timers; + for (int k = 0; k < 10; k++) + times[k] = 0.0; + + timers.Start(); + int initial_step = 0; double time = 0.0; @@ -482,31 +561,78 @@ void UnsteadyNSSolver::SolveROM() } } + timers.Stop(); + times[0] = timers.RealTime(); + timers.Clear(); + for (int step = initial_step; step < nt; step++) { + timers.Start(); + /* set time for forcing/boundary. At this point, time remains at the previous timestep. */ SetTime(time); + timers.Stop(); + times[0] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* copy velocity */ u1 = rsol_view->GetBlock(0); + timers.Stop(); + times[1] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* evaluate nonlinear advection at previous time step */ Hop->Mult(u1, Cu1); + + timers.Stop(); + times[2] += timers.RealTime(); + timers.Clear(); + timers.Start(); + itf_eqp->InterfaceAddMult(u1, Cu1); + timers.Stop(); + times[3] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* Base right-hand side for boundary conditions and forcing */ reduced_rhs = *(rom_handler->GetReducedRHS()); + timers.Stop(); + times[4] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* Add nonlinear convection */ rrhs_view->GetBlock(0).Add(-ab1, Cu1); + timers.Stop(); + times[5] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* Add time derivative term */ // TODO: extend for high order bdf schemes rom_mass->AddMult(u1, rrhs_view->GetBlock(0), -bd1 / dt); + timers.Stop(); + times[6] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* Solve for the next step */ rom_handler->Solve(reduced_rhs, *reduced_sol); + timers.Stop(); + times[7] += timers.RealTime(); + timers.Clear(); + timers.Start(); + /* remove pressure scalar if all dirichlet bc */ if (!pres_dbc) { @@ -515,11 +641,30 @@ void UnsteadyNSSolver::SolveROM() rsol_view->GetBlock(1).Add(-p_const, rom_ones); } + timers.Stop(); + times[8] += timers.RealTime(); + timers.Clear(); + time += dt; } + timers.Start(); + rom_handler->LiftUpGlobal(*reduced_sol, *U); + timers.Stop(); + times[9] += timers.RealTime(); + timers.Clear(); + + double total = 0.0; + for (int k = 0; k < 10; k++) + total += times[k]; + + printf("ROM Solve time: %.3e\n", total); + for (int k = 0; k < 10; k++) + printf("%.3e\t", times[k]); + printf("\n"); + delete rsol_view; delete reduced_sol; return; From c78649ae2bc10bfbd27910d5dc7fcde1e1d331f4 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Thu, 29 Aug 2024 13:46:00 -0700 Subject: [PATCH 04/20] bug fix on SteadyNSSolver::LoadEQPElems --- src/steady_ns_solver.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index 3fff79d7..55ef5e89 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -1184,7 +1184,7 @@ void SteadyNSSolver::LoadEQPElems(const std::string &filename) int num_comp; hdf5_utils::ReadAttribute(grp_id, "number_of_components", num_comp); assert(num_comp >= topol_handler->GetNumComponents()); - assert(comp_eqps.Size() == num_comp); + assert(comp_eqps.Size() == topol_handler->GetNumComponents()); std::string dset_name; for (int c = 0; c < topol_handler->GetNumComponents(); c++) From 348fd1b1161cc1d8415cdb6218c837d0b7aa79ea Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Mon, 13 Jan 2025 13:58:48 -0800 Subject: [PATCH 05/20] default option for number of supremizer --- src/rom_handler.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index a4796b4c..d45ebf3b 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -175,14 +175,17 @@ void ROMHandlerBase::ParseSupremizerInput(Array &num_ref_supreme, Array("basis/number_of_supremizer", -1); BasisTag basis_tag; for (int b = 0; b < num_rom_comp; b++) { basis_tag = GetBasisTagForComponent(b, topol_handler, fom_var_names[1]); - // By default the same number of velocity basis. - num_ref_supreme[b] = config.LookUpFromDict("name", basis_tag.print(), "number_of_supremizer", num_ref_basis[b * num_var], basis_list); + // By default the same number of pressure basis. + const int nsup0 = (nsup_default > 0) ? nsup_default : num_ref_basis[b * num_var] + 1; + + num_ref_supreme[b] = config.LookUpFromDict("name", basis_tag.print(), "number_of_supremizer", nsup0, basis_list); } for (int m = 0; m < numSub; m++) From d59e498f244eba2b97e4b4bb9f05687d2b310d0f Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Tue, 14 Jan 2025 17:37:39 -0800 Subject: [PATCH 06/20] fixed bug for supremizer basis size --- src/rom_handler.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index d45ebf3b..7700ad15 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -183,7 +183,7 @@ void ROMHandlerBase::ParseSupremizerInput(Array &num_ref_supreme, Array 0) ? nsup_default : num_ref_basis[b * num_var] + 1; + const int nsup0 = (nsup_default > 0) ? nsup_default : num_ref_basis[b * num_var + 1]; num_ref_supreme[b] = config.LookUpFromDict("name", basis_tag.print(), "number_of_supremizer", nsup0, basis_list); } From 8f744131d90f4082ba2d8dcfef2cf12dab222642 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Tue, 14 Jan 2025 17:37:57 -0800 Subject: [PATCH 07/20] option for interface eqp tolerance --- src/steady_ns_solver.cpp | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index 55ef5e89..fe314ce0 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -1062,6 +1062,8 @@ void SteadyNSSolver::TrainROMEQPElems(SampleGenerator *sample_generator) if (oper_type != OperType::LF) return; + eqp_tol = config.GetOption("model_reduction/eqp/interface/relative_tolerance", eqp_tol); + /* EQP NNLS for interface ROM, for each reference port */ for (int p = 0; p < topol_handler->GetNumRefPorts(); p++) { From a5e78ac3b445d980418b50742f9eed04a75f6d76 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Tue, 14 Jan 2025 17:38:22 -0800 Subject: [PATCH 08/20] TimeProfiler class --- include/etc.hpp | 77 +++++++++++++++++++++++++++++++++++ include/rom_interfaceform.hpp | 2 + src/rom_interfaceform.cpp | 6 ++- 3 files changed, 84 insertions(+), 1 deletion(-) diff --git a/include/etc.hpp b/include/etc.hpp index f144f62a..bd5d19c0 100644 --- a/include/etc.hpp +++ b/include/etc.hpp @@ -38,4 +38,81 @@ std::string string_format( const std::string& format, Args ... args ) return std::string( buf.get(), buf.get() + size - 1 ); // We don't want the '\0' inside } +class TimeProfiler +{ +private: + Array timers; + + Array calls; + std::vector names; + std::unordered_map indices; + + const MPI_Comm comm; + int rank; + +public: + TimeProfiler(MPI_Comm comm_=MPI_COMM_WORLD) + : comm(comm_), timers(0) + { + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + calls.SetSize(0); + names.clear(); + indices.clear(); + } + + virtual ~TimeProfiler() + { + DeletePointers(timers); + } + + void Start(const std::string &name) + { + if (!indices.count(name)) + { + indices[name] = timers.Size(); + timers.Append(new StopWatch); + names.push_back(name); + calls.Append(0); + } + + assert(indices.count(name)); + int idx = indices[name]; + timers[idx]->Start(); + } + + void Stop(const std::string &name) + { + assert(indices.count(name)); + int idx = indices[name]; + timers[idx]->Stop(); + calls[idx] += 1; + } + + void Print(const std::string &title) + { + int nfunc = timers.Size(); + Array times(nfunc); + for (int k = 0; k < nfunc; k++) + times[k] = timers[k]->RealTime(); + + MPI_Reduce(MPI_IN_PLACE, times.GetData(), nfunc, MPI_DOUBLE, MPI_SUM, 0, comm); + MPI_Reduce(MPI_IN_PLACE, calls.GetData(), nfunc, MPI_INT, MPI_SUM, 0, comm); + + if (rank == 0) + { + printf((title + "\n").c_str()); + + std::string line = std::string(100, '='); + line += "\n"; + printf(line.c_str()); + printf("%20s\t%20s\t%20s\t%20s\n", "Function", "Total time", "Calls", "Time per call"); + for (int k = 0; k < nfunc; k++) + { + printf("%20s\t%20.5e\t%20d\t%20.5e\n", names[k].c_str(), times[k], calls[k], times[k] / calls[k]); + } + printf(line.c_str()); + } + } +}; + #endif diff --git a/include/rom_interfaceform.hpp b/include/rom_interfaceform.hpp index 7e7a033a..ae9963c9 100644 --- a/include/rom_interfaceform.hpp +++ b/include/rom_interfaceform.hpp @@ -10,6 +10,7 @@ #include "hdf5_utils.hpp" #include "hyperreduction_integ.hpp" #include "linalg/NNLS.h" +#include "etc.hpp" namespace mfem { @@ -21,6 +22,7 @@ class ROMInterfaceForm : public InterfaceForm mutable StopWatch *timers[Ntimer]; mutable int topol_call = 0; mutable int assemble_call = 0; + mutable TimeProfiler timer; protected: const int numPorts; diff --git a/src/rom_interfaceform.cpp b/src/rom_interfaceform.cpp index 79ec2535..aaf1c027 100644 --- a/src/rom_interfaceform.cpp +++ b/src/rom_interfaceform.cpp @@ -14,7 +14,7 @@ namespace mfem ROMInterfaceForm::ROMInterfaceForm( Array &meshes_, Array &fes_, Array &comp_fes_, TopologyHandler *topol_) : InterfaceForm(meshes_, fes_, topol_), numPorts(topol_->GetNumPorts()), numRefPorts(topol_->GetNumRefPorts()), - num_comp(topol_->GetNumComponents()), comp_fes(comp_fes_), topol_call(0), assemble_call(0) + num_comp(topol_->GetNumComponents()), comp_fes(comp_fes_), topol_call(0), assemble_call(0), timer() { comp_basis.SetSize(num_comp); comp_basis = NULL; @@ -41,6 +41,8 @@ ROMInterfaceForm::~ROMInterfaceForm() { DeletePointers(fnfi_ref_sample); + timer.Print("ROMInterfaceForm::InterfaceAddMult"); + printf("ROMInterfaceForm::InterfaceAddMult\n"); printf("topol call: %d\n", topol_call); printf("assemble call: %d\n", assemble_call); @@ -96,6 +98,7 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const timers[7]->Start(); timers[0]->Start(); + timer.Start("init"); assert(block_offsets.Min() >= 0); assert(x.Size() == block_offsets.Last()); @@ -121,6 +124,7 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const Vector el_x1, el_x2, el_y1, el_y2; timers[0]->Stop(); + timer.Stop("init"); for (int k = 0; k < fnfi.Size(); k++) { From 8be9eaae396842fd9f5ccb3cc7a7221a6d8abd12 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Tue, 14 Jan 2025 17:56:29 -0800 Subject: [PATCH 09/20] used time profile in rom interface form --- include/etc.hpp | 20 ++++++++++-- include/rom_interfaceform.hpp | 4 --- src/rom_interfaceform.cpp | 58 ++++++++++------------------------- 3 files changed, 35 insertions(+), 47 deletions(-) diff --git a/include/etc.hpp b/include/etc.hpp index bd5d19c0..6e32d868 100644 --- a/include/etc.hpp +++ b/include/etc.hpp @@ -44,6 +44,7 @@ class TimeProfiler Array timers; Array calls; + Array starts; std::vector names; std::unordered_map indices; @@ -52,7 +53,7 @@ class TimeProfiler public: TimeProfiler(MPI_Comm comm_=MPI_COMM_WORLD) - : comm(comm_), timers(0) + : comm(comm_), timers(0), starts(0) { MPI_Comm_rank(MPI_COMM_WORLD, &rank); calls.SetSize(0); @@ -73,27 +74,40 @@ class TimeProfiler timers.Append(new StopWatch); names.push_back(name); calls.Append(0); + starts.Append(false); } assert(indices.count(name)); + int idx = indices[name]; + assert(!starts[idx]); + timers[idx]->Start(); + starts[idx] = true; } void Stop(const std::string &name) { assert(indices.count(name)); + int idx = indices[name]; + assert(starts[idx]); + timers[idx]->Stop(); calls[idx] += 1; + starts[idx] = false; } - void Print(const std::string &title) + void Print(const std::string &title, const bool compute_sum=false) { int nfunc = timers.Size(); + for (int k = 0; k < nfunc; k++) + assert(!starts[k]); + Array times(nfunc); for (int k = 0; k < nfunc; k++) times[k] = timers[k]->RealTime(); + double total = times.Sum(); MPI_Reduce(MPI_IN_PLACE, times.GetData(), nfunc, MPI_DOUBLE, MPI_SUM, 0, comm); MPI_Reduce(MPI_IN_PLACE, calls.GetData(), nfunc, MPI_INT, MPI_SUM, 0, comm); @@ -110,6 +124,8 @@ class TimeProfiler { printf("%20s\t%20.5e\t%20d\t%20.5e\n", names[k].c_str(), times[k], calls[k], times[k] / calls[k]); } + if (compute_sum) + printf("%20s\t%20.5e\t%20d\t%20.5e\n", "Total time", total, 1, total); printf(line.c_str()); } } diff --git a/include/rom_interfaceform.hpp b/include/rom_interfaceform.hpp index ae9963c9..6c5d4fef 100644 --- a/include/rom_interfaceform.hpp +++ b/include/rom_interfaceform.hpp @@ -18,10 +18,6 @@ namespace mfem class ROMInterfaceForm : public InterfaceForm { private: - static const int Ntimer = 9; - mutable StopWatch *timers[Ntimer]; - mutable int topol_call = 0; - mutable int assemble_call = 0; mutable TimeProfiler timer; protected: diff --git a/src/rom_interfaceform.cpp b/src/rom_interfaceform.cpp index aaf1c027..3d28d9e6 100644 --- a/src/rom_interfaceform.cpp +++ b/src/rom_interfaceform.cpp @@ -14,7 +14,7 @@ namespace mfem ROMInterfaceForm::ROMInterfaceForm( Array &meshes_, Array &fes_, Array &comp_fes_, TopologyHandler *topol_) : InterfaceForm(meshes_, fes_, topol_), numPorts(topol_->GetNumPorts()), numRefPorts(topol_->GetNumRefPorts()), - num_comp(topol_->GetNumComponents()), comp_fes(comp_fes_), topol_call(0), assemble_call(0), timer() + num_comp(topol_->GetNumComponents()), comp_fes(comp_fes_), timer() { comp_basis.SetSize(num_comp); comp_basis = NULL; @@ -32,9 +32,6 @@ ROMInterfaceForm::ROMInterfaceForm( else if (nnls_str == "linf") nnls_criterion = CAROM::NNLS_termination::LINF; else mfem_error("ROMNonlinearForm: unknown NNLS criterion!\n"); - - for (int k = 0; k < Ntimer; k++) - timers[k] = new StopWatch; } ROMInterfaceForm::~ROMInterfaceForm() @@ -42,23 +39,6 @@ ROMInterfaceForm::~ROMInterfaceForm() DeletePointers(fnfi_ref_sample); timer.Print("ROMInterfaceForm::InterfaceAddMult"); - - printf("ROMInterfaceForm::InterfaceAddMult\n"); - printf("topol call: %d\n", topol_call); - printf("assemble call: %d\n", assemble_call); - printf("%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\n", - "init", "itg-init", "port-init", "elem-init", "fast-eqp", "slow-eqp", "final", "sum", "total", "topol"); - double sum = 0.0; - for (int k = 0; k < Ntimer-1; k++) - { - printf("%.3E\t", timers[k]->RealTime()); - sum += timers[k]->RealTime(); - } - printf("%.3E\t", sum); - printf("%.3E\n", timers[Ntimer-1]->RealTime()); - - for (int k = 0; k < Ntimer; k++) - delete timers[k]; } void ROMInterfaceForm::SetBasisAtComponent(const int c, DenseMatrix &basis_, const int offset) @@ -95,9 +75,8 @@ void ROMInterfaceForm::UpdateBlockOffsets() void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const { - timers[7]->Start(); + timer.Start("Total"); - timers[0]->Start(); timer.Start("init"); assert(block_offsets.Min() >= 0); @@ -123,23 +102,22 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const Array vdofs1, vdofs2; Vector el_x1, el_x2, el_y1, el_y2; - timers[0]->Stop(); timer.Stop("init"); for (int k = 0; k < fnfi.Size(); k++) { - timers[1]->Start(); + timer.Start("itg-init"); assert(fnfi[k]); const IntegrationRule *ir = fnfi[k]->GetIntegrationRule(); assert(ir); // we enforce all integrators to set the IntegrationRule a priori. - timers[1]->Stop(); + timer.Stop("itg-init"); for (int p = 0; p < numPorts; p++) { - timers[2]->Start(); + timer.Start("port-init"); const PortInfo *pInfo = topol_handler->GetPortInfo(p); @@ -168,14 +146,14 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const eqp_elem = fnfi_sample[p + k * numPorts]; - timers[2]->Stop(); + timer.Stop("port-init"); if (!eqp_elem) continue; int prev_itf = -1; for (int i = 0; i < eqp_elem->Size(); i++) { - timers[3]->Start(); + timer.Start("elem-init"); EQPSample *sample = eqp_elem->GetSample(i); @@ -183,11 +161,9 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const if (itf != prev_itf) { InterfaceInfo *if_info = &((*interface_infos)[itf]); - timers[8]->Start(); + timer.Start("topol"); topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); - topol_call++; - assemble_call++; - timers[8]->Stop(); + timer.Stop("topol"); if (!precompute) { @@ -218,17 +194,17 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const const IntegrationPoint &ip = ir->IntPoint(sample->info.qp); - timers[3]->Stop(); + timer.Stop("elem-init"); if (precompute) { - timers[4]->Start(); + timer.Start("fast-eqp"); fnfi[k]->AddAssembleVector_Fast(*sample, *tr1, *tr2, x1, x2, y1, y2); - timers[4]->Stop(); + timer.Stop("fast-eqp"); continue; } - timers[5]->Start(); + timer.Start("slow-eqp"); fnfi[k]->AssembleQuadratureVector( *fe1, *fe2, *tr1, *tr2, ip, sample->info.qw, el_x1, el_x2, el_y1, el_y2); @@ -236,19 +212,19 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const AddMultTransposeSubMatrix(*basis1, vdofs1, el_y1, y1); AddMultTransposeSubMatrix(*basis2, vdofs2, el_y2, y2); - timers[5]->Stop(); + timer.Stop("slow-eqp"); } // for (int i = 0; i < sample_info->Size(); i++, sample++) } // for (int p = 0; p < numPorts; p++) } // for (int k = 0; k < fnfi.Size(); k++) - timers[6]->Start(); + timer.Start("final"); for (int i=0; i < y_tmp.NumBlocks(); ++i) y_tmp.GetBlock(i).SyncAliasMemory(y); - timers[6]->Stop(); + timer.Stop("final"); - timers[7]->Stop(); + timer.Stop("Total"); } void ROMInterfaceForm::InterfaceGetGradient(const Vector &x, Array2D &mats) const From e36deb0fdf64d3e2475e789665b0238c45fb5a1b Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Tue, 14 Jan 2025 18:11:47 -0800 Subject: [PATCH 10/20] used time profiler in interface form --- include/interface_form.hpp | 5 +--- src/interface_form.cpp | 52 ++++++++++++-------------------------- 2 files changed, 17 insertions(+), 40 deletions(-) diff --git a/include/interface_form.hpp b/include/interface_form.hpp index ed1cfcbb..9a2185b9 100644 --- a/include/interface_form.hpp +++ b/include/interface_form.hpp @@ -16,10 +16,7 @@ namespace mfem class InterfaceForm { private: - static const int Ntimer = 6; - mutable StopWatch *timers[Ntimer]; - mutable int topol_call = 0; - mutable int assemble_call = 0; + mutable TimeProfiler timer; protected: int numSub = -1; diff --git a/src/interface_form.cpp b/src/interface_form.cpp index 845cf19d..a2cf7c24 100644 --- a/src/interface_form.cpp +++ b/src/interface_form.cpp @@ -12,7 +12,7 @@ namespace mfem InterfaceForm::InterfaceForm( Array &meshes_, Array &fes_, TopologyHandler *topol_) - : meshes(meshes_), fes(fes_), topol_handler(topol_), numSub(meshes_.Size()), topol_call(0), assemble_call(0) + : meshes(meshes_), fes(fes_), topol_handler(topol_), numSub(meshes_.Size()), timer() { assert(fes_.Size() == numSub); @@ -21,31 +21,13 @@ InterfaceForm::InterfaceForm( for (int i = 1; i < numSub + 1; i++) block_offsets[i] = fes[i-1]->GetTrueVSize(); block_offsets.PartialSum(); - - for (int k = 0; k < Ntimer; k++) - timers[k] = new StopWatch; } InterfaceForm::~InterfaceForm() { DeletePointers(fnfi); - printf("InterfaceForm::InterfaceAddMult\n"); - printf("topol call: %d\n", topol_call); - printf("assemble call: %d\n", assemble_call); - printf("%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\n", - "init", "port-init", "topol", "assemble", "final", "sum", "total"); - double sum = 0.0; - for (int k = 0; k < Ntimer-1; k++) - { - printf("%.3E\t", timers[k]->RealTime()); - sum += timers[k]->RealTime(); - } - printf("%.3E\t", sum); - printf("%.3E\n", timers[Ntimer-1]->RealTime()); - - for (int k = 0; k < Ntimer; k++) - delete timers[k]; + timer.Print("InterfaceForm::InterfaceAddMult"); } void InterfaceForm::AssembleInterfaceMatrices(Array2D &mats) const @@ -85,9 +67,8 @@ void InterfaceForm::AssembleInterfaceMatrices(Array2D &mats) con void InterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const { -printf("InterfaceForm::InterfaceAddMult\n"); - timers[5]->Start(); - timers[0]->Start(); + timer.Start("Total"); + timer.Start("init"); x_tmp.Update(const_cast(x), block_offsets); y_tmp.Update(y, block_offsets); @@ -96,11 +77,11 @@ printf("InterfaceForm::InterfaceAddMult\n"); Mesh *mesh1, *mesh2; FiniteElementSpace *fes1, *fes2; - timers[0]->Stop(); + timer.Stop("init"); for (int p = 0; p < topol_handler->GetNumPorts(); p++) { - timers[1]->Start(); + timer.Start("port-init"); const PortInfo *pInfo = topol_handler->GetPortInfo(p); @@ -115,21 +96,20 @@ printf("InterfaceForm::InterfaceAddMult\n"); Array* const interface_infos = topol_handler->GetInterfaceInfos(p); - timers[1]->Stop(); + timer.Stop("port-init"); AssembleInterfaceVector(mesh1, mesh2, fes1, fes2, interface_infos, x_tmp.GetBlock(midx[0]), x_tmp.GetBlock(midx[1]), y_tmp.GetBlock(midx[0]), y_tmp.GetBlock(midx[1])); } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) - timers[4]->Start(); + timer.Start("final"); for (int i=0; i < y_tmp.NumBlocks(); ++i) y_tmp.GetBlock(i).SyncAliasMemory(y); - timers[4]->Stop(); - - timers[5]->Stop(); + timer.Stop("final"); + timer.Stop("Total"); } void InterfaceForm::InterfaceGetGradient(const Vector &x, Array2D &mats) const @@ -274,16 +254,15 @@ void InterfaceForm::AssembleInterfaceVector(Mesh *mesh1, Mesh *mesh2, for (int bn = 0; bn < interface_infos->Size(); bn++) { - timers[2]->Start(); + timer.Start("topol"); InterfaceInfo *if_info = &((*interface_infos)[bn]); topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); - topol_call++; - timers[2]->Stop(); + timer.Stop("topol"); - timers[3]->Start(); + timer.Start("assemble"); if ((tr1 != NULL) && (tr2 != NULL)) { @@ -301,15 +280,16 @@ void InterfaceForm::AssembleInterfaceVector(Mesh *mesh1, Mesh *mesh2, { assert(fnfi[itg]); + timer.Start("assemble-itf-vec"); fnfi[itg]->AssembleInterfaceVector(*fe1, *fe2, *tr1, *tr2, el_x1, el_x2, el_y1, el_y2); - assemble_call += fnfi[itg]->GetIntegrationRule()->GetNPoints(); + timer.Stop("assemble-itf-vec"); y1.AddElementVector(vdofs1, el_y1); y2.AddElementVector(vdofs2, el_y2); } } // if ((tr1 != NULL) && (tr2 != NULL)) - timers[3]->Stop(); + timer.Stop("assemble"); } // for (int bn = 0; bn < interface_infos.Size(); bn++) } From b79a16d097cc8726fcced636ff3379905747eef4 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Fri, 17 Jan 2025 11:55:30 -0800 Subject: [PATCH 11/20] ROMNonlinearForm::GetGradient uses TimeProfiler. --- include/etc.hpp | 12 ++++++++-- include/rom_nonlinearform.hpp | 6 +++-- src/rom_nonlinearform.cpp | 43 +++++++++++++---------------------- 3 files changed, 30 insertions(+), 31 deletions(-) diff --git a/include/etc.hpp b/include/etc.hpp index 6e32d868..5b05d822 100644 --- a/include/etc.hpp +++ b/include/etc.hpp @@ -80,7 +80,11 @@ class TimeProfiler assert(indices.count(name)); int idx = indices[name]; - assert(!starts[idx]); + if (starts[idx]) + { + printf("TimeProfiler: %s is already started!\n", name.c_str()); + assert(!starts[idx]); + } timers[idx]->Start(); starts[idx] = true; @@ -91,7 +95,11 @@ class TimeProfiler assert(indices.count(name)); int idx = indices[name]; - assert(starts[idx]); + if (!starts[idx]) + { + printf("TimeProfiler: %s is not started!\n", name.c_str()); + assert(starts[idx]); + } timers[idx]->Stop(); calls[idx] += 1; diff --git a/include/rom_nonlinearform.hpp b/include/rom_nonlinearform.hpp index 0ec1c360..20db9d03 100644 --- a/include/rom_nonlinearform.hpp +++ b/include/rom_nonlinearform.hpp @@ -16,9 +16,11 @@ namespace mfem class ROMNonlinearForm : public NonlinearForm { +// private: +// static const int Nt = 6; +// mutable StopWatch *jac_timers[Nt]; private: - static const int Nt = 6; - mutable StopWatch *jac_timers[Nt]; + mutable TimeProfiler timer; protected: /// ROM basis for projection. diff --git a/src/rom_nonlinearform.cpp b/src/rom_nonlinearform.cpp index 61651773..511fba9a 100644 --- a/src/rom_nonlinearform.cpp +++ b/src/rom_nonlinearform.cpp @@ -12,7 +12,7 @@ namespace mfem { ROMNonlinearForm::ROMNonlinearForm(const int num_basis, FiniteElementSpace *f, const bool reference_) - : NonlinearForm(f), reference(reference_) + : NonlinearForm(f), reference(reference_), timer() { height = width = num_basis; @@ -21,8 +21,6 @@ ROMNonlinearForm::ROMNonlinearForm(const int num_basis, FiniteElementSpace *f, c else if (nnls_str == "linf") nnls_criterion = CAROM::NNLS_termination::LINF; else mfem_error("ROMNonlinearForm: unknown NNLS criterion!\n"); - - for (int k = 0; k < Nt; k++) jac_timers[k] = new StopWatch; } ROMNonlinearForm::~ROMNonlinearForm() @@ -44,18 +42,7 @@ ROMNonlinearForm::~ROMNonlinearForm() DeletePointers(bfnfi_sample); } - printf("%10s\t%10s\t%10s\t%10s\t%10s\t%10s\t%10s\n", - "init", "sample-load", "elem-load", "assem-grad", "others", "sum", "total"); - double sum = 0.0; - for (int k = 0; k < Nt-1; k++) - { - printf("%.3E\t", jac_timers[k]->RealTime()); - sum += jac_timers[k]->RealTime(); - } - printf("%.3E\t", sum); - printf("%.3E\n", jac_timers[Nt-1]->RealTime()); - - for (int k = 0; k < Nt; k++) delete jac_timers[k]; + timer.Print("ROMNonlinearForm::GetGradient"); } void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const @@ -312,9 +299,10 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const Operator& ROMNonlinearForm::GetGradient(const Vector &x) const { - jac_timers[5]->Start(); + timer.Start("Total"); - jac_timers[0]->Start(); + // jac_timers[0]->Start(); + timer.Start("init"); assert(x.Size() == width); // if (ext) // { @@ -356,13 +344,13 @@ Operator& ROMNonlinearForm::GetGradient(const Vector &x) const { *Grad = 0.0; } - jac_timers[0]->Stop(); + timer.Stop("init"); if (dnfi.Size()) { for (int k = 0; k < dnfi.Size(); k++) { - jac_timers[1]->Start(); + timer.Start("sample-load"); const IntegrationRule *ir = dnfi[k]->GetIntegrationRule(); assert(ir); // we enforce all integrators to set the IntegrationRule a priori. @@ -370,19 +358,20 @@ Operator& ROMNonlinearForm::GetGradient(const Vector &x) const assert(eqp_elem); int prev_el = -1; - jac_timers[1]->Stop(); + timer.Stop("sample-load"); for (int i = 0; i < eqp_elem->Size(); i++) { - jac_timers[2]->Start(); + timer.Start("elem-load"); EQPSample *sample = eqp_elem->GetSample(i); int el = sample->info.el; T = fes->GetElementTransformation(el); - jac_timers[2]->Stop(); + timer.Stop("elem-load"); - jac_timers[3]->Start(); + timer.Start("assem-grad"); if (precompute) { dnfi[k]->AddAssembleGrad_Fast(*sample, *T, x, *Grad); + timer.Stop("assem-grad"); continue; } @@ -403,12 +392,12 @@ Operator& ROMNonlinearForm::GetGradient(const Vector &x) const AddSubMatrixRtAP(*basis, vdofs, elmat, *basis, vdofs, *Grad); // Grad->AddSubMatrix(rom_vdofs, rom_vdofs, quadmat, skip_zeros); - jac_timers[3]->Stop(); + timer.Stop("assem-grad"); } // for (int i = 0; i < el_samples->Size(); i++) } // for (int k = 0; k < dnfi.Size(); k++) } // if (dnfi.Size()) - jac_timers[4]->Start(); + timer.Start("other-integ"); if (fnfi.Size()) { FaceElementTransformations *tr = NULL; @@ -566,8 +555,8 @@ Operator& ROMNonlinearForm::GetGradient(const Vector &x) const DenseMatrix *mGrad = Grad; - jac_timers[4]->Stop(); - jac_timers[5]->Stop(); + timer.Stop("other-integ"); + timer.Stop("Total"); // TODO(kevin): will need to consider how we should address this case. // do we really need prolongation when we lift up from ROM? // we might need a similar operation on the ROM basis DenseMatrix. From f81d194d6742780cfd4d076f35b7a90910e3e3d7 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Fri, 17 Jan 2025 16:33:01 -0800 Subject: [PATCH 12/20] time profiling on rom nonlinear form --- include/steady_ns_solver.hpp | 9 ++++- src/rom_interfaceform.cpp | 49 +++++++++++----------- src/rom_nonlinearform.cpp | 78 +++++++++++++++++++++++++++++------- src/steady_ns_solver.cpp | 22 +++++++++- 4 files changed, 115 insertions(+), 43 deletions(-) diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 13ee520d..2531ebe3 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -8,6 +8,7 @@ #include "stokes_solver.hpp" #include "rom_nonlinearform.hpp" #include "rom_interfaceform.hpp" +#include "etc.hpp" // By convention we only use mfem namespace as default, not CAROM. using namespace mfem; @@ -99,6 +100,9 @@ class SteadyNSTensorROM : public SteadyNSROM class SteadyNSEQPROM : public SteadyNSROM { +private: + mutable TimeProfiler timer; + protected: Array hs; // not owned by SteadyNSEQPROM. ROMInterfaceForm *itf = NULL; // not owned by SteadyNSEQPROM. @@ -111,7 +115,10 @@ class SteadyNSEQPROM : public SteadyNSROM SteadyNSEQPROM(ROMHandlerBase *rom_handler, Array &hs_, ROMInterfaceForm *itf_, const bool direct_solve_=true); - virtual ~SteadyNSEQPROM() {} + virtual ~SteadyNSEQPROM() + { + // timer.Print("SteadyNSEQPROM"); + } virtual void Mult(const Vector &x, Vector &y) const; virtual Operator &GetGradient(const Vector &x) const; diff --git a/src/rom_interfaceform.cpp b/src/rom_interfaceform.cpp index 3d28d9e6..de46149c 100644 --- a/src/rom_interfaceform.cpp +++ b/src/rom_interfaceform.cpp @@ -165,35 +165,9 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); timer.Stop("topol"); - if (!precompute) - { - if ((tr1 == NULL) || (tr2 == NULL)) - mfem_error("InterfaceTransformation of the sampled face is NULL,\n" - " indicating that an empty quadrature point is sampled.\n"); - - fes1->GetElementVDofs(tr1->Elem1No, vdofs1); - fes2->GetElementVDofs(tr2->Elem1No, vdofs2); - - if (basis1_offset > 0) - for (int v = 0; v < vdofs1.Size(); v++) - vdofs1[v] += basis1_offset; - if (basis2_offset > 0) - for (int v = 0; v < vdofs2.Size(); v++) - vdofs2[v] += basis2_offset; - - // Both domains will have the adjacent element as Elem1. - fe1 = fes1->GetFE(tr1->Elem1No); - fe2 = fes2->GetFE(tr2->Elem1No); - - MultSubMatrix(*basis1, vdofs1, x1, el_x1); - MultSubMatrix(*basis2, vdofs2, x2, el_x2); - } // if (!precompute) - prev_itf = itf; } // if (itf != prev_itf) - const IntegrationPoint &ip = ir->IntPoint(sample->info.qp); - timer.Stop("elem-init"); if (precompute) @@ -206,6 +180,29 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const timer.Start("slow-eqp"); + if ((tr1 == NULL) || (tr2 == NULL)) + mfem_error("InterfaceTransformation of the sampled face is NULL,\n" + " indicating that an empty quadrature point is sampled.\n"); + + fes1->GetElementVDofs(tr1->Elem1No, vdofs1); + fes2->GetElementVDofs(tr2->Elem1No, vdofs2); + + if (basis1_offset > 0) + for (int v = 0; v < vdofs1.Size(); v++) + vdofs1[v] += basis1_offset; + if (basis2_offset > 0) + for (int v = 0; v < vdofs2.Size(); v++) + vdofs2[v] += basis2_offset; + + // Both domains will have the adjacent element as Elem1. + fe1 = fes1->GetFE(tr1->Elem1No); + fe2 = fes2->GetFE(tr2->Elem1No); + + MultSubMatrix(*basis1, vdofs1, x1, el_x1); + MultSubMatrix(*basis2, vdofs2, x2, el_x2); + + const IntegrationPoint &ip = ir->IntPoint(sample->info.qp); + fnfi[k]->AssembleQuadratureVector( *fe1, *fe2, *tr1, *tr2, ip, sample->info.qw, el_x1, el_x2, el_y1, el_y2); diff --git a/src/rom_nonlinearform.cpp b/src/rom_nonlinearform.cpp index 511fba9a..391e1698 100644 --- a/src/rom_nonlinearform.cpp +++ b/src/rom_nonlinearform.cpp @@ -42,11 +42,15 @@ ROMNonlinearForm::~ROMNonlinearForm() DeletePointers(bfnfi_sample); } - timer.Print("ROMNonlinearForm::GetGradient"); + timer.Print("ROMNonlinearForm"); } void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const { + timer.Start("Mult"); + + timer.Start("Mult/init"); + assert(x.Size() == width); assert(y.Size() == height); @@ -90,29 +94,43 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const // py = 0.0; y = 0.0; + timer.Stop("Mult/init"); + if (dnfi.Size()) { for (int k = 0; k < dnfi.Size(); k++) { + timer.Start("Mult/dnfi-init"); + const IntegrationRule *ir = dnfi[k]->GetIntegrationRule(); assert(ir); // we enforce all integrators to set the IntegrationRule a priori. EQPElement *eqp_elem = dnfi_sample[k]; assert(eqp_elem); + timer.Stop("Mult/dnfi-init"); + int prev_el = -1; for (int i = 0; i < eqp_elem->Size(); i++) { + timer.Start("Mult/dnfi-elem-init"); + EQPSample *sample = eqp_elem->GetSample(i); int el = sample->info.el; T = fes->GetElementTransformation(el); + timer.Stop("Mult/dnfi-elem-init"); + if (precompute) { + timer.Start("Mult/dnfi-fast-eqp"); dnfi[k]->AddAssembleVector_Fast(*sample, *T, x, y); + timer.Stop("Mult/dnfi-fast-eqp"); continue; } + timer.Start("Mult/dnfi-slow-eqp"); + const IntegrationPoint &ip = ir->IntPoint(sample->info.qp); if (el != prev_el) { @@ -128,6 +146,8 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const if (doftrans) { doftrans->TransformDual(el_y); } AddMultTransposeSubMatrix(*basis, vdofs, el_y, y); + + timer.Stop("Mult/dnfi-slow-eqp"); } // for (int i = 0; i < el_samples->Size(); i++) } // for (int k = 0; k < dnfi.Size(); k++) } // if (dnfi.Size()) @@ -140,26 +160,38 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const for (int k = 0; k < fnfi.Size(); k++) { + timer.Start("Mult/fnfi-init"); + const IntegrationRule *ir = fnfi[k]->GetIntegrationRule(); assert(ir); // we enforce all integrators to set the IntegrationRule a priori. EQPElement *eqp_elem = fnfi_sample[k]; assert(eqp_elem); + timer.Stop("Mult/fnfi-init"); + int prev_face = -1; for (int i = 0; i < eqp_elem->Size(); i++) { + timer.Start("Mult/fnfi-elem-init"); + EQPSample *sample = eqp_elem->GetSample(i); int face = sample->info.el; tr = mesh->GetInteriorFaceTransformations(face); + timer.Stop("Mult/fnfi-elem-init"); + if (precompute) { + timer.Start("Mult/fnfi-fast-eqp"); fnfi[k]->AddAssembleVector_Fast(*sample, *tr, x, y); + timer.Stop("Mult/fnfi-fast-eqp"); continue; } + timer.Start("Mult/fnfi-slow-eqp"); + const IntegrationPoint &ip = ir->IntPoint(sample->info.qp); if (face != prev_face) @@ -187,6 +219,8 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const fnfi[k]->AssembleQuadratureVector(*fe1, *fe2, *tr, ip, sample->info.qw, el_x, el_y); AddMultTransposeSubMatrix(*basis, vdofs, el_y, y); + + timer.Stop("Mult/fnfi-slow-eqp"); } // for (int i = 0; i < sample_info->Size(); i++, sample++) } // for (int k = 0; k < fnfi.Size(); k++) } // if (fnfi.Size()) @@ -219,15 +253,21 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const for (int k = 0; k < bfnfi.Size(); k++) { + timer.Start("Mult/bfnfi-init"); + const IntegrationRule *ir = bfnfi[k]->GetIntegrationRule(); assert(ir); // we enforce all integrators to set the IntegrationRule a priori. EQPElement *eqp_elem = bfnfi_sample[k]; assert(eqp_elem); + timer.Stop("Mult/bfnfi-init"); + int prev_be = -1; for (int i = 0; i < eqp_elem->Size(); i++) { + timer.Start("Mult/bfnfi-elem-init"); + EQPSample *sample = eqp_elem->GetSample(i); int be = sample->info.el; const int bdr_attr = mesh->GetBdrAttribute(be); @@ -243,12 +283,18 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const tr = mesh->GetBdrFaceTransformations (be); + timer.Stop("Mult/bfnfi-elem-init"); + if (precompute) { + timer.Start("Mult/bfnfi-fast-eqp"); bfnfi[k]->AddAssembleVector_Fast(*sample, *tr, x, y); + timer.Stop("Mult/bfnfi-fast-eqp"); continue; } + timer.Start("Mult/bfnfi-slow-eqp"); + const IntegrationPoint &ip = ir->IntPoint(sample->info.qp); if (be != prev_be) @@ -277,6 +323,8 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const bfnfi[k]->AssembleQuadratureVector(*fe1, *fe2, *tr, ip, sample->info.qw, el_x, el_y); AddMultTransposeSubMatrix(*basis, vdofs, el_y, y); + + timer.Stop("Mult/bfnfi-slow-eqp"); } // for (int i = 0; i < sample_info->Size(); i++, sample++) } // for (int k = 0; k < bfnfi.Size(); k++) } // if (bfnfi.Size()) @@ -295,14 +343,14 @@ void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const // // y(ess_tdof_list[i]) = x(ess_tdof_list[i]); // } // // In parallel, the result is in 'py' which is an alias for 'aux2'. + timer.Stop("Mult"); } Operator& ROMNonlinearForm::GetGradient(const Vector &x) const { - timer.Start("Total"); + timer.Start("GetGradient"); - // jac_timers[0]->Start(); - timer.Start("init"); + timer.Start("Grad/init"); assert(x.Size() == width); // if (ext) // { @@ -344,13 +392,13 @@ Operator& ROMNonlinearForm::GetGradient(const Vector &x) const { *Grad = 0.0; } - timer.Stop("init"); + timer.Stop("Grad/init"); if (dnfi.Size()) { for (int k = 0; k < dnfi.Size(); k++) { - timer.Start("sample-load"); + timer.Start("Grad/sample-load"); const IntegrationRule *ir = dnfi[k]->GetIntegrationRule(); assert(ir); // we enforce all integrators to set the IntegrationRule a priori. @@ -358,20 +406,20 @@ Operator& ROMNonlinearForm::GetGradient(const Vector &x) const assert(eqp_elem); int prev_el = -1; - timer.Stop("sample-load"); + timer.Stop("Grad/sample-load"); for (int i = 0; i < eqp_elem->Size(); i++) { - timer.Start("elem-load"); + timer.Start("Grad/elem-load"); EQPSample *sample = eqp_elem->GetSample(i); int el = sample->info.el; T = fes->GetElementTransformation(el); - timer.Stop("elem-load"); + timer.Stop("Grad/elem-load"); - timer.Start("assem-grad"); + timer.Start("Grad/assem-grad"); if (precompute) { dnfi[k]->AddAssembleGrad_Fast(*sample, *T, x, *Grad); - timer.Stop("assem-grad"); + timer.Stop("Grad/assem-grad"); continue; } @@ -392,12 +440,12 @@ Operator& ROMNonlinearForm::GetGradient(const Vector &x) const AddSubMatrixRtAP(*basis, vdofs, elmat, *basis, vdofs, *Grad); // Grad->AddSubMatrix(rom_vdofs, rom_vdofs, quadmat, skip_zeros); - timer.Stop("assem-grad"); + timer.Stop("Grad/assem-grad"); } // for (int i = 0; i < el_samples->Size(); i++) } // for (int k = 0; k < dnfi.Size(); k++) } // if (dnfi.Size()) - timer.Start("other-integ"); + timer.Start("Grad/other-integ"); if (fnfi.Size()) { FaceElementTransformations *tr = NULL; @@ -555,8 +603,8 @@ Operator& ROMNonlinearForm::GetGradient(const Vector &x) const DenseMatrix *mGrad = Grad; - timer.Stop("other-integ"); - timer.Stop("Total"); + timer.Stop("Grad/other-integ"); + timer.Stop("GetGradient"); // TODO(kevin): will need to consider how we should address this case. // do we really need prolongation when we lift up from ROM? // we might need a similar operation on the ROM basis DenseMatrix. diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index fe314ce0..d711c005 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -222,7 +222,7 @@ Operator& SteadyNSTensorROM::GetGradient(const Vector &x) const SteadyNSEQPROM::SteadyNSEQPROM( ROMHandlerBase *rom_handler, Array &hs_, ROMInterfaceForm *itf_, const bool direct_solve_) - : SteadyNSROM(hs_.Size(), rom_handler, direct_solve_), hs(hs_), itf(itf_) + : SteadyNSROM(hs_.Size(), rom_handler, direct_solve_), hs(hs_), itf(itf_), timer() { if (!separate_variable) { @@ -246,6 +246,8 @@ SteadyNSEQPROM::SteadyNSEQPROM( void SteadyNSEQPROM::Mult(const Vector &x, Vector &y) const { + timer.Start("Mult"); + y = 0.0; linearOp->Mult(x, y); @@ -255,19 +257,29 @@ void SteadyNSEQPROM::Mult(const Vector &x, Vector &y) const x_comp.MakeRef(const_cast(x), block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); y_comp.MakeRef(y, block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); + timer.Start("Mult_nlin_domain"); hs[m]->AddMult(x_comp, y_comp); + timer.Stop("Mult_nlin_domain"); } if (!itf) return; GetVel(x, x_u); y_u = 0.0; + + timer.Start("Mult_nlin_itf"); itf->InterfaceAddMult(x_u, y_u); + timer.Stop("Mult_nlin_itf"); + AddVel(y_u, y); + + timer.Stop("Mult"); } Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const { + timer.Start("GetGradient"); + delete jac_mono; delete jac_hypre; jac_mono = new SparseMatrix(*linearOp); @@ -288,7 +300,10 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const } GetVel(x, x_u); + + timer.Start("Grad_nlin_itf"); itf->InterfaceGetGradient(x_u, mats); + timer.Stop("Grad_nlin_itf"); MatrixBlocks u_mats(mats); @@ -307,10 +322,14 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const x_comp.MakeRef(const_cast(x), block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); // NOTE(kevin): jac_comp is owned by hs[m]. No need of deleting it. + timer.Start("Grad_nlin_domain"); jac_comp = dynamic_cast(&hs[m]->GetGradient(x_comp)); + timer.Stop("Grad_nlin_domain"); jac_mono->AddSubMatrix(*block_idxs[m], *block_idxs[m], *jac_comp); } jac_mono->Finalize(); + + timer.Stop("GetGradient"); if (direct_solve) { @@ -424,6 +443,7 @@ SteadyNSSolver::~SteadyNSSolver() else if (rom_handler->GetNonlinearHandling() == NonlinearHandling::EQP) { DeletePointers(comp_eqps); + DeletePointers(subdomain_eqps); delete itf_eqp; } } From c856979c99351da6e9474b677c6cc7485e87a3c7 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Fri, 17 Jan 2025 16:45:45 -0800 Subject: [PATCH 13/20] rom interface form time profiling --- src/interface_form.cpp | 3 +- src/rom_interfaceform.cpp | 66 ++++++++++++++++++++++++++++----------- src/rom_nonlinearform.cpp | 3 +- 3 files changed, 51 insertions(+), 21 deletions(-) diff --git a/src/interface_form.cpp b/src/interface_form.cpp index a2cf7c24..db8254a8 100644 --- a/src/interface_form.cpp +++ b/src/interface_form.cpp @@ -27,7 +27,8 @@ InterfaceForm::~InterfaceForm() { DeletePointers(fnfi); - timer.Print("InterfaceForm::InterfaceAddMult"); + if (config.GetOption("time_profile/enabled", false)) + timer.Print("InterfaceForm::InterfaceAddMult"); } void InterfaceForm::AssembleInterfaceMatrices(Array2D &mats) const diff --git a/src/rom_interfaceform.cpp b/src/rom_interfaceform.cpp index de46149c..e304c210 100644 --- a/src/rom_interfaceform.cpp +++ b/src/rom_interfaceform.cpp @@ -38,7 +38,8 @@ ROMInterfaceForm::~ROMInterfaceForm() { DeletePointers(fnfi_ref_sample); - timer.Print("ROMInterfaceForm::InterfaceAddMult"); + if (config.GetOption("time_profile/enabled", false)) + timer.Print("ROMInterfaceForm"); } void ROMInterfaceForm::SetBasisAtComponent(const int c, DenseMatrix &basis_, const int offset) @@ -75,9 +76,9 @@ void ROMInterfaceForm::UpdateBlockOffsets() void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const { - timer.Start("Total"); + timer.Start("Mult"); - timer.Start("init"); + timer.Start("Mult/init"); assert(block_offsets.Min() >= 0); assert(x.Size() == block_offsets.Last()); @@ -102,22 +103,22 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const Array vdofs1, vdofs2; Vector el_x1, el_x2, el_y1, el_y2; - timer.Stop("init"); + timer.Stop("Mult/init"); for (int k = 0; k < fnfi.Size(); k++) { - timer.Start("itg-init"); + timer.Start("Mult/itg-init"); assert(fnfi[k]); const IntegrationRule *ir = fnfi[k]->GetIntegrationRule(); assert(ir); // we enforce all integrators to set the IntegrationRule a priori. - timer.Stop("itg-init"); + timer.Stop("Mult/itg-init"); for (int p = 0; p < numPorts; p++) { - timer.Start("port-init"); + timer.Start("Mult/port-init"); const PortInfo *pInfo = topol_handler->GetPortInfo(p); @@ -146,14 +147,14 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const eqp_elem = fnfi_sample[p + k * numPorts]; - timer.Stop("port-init"); + timer.Stop("Mult/port-init"); if (!eqp_elem) continue; int prev_itf = -1; for (int i = 0; i < eqp_elem->Size(); i++) { - timer.Start("elem-init"); + timer.Start("Mult/elem-init"); EQPSample *sample = eqp_elem->GetSample(i); @@ -161,24 +162,24 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const if (itf != prev_itf) { InterfaceInfo *if_info = &((*interface_infos)[itf]); - timer.Start("topol"); + timer.Start("Mult/topol"); topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); - timer.Stop("topol"); + timer.Stop("Mult/topol"); prev_itf = itf; } // if (itf != prev_itf) - timer.Stop("elem-init"); + timer.Stop("Mult/elem-init"); if (precompute) { - timer.Start("fast-eqp"); + timer.Start("Mult/fast-eqp"); fnfi[k]->AddAssembleVector_Fast(*sample, *tr1, *tr2, x1, x2, y1, y2); - timer.Stop("fast-eqp"); + timer.Stop("Mult/fast-eqp"); continue; } - timer.Start("slow-eqp"); + timer.Start("Mult/slow-eqp"); if ((tr1 == NULL) || (tr2 == NULL)) mfem_error("InterfaceTransformation of the sampled face is NULL,\n" @@ -209,23 +210,27 @@ void ROMInterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const AddMultTransposeSubMatrix(*basis1, vdofs1, el_y1, y1); AddMultTransposeSubMatrix(*basis2, vdofs2, el_y2, y2); - timer.Stop("slow-eqp"); + timer.Stop("Mult/slow-eqp"); } // for (int i = 0; i < sample_info->Size(); i++, sample++) } // for (int p = 0; p < numPorts; p++) } // for (int k = 0; k < fnfi.Size(); k++) - timer.Start("final"); + timer.Start("Mult/final"); for (int i=0; i < y_tmp.NumBlocks(); ++i) y_tmp.GetBlock(i).SyncAliasMemory(y); - timer.Stop("final"); + timer.Stop("Mult/final"); - timer.Stop("Total"); + timer.Stop("Mult"); } void ROMInterfaceForm::InterfaceGetGradient(const Vector &x, Array2D &mats) const { + timer.Start("GetGradient"); + + timer.Start("Grad/init"); + assert(mats.NumRows() == numSub); assert(mats.NumCols() == numSub); for (int i = 0; i < numSub; i++) @@ -255,15 +260,23 @@ void ROMInterfaceForm::InterfaceGetGradient(const Vector &x, Array2DGetIntegrationRule(); assert(ir); // we enforce all integrators to set the IntegrationRule a priori. + timer.Stop("Grad/fnfi-init"); + for (int p = 0; p < numPorts; p++) { + timer.Start("Grad/port-init"); + const PortInfo *pInfo = topol_handler->GetPortInfo(p); midx[0] = pInfo->Mesh1; @@ -292,24 +305,35 @@ void ROMInterfaceForm::InterfaceGetGradient(const Vector &x, Array2DSize(); i++) { + timer.Start("Grad/elem-init"); + EQPSample *sample = eqp_elem->GetSample(i); int itf = sample->info.el; InterfaceInfo *if_info = &((*interface_infos)[itf]); + timer.Start("Grad/topol"); topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); + timer.Stop("Grad/topol"); const IntegrationPoint &ip = ir->IntPoint(sample->info.qp); + timer.Stop("Grad/elem-init"); + if (precompute) { + timer.Start("Grad/fast-eqp"); fnfi[k]->AddAssembleGrad_Fast(*sample, *tr1, *tr2, x1, x2, mats_p); + timer.Stop("Grad/fast-eqp"); continue; } + timer.Start("Grad/slow-eqp"); + if (itf != prev_itf) { if ((tr1 == NULL) || (tr2 == NULL)) @@ -343,11 +367,15 @@ void ROMInterfaceForm::InterfaceGetGradient(const Vector &x, Array2DSize(); i++, sample++) } // for (int p = 0; p < numPorts; p++) } // for (int k = 0; k < fnfi.Size(); k++) DeletePointers(quadmats); + + timer.Stop("GetGradient"); } void ROMInterfaceForm::TrainEQPForRefPort(const int p, diff --git a/src/rom_nonlinearform.cpp b/src/rom_nonlinearform.cpp index 391e1698..f1fc6e71 100644 --- a/src/rom_nonlinearform.cpp +++ b/src/rom_nonlinearform.cpp @@ -42,7 +42,8 @@ ROMNonlinearForm::~ROMNonlinearForm() DeletePointers(bfnfi_sample); } - timer.Print("ROMNonlinearForm"); + if (config.GetOption("time_profile/enabled", false)) + timer.Print("ROMNonlinearForm"); } void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const From 087f35e5d29816a3a8460246e4e065c3554d4d07 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Wed, 22 Jan 2025 14:21:47 -0800 Subject: [PATCH 14/20] time profiling in DGLaxFriedrichsFluxIntegrator --- include/interfaceinteg.hpp | 5 +++ src/interface_form.cpp | 41 ++++++++++++++++----- src/interfaceinteg.cpp | 74 ++++++++++++++++++++++++++++++++++++++ src/rom_interfaceform.cpp | 4 ++- 4 files changed, 114 insertions(+), 10 deletions(-) diff --git a/include/interfaceinteg.hpp b/include/interfaceinteg.hpp index a48f69c7..6b8372b5 100644 --- a/include/interfaceinteg.hpp +++ b/include/interfaceinteg.hpp @@ -7,6 +7,7 @@ #include "mfem.hpp" #include "hyperreduction_integ.hpp" +#include "etc.hpp" namespace mfem { @@ -298,10 +299,14 @@ class DGLaxFriedrichsFluxIntegrator : public InterfaceNonlinearFormIntegrator DenseMatrix udof1, udof2, elv1, elv2; DenseMatrix elmat_comp11, elmat_comp12, elmat_comp21, elmat_comp22; + mutable TimeProfiler timer; + public: DGLaxFriedrichsFluxIntegrator(Coefficient &q, VectorCoefficient *ud = NULL, const IntegrationRule *ir = NULL) : InterfaceNonlinearFormIntegrator(ir), Q(&q), UD(ud) {} + virtual ~DGLaxFriedrichsFluxIntegrator(); + void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, diff --git a/src/interface_form.cpp b/src/interface_form.cpp index db8254a8..4bf1cc12 100644 --- a/src/interface_form.cpp +++ b/src/interface_form.cpp @@ -27,7 +27,7 @@ InterfaceForm::~InterfaceForm() { DeletePointers(fnfi); - if (config.GetOption("time_profile/enabled", false)) + if (config.GetOption("time_profile/InterfaceForm", false)) timer.Print("InterfaceForm::InterfaceAddMult"); } @@ -68,8 +68,8 @@ void InterfaceForm::AssembleInterfaceMatrices(Array2D &mats) con void InterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const { - timer.Start("Total"); - timer.Start("init"); + timer.Start("Mult"); + timer.Start("Mult/init"); x_tmp.Update(const_cast(x), block_offsets); y_tmp.Update(y, block_offsets); @@ -78,11 +78,11 @@ void InterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const Mesh *mesh1, *mesh2; FiniteElementSpace *fes1, *fes2; - timer.Stop("init"); + timer.Stop("Mult/init"); for (int p = 0; p < topol_handler->GetNumPorts(); p++) { - timer.Start("port-init"); + timer.Start("Mult/port-init"); const PortInfo *pInfo = topol_handler->GetPortInfo(p); @@ -97,24 +97,28 @@ void InterfaceForm::InterfaceAddMult(const Vector &x, Vector &y) const Array* const interface_infos = topol_handler->GetInterfaceInfos(p); - timer.Stop("port-init"); + timer.Stop("Mult/port-init"); AssembleInterfaceVector(mesh1, mesh2, fes1, fes2, interface_infos, x_tmp.GetBlock(midx[0]), x_tmp.GetBlock(midx[1]), y_tmp.GetBlock(midx[0]), y_tmp.GetBlock(midx[1])); } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) - timer.Start("final"); + timer.Start("Mult/final"); for (int i=0; i < y_tmp.NumBlocks(); ++i) y_tmp.GetBlock(i).SyncAliasMemory(y); - timer.Stop("final"); - timer.Stop("Total"); + timer.Stop("Mult/final"); + timer.Stop("Mult"); } void InterfaceForm::InterfaceGetGradient(const Vector &x, Array2D &mats) const { + timer.Start("GetGradient"); + + timer.Start("Grad/init"); + assert(mats.NumRows() == numSub); assert(mats.NumCols() == numSub); for (int i = 0; i < numSub; i++) @@ -128,8 +132,12 @@ void InterfaceForm::InterfaceGetGradient(const Vector &x, Array2DGetNumPorts(); p++) { + timer.Start("Grad/port-init"); + const PortInfo *pInfo = topol_handler->GetPortInfo(p); midx[0] = pInfo->Mesh1; @@ -146,9 +154,12 @@ void InterfaceForm::InterfaceGetGradient(const Vector &x, Array2D* const interface_infos = topol_handler->GetInterfaceInfos(p); + timer.Stop("Grad/port-init"); AssembleInterfaceGrad(mesh1, mesh2, fes1, fes2, interface_infos, x_tmp.GetBlock(midx[0]), x_tmp.GetBlock(midx[1]), mats_p); } // for (int p = 0; p < topol_handler->GetNumPorts(); p++) + + timer.Stop("GetGradient"); } void InterfaceForm::AssembleInterfaceMatrix( @@ -298,6 +309,8 @@ void InterfaceForm::AssembleInterfaceGrad(Mesh *mesh1, Mesh *mesh2, FiniteElementSpace *fes1, FiniteElementSpace *fes2, Array *interface_infos, const Vector &x1, const Vector &x2, Array2D &mats) const { + timer.Start("Grad/assemble-init"); + assert(x1.Size() == fes1->GetTrueVSize()); assert(x2.Size() == fes2->GetTrueVSize()); @@ -309,14 +322,22 @@ void InterfaceForm::AssembleInterfaceGrad(Mesh *mesh1, Mesh *mesh2, Array vdofs1, vdofs2; Vector el_x1, el_x2, el_y1, el_y2; + timer.Stop("Grad/assemble-init"); + for (int bn = 0; bn < interface_infos->Size(); bn++) { + timer.Start("Grad/topol"); + InterfaceInfo *if_info = &((*interface_infos)[bn]); topol_handler->GetInterfaceTransformations(mesh1, mesh2, if_info, tr1, tr2); + timer.Stop("Grad/topol"); + if ((tr1 != NULL) && (tr2 != NULL)) { + timer.Start("Grad/assemble"); + fes1->GetElementVDofs(tr1->Elem1No, vdofs1); fes2->GetElementVDofs(tr2->Elem1No, vdofs2); @@ -338,6 +359,8 @@ void InterfaceForm::AssembleInterfaceGrad(Mesh *mesh1, Mesh *mesh2, mats(1, 0)->AddSubMatrix(vdofs2, vdofs1, *elemmats(1, 0), skip_zeros); mats(1, 1)->AddSubMatrix(vdofs2, vdofs2, *elemmats(1, 1), skip_zeros); } + + timer.Stop("Grad/assemble"); } // if ((tr1 != NULL) && (tr2 != NULL)) } // for (int bn = 0; bn < interface_infos.Size(); bn++) diff --git a/src/interfaceinteg.cpp b/src/interfaceinteg.cpp index ae9d7033..50059b79 100644 --- a/src/interfaceinteg.cpp +++ b/src/interfaceinteg.cpp @@ -5,6 +5,7 @@ #include "interfaceinteg.hpp" #include "etc.hpp" #include "linalg_utils.hpp" +#include "input_parser.hpp" // #include // #include @@ -1307,6 +1308,12 @@ void InterfaceDGElasticityIntegrator::AssembleBlock( DGLaxFriedrichsFluxIntegrator */ +DGLaxFriedrichsFluxIntegrator::~DGLaxFriedrichsFluxIntegrator() +{ + if (config.GetOption("time_profile/DGLaxFriedrichsFluxIntegrator", false)) + timer.Print("DGLaxFriedrichsFluxIntegrator"); +} + void DGLaxFriedrichsFluxIntegrator::ComputeFluxDotN( const Vector &u1, const Vector &u2, const Vector &nor, const bool &eval2, Vector &flux) @@ -1527,6 +1534,8 @@ void DGLaxFriedrichsFluxIntegrator::AssembleFaceVector( const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, Vector &elvect) { + timer.Start("AssembleFaceVector_interior"); + dim = el1.GetDim(); ndofs1 = el1.GetDof(); ndofs2 = (Tr.Elem2No >= 0) ? el2.GetDof() : 0; @@ -1564,12 +1573,16 @@ void DGLaxFriedrichsFluxIntegrator::AssembleFaceVector( AssembleQuadVectorBase(el1, el2, &Tr, NULL, ip, ip.weight, ndofs2, udof1, udof2, elv1, elv2); } + + timer.Stop("AssembleFaceVector_interior"); } void DGLaxFriedrichsFluxIntegrator::AssembleFaceGrad( const FiniteElement &el1, const FiniteElement &el2, FaceElementTransformations &Tr, const Vector &elfun, DenseMatrix &elmat) { + timer.Start("AssembleFaceGrad_interior"); + dim = el1.GetDim(); ndofs1 = el1.GetDof(); ndofs2 = (Tr.Elem2No >= 0) ? el2.GetDof() : 0; @@ -1627,6 +1640,8 @@ void DGLaxFriedrichsFluxIntegrator::AssembleFaceGrad( } } } // for (int pind = 0; pind < ir->GetNPoints(); ++pind) + + timer.Stop("AssembleFaceGrad_interior"); } void DGLaxFriedrichsFluxIntegrator::AssembleQuadratureVector( @@ -1712,6 +1727,8 @@ void DGLaxFriedrichsFluxIntegrator::AssembleQuadratureGrad( void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( const EQPSample &eqp_sample, FaceElementTransformations &T, const Vector &x, Vector &y) { + timer.Start("eqpvec_Fast_interior"); + const IntegrationPoint &ip = GetIntegrationRule()->IntPoint(eqp_sample.info.qp); const double qw = eqp_sample.info.qw; DenseMatrix *shapes1 = eqp_sample.shape1; @@ -1756,11 +1773,15 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( assert(y.Size() == x.Size()); shapes1->AddMultTranspose(flux, y, -w); if (el2) shapes2->AddMultTranspose(flux, y, w); + + timer.Stop("eqpvec_Fast_interior"); } void DGLaxFriedrichsFluxIntegrator::AddAssembleGrad_Fast( const EQPSample &eqp_sample, FaceElementTransformations &T, const Vector &x, DenseMatrix &jac) { + timer.Start("eqpgrad_Fast_interior"); + const IntegrationPoint &ip = GetIntegrationRule()->IntPoint(eqp_sample.info.qp); const double qw = eqp_sample.info.qw; DenseMatrix *shapes1 = eqp_sample.shape1; @@ -1812,6 +1833,8 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleGrad_Fast( AddwRtAP(*shapes1, gradu2, *shapes2, jac, -w); AddwRtAP(*shapes2, gradu2, *shapes2, jac, w); } + + timer.Stop("eqpgrad_Fast_interior"); } void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceVector( @@ -1820,6 +1843,8 @@ void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceVector( const Vector &elfun1, const Vector &elfun2, Vector &elvect1, Vector &elvect2) { + timer.Start("fomvec_interface"); + dim = el1.GetDim(); ndofs1 = el1.GetDof(); ndofs2 = el2.GetDof(); @@ -1855,6 +1880,8 @@ void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceVector( AssembleQuadVectorBase(el1, el2, &Tr1, &Tr2, ip, ip.weight, ndofs2, udof1, udof2, elv1, elv2); } + + timer.Stop("fomvec_interface"); } void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceGrad( @@ -1862,6 +1889,10 @@ void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceGrad( FaceElementTransformations &Tr1, FaceElementTransformations &Tr2, const Vector &elfun1, const Vector &elfun2, Array2D &elmats) { + timer.Start("fomgrad_interface"); + + timer.Start("fomgrad_interface/init"); + assert(elmats.NumRows() == 2); assert(elmats.NumCols() == 2); for (int i = 0; i < 2; i++) @@ -1905,13 +1936,20 @@ void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceGrad( for (int i = 0; i < 2; i++) for (int j = 0; j < 2; j++) *(elmats(i, j)) = 0.0; + timer.Stop("fomgrad_interface/init"); + for (int pind = 0; pind < ir->GetNPoints(); ++pind) { + timer.Start("fomgrad_interface/eval"); + const IntegrationPoint &ip = ir->IntPoint(pind); AssembleQuadGradBase(el1, el2, &Tr1, &Tr2, ip, ip.weight, ndofs2, udof1, udof2, w, gradu1, gradu2, elmat_comp11, elmat_comp12, elmat_comp21, elmat_comp22); + timer.Stop("fomgrad_interface/eval"); + timer.Start("fomgrad_interface/assemble"); + for (int di = 0; di < dim; di++) { for (int dj = 0; dj < dim; dj++) @@ -1925,7 +1963,11 @@ void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceGrad( } } } + + timer.Stop("fomgrad_interface/assemble"); } // for (int pind = 0; pind < ir->GetNPoints(); ++pind) + + timer.Stop("fomgrad_interface"); } void DGLaxFriedrichsFluxIntegrator::AssembleQuadratureVector( @@ -2023,6 +2065,8 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( FaceElementTransformations &Tr2, const Vector &x1, const Vector &x2, Vector &y1, Vector &y2) { + timer.Start("eqpvec_Fast_interface"); + const IntegrationRule *ir = GetIntegrationRule(); const IntegrationPoint &ip = ir->IntPoint(eqp_sample.info.qp); const double qw = eqp_sample.info.qw; @@ -2062,6 +2106,8 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( assert(y2.Size() == x2.Size()); shapes1->AddMultTranspose(flux, y1, -w); shapes2->AddMultTranspose(flux, y2, w); + + timer.Stop("eqpvec_Fast_interface"); } void DGLaxFriedrichsFluxIntegrator::AddAssembleGrad_Fast( @@ -2069,6 +2115,9 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleGrad_Fast( FaceElementTransformations &Tr2, const Vector &x1, const Vector &x2, Array2D &jac) { + timer.Start("eqpgrad_Fast_interface"); + timer.Start("eqpgrad_Fast_interface/init"); + const IntegrationRule *ir = GetIntegrationRule(); const IntegrationPoint &ip = ir->IntPoint(eqp_sample.info.qp); const double qw = eqp_sample.info.qw; @@ -2089,6 +2138,9 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleGrad_Fast( const IntegrationPoint &eip1 = Tr1.GetElement1IntPoint(); // const IntegrationPoint &eip2 = T.GetElement2IntPoint(); + timer.Stop("eqpgrad_Fast_interface/init"); + timer.Start("eqpgrad_Fast_interface/eval"); + shapes1->Mult(x1, u1); shapes2->Mult(x2, u2); @@ -2108,10 +2160,32 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleGrad_Fast( if (Q) w *= Q->Eval(Tr1, ip); + timer.Stop("eqpgrad_Fast_interface/eval"); +{ + Array2D test_jac(2, 2); + for (int i = 0; i < 2; i++) + for (int j = 0; j < 2; j++) + test_jac(i, j) = new DenseMatrix(jac(i, j)->Height(), jac(i, j)->Width()); + + timer.Start("eqpgrad_Fast_interface/assemble_test"); + AddwRtAP(*shapes1, gradu1, *shapes1, *test_jac(0, 0), -w); + AddwRtAP(*shapes2, gradu1, *shapes1, *test_jac(1, 0), w); + AddwRtAP(*shapes1, gradu2, *shapes2, *test_jac(0, 1), -w); + AddwRtAP(*shapes2, gradu2, *shapes2, *test_jac(1, 1), w); + timer.Stop("eqpgrad_Fast_interface/assemble_test"); + + DeletePointers(test_jac); +} + timer.Start("eqpgrad_Fast_interface/assemble"); + AddwRtAP(*shapes1, gradu1, *shapes1, *jac(0, 0), -w); AddwRtAP(*shapes2, gradu1, *shapes1, *jac(1, 0), w); AddwRtAP(*shapes1, gradu2, *shapes2, *jac(0, 1), -w); AddwRtAP(*shapes2, gradu2, *shapes2, *jac(1, 1), w); + + timer.Stop("eqpgrad_Fast_interface/assemble"); + + timer.Stop("eqpgrad_Fast_interface"); } } diff --git a/src/rom_interfaceform.cpp b/src/rom_interfaceform.cpp index e304c210..c85625c7 100644 --- a/src/rom_interfaceform.cpp +++ b/src/rom_interfaceform.cpp @@ -38,7 +38,7 @@ ROMInterfaceForm::~ROMInterfaceForm() { DeletePointers(fnfi_ref_sample); - if (config.GetOption("time_profile/enabled", false)) + if (config.GetOption("time_profile/ROMInterfaceForm", false)) timer.Print("ROMInterfaceForm"); } @@ -308,6 +308,8 @@ void ROMInterfaceForm::InterfaceGetGradient(const Vector &x, Array2DSize(), eqp_elem->Size()); + int prev_itf = -1; for (int i = 0; i < eqp_elem->Size(); i++) { From 46cdccbafd636db33d13ab66f2eb03c5453992b4 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Fri, 7 Feb 2025 16:41:27 -0800 Subject: [PATCH 15/20] PrintEQPCoords mode --- bin/main.cpp | 1 + include/main_workflow.hpp | 3 +++ include/multiblock_solver.hpp | 2 ++ include/rom_nonlinearform.hpp | 2 ++ include/steady_ns_solver.hpp | 2 ++ src/main_workflow.cpp | 49 +++++++++++++++++++++++++++++++++++ src/rom_interfaceform.cpp | 2 +- src/rom_nonlinearform.cpp | 45 +++++++++++++++++++++++++++++++- src/steady_ns_solver.cpp | 49 +++++++++++++++++++++++++++++++++++ 9 files changed, 153 insertions(+), 2 deletions(-) diff --git a/bin/main.cpp b/bin/main.cpp index 2977c26f..ced806aa 100644 --- a/bin/main.cpp +++ b/bin/main.cpp @@ -34,6 +34,7 @@ int main(int argc, char *argv[]) else if (mode == "aux_train_rom") AuxiliaryTrainROM(MPI_COMM_WORLD); else if (mode == "train_eqp") TrainEQP(MPI_COMM_WORLD); else if (mode == "single_run") double dump = SingleRun(MPI_COMM_WORLD, output_file); + else if (mode == "print_eqp") PrintEQPCoords(MPI_COMM_WORLD); else { if (rank == 0) printf("Unknown mode %s!\n", mode.c_str()); diff --git a/include/main_workflow.hpp b/include/main_workflow.hpp index e903b75d..c08684b7 100644 --- a/include/main_workflow.hpp +++ b/include/main_workflow.hpp @@ -32,4 +32,7 @@ void FindSnapshotFilesForBasis(const BasisTag &basis_tag, const std::string &def // return relative error if comparing solution. double SingleRun(MPI_Comm comm, const std::string output_file = ""); +// Auxiliary function to print out EQP point coordinates +void PrintEQPCoords(MPI_Comm comm); + #endif diff --git a/include/multiblock_solver.hpp b/include/multiblock_solver.hpp index c215ef49..ea327e0d 100644 --- a/include/multiblock_solver.hpp +++ b/include/multiblock_solver.hpp @@ -260,6 +260,8 @@ friend class ParameterizedProblem; void ComputeRelativeError(Array fom_sols, Array rom_sols, Vector &error); void CompareSolution(BlockVector &test_U, Vector &error); + virtual void SaveEQPCoords(const std::string &filename) {} + protected: virtual void AssembleROMMat(BlockMatrix &romMat); }; diff --git a/include/rom_nonlinearform.hpp b/include/rom_nonlinearform.hpp index 20db9d03..2bbe8850 100644 --- a/include/rom_nonlinearform.hpp +++ b/include/rom_nonlinearform.hpp @@ -233,6 +233,8 @@ class ROMNonlinearForm : public NonlinearForm The state @a x must be a true-dof vector. */ virtual Operator &GetGradient(const Vector &x) const; + void SaveDomainEQPCoords(const int k, hid_t file_id, const std::string &dsetname); + private: void PrecomputeDomainEQPSample(const IntegrationRule &ir, const DenseMatrix &basis, EQPSample &eqp_sample); void PrecomputeFaceEQPSample(const IntegrationRule &ir, const DenseMatrix &basis, diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 2531ebe3..68962c38 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -196,6 +196,8 @@ friend class SteadyNSOperator; virtual void LoadROMNlinElems(const std::string &input_prefix) override; virtual void AssembleROMNlinOper() override; + void SaveEQPCoords(const std::string &filename) override; + private: DenseTensor* GetReducedTensor(DenseMatrix *basis, FiniteElementSpace *fespace); diff --git a/src/main_workflow.cpp b/src/main_workflow.cpp index 57032254..cec992a6 100644 --- a/src/main_workflow.cpp +++ b/src/main_workflow.cpp @@ -702,3 +702,52 @@ double SingleRun(MPI_Comm comm, const std::string output_file) // return the maximum error over all variables. return error.Max(); } + +void PrintEQPCoords(MPI_Comm comm) +{ + ParameterizedProblem *problem = InitParameterizedProblem(); + MultiBlockSolver *test = InitSolver(); + test->InitVariables(); + + assert(test->UseRom()); + if (!test->IsNonlinear()) + { + printf("PrintEQPPoints: given solver is not nonlinear. Exiting\n"); + return; + } + + test->InitROMHandler(); + + StopWatch solveTimer; + + problem->SetSingleRun(); + test->SetParameterizedProblem(problem); + + // TODO: there are skippable operations depending on rom/fom mode. + test->BuildRHSOperators(); + test->SetupRHSBCOperators(); + test->AssembleRHS(); + + const int num_var = test->GetNumVar(); + + ROMHandlerBase *rom = test->GetROMHandler(); + test->LoadReducedBasis(); + test->AllocateROMNlinElems(); + + ROMBuildingLevel save_operator = rom->GetBuildingLevel(); + TopologyHandlerMode topol_mode = test->GetTopologyMode(); + assert(save_operator == ROMBuildingLevel::COMPONENT); + assert(topol_mode != TopologyHandlerMode::SUBMESH); + assert(rom->GetNonlinearHandling() == NonlinearHandling::EQP); + + test->LoadROMNlinElems(rom->GetOperatorPrefix()); + + const std::string filename = config.GetRequiredOption("model_reduction/eqp/save_coords"); + test->SaveEQPCoords(filename); + + delete test; + delete problem; + + // return the maximum error over all variables. + return; +} diff --git a/src/rom_interfaceform.cpp b/src/rom_interfaceform.cpp index c85625c7..06359833 100644 --- a/src/rom_interfaceform.cpp +++ b/src/rom_interfaceform.cpp @@ -696,7 +696,7 @@ void ROMInterfaceForm::TrainEQPForIntegrator( samples.Append({.el = i / nqe, .qp = i % nqe, .qw = eqpSol_global(i)}); } } - printf("Size of sampled qp: %d\n", samples.Size()); + printf("Size of sampled qp: %d/%d\n", samples.Size(), eqpSol_global.dim()); if (nnz != samples.Size()) printf("Sample quadrature points with weight < 1.0e-12 are neglected.\n"); } diff --git a/src/rom_nonlinearform.cpp b/src/rom_nonlinearform.cpp index f1fc6e71..bf8a09c0 100644 --- a/src/rom_nonlinearform.cpp +++ b/src/rom_nonlinearform.cpp @@ -891,7 +891,7 @@ void ROMNonlinearForm::TrainEQPForIntegrator( samples.Append({.el = i / nqe, .qp = i % nqe, .qw = eqpSol_global(i)}); } } - printf("Size of sampled qp: %d\n", samples.Size()); + printf("Size of sampled qp: %d/%d\n", samples.Size(), eqpSol_global.dim()); if (nnz != samples.Size()) printf("Sample quadrature points with weight < 1.0e-12 are neglected.\n"); } @@ -1521,4 +1521,47 @@ void ROMNonlinearForm::PrecomputeBdrFaceEQPSample( PrecomputeFaceEQPSample(ir, basis, T, eqp_sample); } +void ROMNonlinearForm::SaveDomainEQPCoords(const int k, hid_t file_id, const std::string &dsetname) +{ + if (!dnfi.Size()) + return; + + assert((k >= 0) && (k < dnfi.Size())); + + const IntegrationRule *ir = dnfi[k]->GetIntegrationRule(); + EQPElement *eqp_elem = dnfi_sample[k]; + assert(ir); + assert(eqp_elem); + + DenseMatrix eqp_coords(eqp_elem->Size(), 3); + + EQPSample *eqp_sample; + for (int i = 0; i < eqp_elem->Size(); i++) + { + eqp_sample = eqp_elem->GetSample(i); + const int el = eqp_sample->info.el; + const FiniteElement *fe = fes->GetFE(el); + ElementTransformation *T = fes->GetElementTransformation(el); + const IntegrationPoint &ip = ir->IntPoint(eqp_sample->info.qp); + + double x[3]; + Vector transip(x, 3); + T->Transform(ip, transip); + + eqp_coords.SetRow(i, transip); + } // for (int i = 0; i < eqp_elem->Size(); i++) + + assert(file_id >= 0); + hid_t grp_id; + herr_t errf; + + grp_id = H5Gcreate(file_id, dsetname.c_str(), H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(grp_id >= 0); + + hdf5_utils::WriteDataset(grp_id, "coords", eqp_coords); + + errf = H5Gclose(grp_id); + assert(errf >= 0); +} + } diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index d711c005..fa66c877 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -1388,3 +1388,52 @@ DenseTensor* SteadyNSSolver::GetReducedTensor(DenseMatrix *basis, FiniteElementS return tensor; } + +void SteadyNSSolver::SaveEQPCoords(const std::string &filename) +{ + assert(topol_mode == TopologyHandlerMode::COMPONENT); + + /* + TODO(kevin): this is a boilerplate for parallel POD/EQP training. + Full parallelization will save EQ points/weights in a parallel way. + */ + if (rank == 0) + { + hid_t file_id; + herr_t errf = 0; + file_id = H5Fcreate(filename.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT); + assert(file_id >= 0); + + hid_t grp_id; + + grp_id = H5Gcreate(file_id, "components", H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT); + assert(grp_id >= 0); + + const int num_comp = topol_handler->GetNumComponents(); + assert(comp_eqps.Size() == num_comp); + + hdf5_utils::WriteAttribute(grp_id, "number_of_components", num_comp); + + std::string dset_name; + for (int c = 0; c < num_comp; c++) + { + assert(comp_eqps[c]); + dset_name = topol_handler->GetComponentName(c); + + comp_eqps[c]->SaveDomainEQPCoords(0, grp_id, dset_name + "_integ0"); + if (oper_type == OperType::LF) + { + assert(!full_dg); + // comp_eqps[c]->SaveEQPForIntegrator(IntegratorType::INTERIORFACE, 0, grp_id, dset_name + "_integ1"); + } // if (oper_type == OperType::LF) + } // for (int c = 0; c < num_comp; c++) + + errf = H5Gclose(grp_id); + assert(errf >= 0); + + errf = H5Fclose(file_id); + assert(errf >= 0); + } + MPI_Barrier(MPI_COMM_WORLD); + return; +} From 2c7d904ff7fa2158ca09ed0e4556d054a77f38a9 Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Tue, 18 Feb 2025 11:25:18 -0800 Subject: [PATCH 16/20] interface time profiling. --- src/interfaceinteg.cpp | 54 ++++++++++++++++++++++++++++++------------ 1 file changed, 39 insertions(+), 15 deletions(-) diff --git a/src/interfaceinteg.cpp b/src/interfaceinteg.cpp index 50059b79..dacd34d3 100644 --- a/src/interfaceinteg.cpp +++ b/src/interfaceinteg.cpp @@ -1404,6 +1404,8 @@ void DGLaxFriedrichsFluxIntegrator::AssembleQuadVectorBase( const DenseMatrix &elfun1, const DenseMatrix &elfun2, DenseMatrix &elvect1, DenseMatrix &elvect2) { + timer.Start("fomvec_interface/eval"); + assert(Tr1); assert(elfun1.NumRows() == shape1.Size()); assert(elfun1.NumCols() == u1.Size()); @@ -1455,12 +1457,17 @@ void DGLaxFriedrichsFluxIntegrator::AssembleQuadVectorBase( ComputeFluxDotN(u1, u2, nor, eval2, flux); + timer.Stop("fomvec_interface/eval"); + timer.Start("fomvec_interface/assemble"); + w = iw; if (Q) { w *= Q->Eval(*Tr1, ip); } AddMult_a_VWt(-w, shape1, flux, elvect1); if (ndofs2) AddMult_a_VWt(w, shape2, flux, elvect2); + + timer.Stop("fomvec_interface/assemble"); } void DGLaxFriedrichsFluxIntegrator::AssembleQuadGradBase( @@ -1728,6 +1735,7 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( const EQPSample &eqp_sample, FaceElementTransformations &T, const Vector &x, Vector &y) { timer.Start("eqpvec_Fast_interior"); + timer.Start("eqpvec_Fast_interior/eval"); const IntegrationPoint &ip = GetIntegrationRule()->IntPoint(eqp_sample.info.qp); const double qw = eqp_sample.info.qw; @@ -1767,6 +1775,9 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( ComputeFluxDotN(u1, u2, nor, eval2, flux); + timer.Stop("eqpvec_Fast_interior/eval"); + timer.Start("eqpvec_Fast_interior/assemble"); + w = qw; if (Q) { w *= Q->Eval(T, ip); } @@ -1774,6 +1785,7 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( shapes1->AddMultTranspose(flux, y, -w); if (el2) shapes2->AddMultTranspose(flux, y, w); + timer.Stop("eqpvec_Fast_interior/assemble"); timer.Stop("eqpvec_Fast_interior"); } @@ -1845,6 +1857,8 @@ void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceVector( { timer.Start("fomvec_interface"); + timer.Start("fomvec_interface/init"); + dim = el1.GetDim(); ndofs1 = el1.GetDof(); ndofs2 = el2.GetDof(); @@ -1873,6 +1887,8 @@ void DGLaxFriedrichsFluxIntegrator::AssembleInterfaceVector( ir = &IntRules.Get(Tr1.GetGeometryType(), order); } + timer.Stop("fomvec_interface/init"); + elvect1 = 0.0; elvect2 = 0.0; for (int pind = 0; pind < ir->GetNPoints(); ++pind) { @@ -2066,6 +2082,7 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( const Vector &x1, const Vector &x2, Vector &y1, Vector &y2) { timer.Start("eqpvec_Fast_interface"); + timer.Start("eqpvec_Fast_interior/init"); const IntegrationRule *ir = GetIntegrationRule(); const IntegrationPoint &ip = ir->IntPoint(eqp_sample.info.qp); @@ -2085,6 +2102,9 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( const IntegrationPoint &eip1 = Tr1.GetElement1IntPoint(); // const IntegrationPoint &eip2 = T.GetElement2IntPoint(); + timer.Stop("eqpvec_Fast_interior/init"); + timer.Start("eqpvec_Fast_interior/eval"); + shapes1->Mult(x1, u1); shapes2->Mult(x2, u2); @@ -2099,6 +2119,9 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( ComputeFluxDotN(u1, u2, nor, true, flux); + timer.Stop("eqpvec_Fast_interior/eval"); + timer.Start("eqpvec_Fast_interior/assemble"); + w = qw; if (Q) { w *= Q->Eval(Tr1, ip); } @@ -2107,6 +2130,7 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleVector_Fast( shapes1->AddMultTranspose(flux, y1, -w); shapes2->AddMultTranspose(flux, y2, w); + timer.Stop("eqpvec_Fast_interior/assemble"); timer.Stop("eqpvec_Fast_interface"); } @@ -2161,21 +2185,21 @@ void DGLaxFriedrichsFluxIntegrator::AddAssembleGrad_Fast( w *= Q->Eval(Tr1, ip); timer.Stop("eqpgrad_Fast_interface/eval"); -{ - Array2D test_jac(2, 2); - for (int i = 0; i < 2; i++) - for (int j = 0; j < 2; j++) - test_jac(i, j) = new DenseMatrix(jac(i, j)->Height(), jac(i, j)->Width()); - - timer.Start("eqpgrad_Fast_interface/assemble_test"); - AddwRtAP(*shapes1, gradu1, *shapes1, *test_jac(0, 0), -w); - AddwRtAP(*shapes2, gradu1, *shapes1, *test_jac(1, 0), w); - AddwRtAP(*shapes1, gradu2, *shapes2, *test_jac(0, 1), -w); - AddwRtAP(*shapes2, gradu2, *shapes2, *test_jac(1, 1), w); - timer.Stop("eqpgrad_Fast_interface/assemble_test"); - - DeletePointers(test_jac); -} +// { +// Array2D test_jac(2, 2); +// for (int i = 0; i < 2; i++) +// for (int j = 0; j < 2; j++) +// test_jac(i, j) = new DenseMatrix(jac(i, j)->Height(), jac(i, j)->Width()); + +// timer.Start("eqpgrad_Fast_interface/assemble_test"); +// AddwRtAP(*shapes1, gradu1, *shapes1, *test_jac(0, 0), -w); +// AddwRtAP(*shapes2, gradu1, *shapes1, *test_jac(1, 0), w); +// AddwRtAP(*shapes1, gradu2, *shapes2, *test_jac(0, 1), -w); +// AddwRtAP(*shapes2, gradu2, *shapes2, *test_jac(1, 1), w); +// timer.Stop("eqpgrad_Fast_interface/assemble_test"); + +// DeletePointers(test_jac); +// } timer.Start("eqpgrad_Fast_interface/assemble"); AddwRtAP(*shapes1, gradu1, *shapes1, *jac(0, 0), -w); From 87e95ea15a2ca42cebd78785a4e0e24957cb41ed Mon Sep 17 00:00:00 2001 From: Seung Whan Chung Date: Tue, 18 Feb 2025 12:53:49 -0800 Subject: [PATCH 17/20] steady NS time profiling --- include/steady_ns_solver.hpp | 6 +++- src/rom_nonlinearform.cpp | 2 +- src/steady_ns_solver.cpp | 56 +++++++++++++++++++++++++++++------- 3 files changed, 52 insertions(+), 12 deletions(-) diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 68962c38..1fe70c98 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -18,6 +18,9 @@ using namespace mfem; // Ultimately, we should implement InterfaceForm to pass, not the MultiBlockSolver itself. class SteadyNSOperator : public Operator { +private: + mutable TimeProfiler timer; + protected: bool direct_solve; @@ -117,7 +120,8 @@ class SteadyNSEQPROM : public SteadyNSROM virtual ~SteadyNSEQPROM() { - // timer.Print("SteadyNSEQPROM"); + if (config.GetOption("time_profile/SteadyNSEQPROM", false)) + timer.Print("SteadyNSEQPROM"); } virtual void Mult(const Vector &x, Vector &y) const; diff --git a/src/rom_nonlinearform.cpp b/src/rom_nonlinearform.cpp index bf8a09c0..75342753 100644 --- a/src/rom_nonlinearform.cpp +++ b/src/rom_nonlinearform.cpp @@ -42,7 +42,7 @@ ROMNonlinearForm::~ROMNonlinearForm() DeletePointers(bfnfi_sample); } - if (config.GetOption("time_profile/enabled", false)) + if (config.GetOption("time_profile/ROMNonlinearForm", false)) timer.Print("ROMNonlinearForm"); } diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index fa66c877..5680f577 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -57,23 +57,41 @@ SteadyNSOperator::~SteadyNSOperator() delete uu_mono; delete mono_jac; delete jac_hypre; + + if (config.GetOption("time_profile/SteadyNSOperator", false)) + timer.Print("SteadyNSOperator"); } void SteadyNSOperator::Mult(const Vector &x, Vector &y) const { + timer.Start("Mult"); + assert(linearOp); x_u.MakeRef(const_cast(x), 0, u_offsets.Last()); y_u.MakeRef(y, 0, u_offsets.Last()); y = 0.0; + timer.Start("Mult/nlin_domain"); Hop->Mult(x_u, y_u); - if (nl_itf) nl_itf->InterfaceAddMult(x_u, y_u); + timer.Stop("Mult/nlin_domain"); + if (nl_itf) + { + timer.Start("Mult/itf"); + nl_itf->InterfaceAddMult(x_u, y_u); + timer.Stop("Mult/itf"); + } + timer.Start("Mult/lin"); linearOp->AddMult(x, y); + timer.Stop("Mult/lin"); + + timer.Stop("Mult"); } Operator& SteadyNSOperator::GetGradient(const Vector &x) const { + timer.Start("GetGradient"); + // NonlinearForm owns the gradient operator. // DeletePointers(hs_mats); delete hs_jac; @@ -82,6 +100,7 @@ Operator& SteadyNSOperator::GetGradient(const Vector &x) const delete mono_jac; delete jac_hypre; + timer.Start("GetGradient/nlin-jac"); hs_jac = new BlockMatrix(u_offsets); for (int i = 0; i < hs.Size(); i++) for (int j = 0; j < hs.Size(); j++) @@ -97,9 +116,12 @@ Operator& SteadyNSOperator::GetGradient(const Vector &x) const hs_mats(i, j) = new SparseMatrix(u_offsets[i+1] - u_offsets[i], u_offsets[j+1] - u_offsets[j]); } } + timer.Stop("GetGradient/nlin-jac"); + timer.Start("GetGradient/itf-jac"); x_u.MakeRef(const_cast(x), 0, u_offsets.Last()); if (nl_itf) nl_itf->InterfaceGetGradient(x_u, hs_mats); + timer.Stop("GetGradient/itf-jac"); for (int i = 0; i < hs.Size(); i++) for (int j = 0; j < hs.Size(); j++) @@ -108,9 +130,11 @@ Operator& SteadyNSOperator::GetGradient(const Vector &x) const hs_jac->SetBlock(i, j, hs_mats(i, j)); } + timer.Start("GetGradient/lin-jac"); SparseMatrix *hs_jac_mono = hs_jac->CreateMonolithic(); uu_mono = Add(*M, *hs_jac_mono); delete hs_jac_mono; + timer.Stop("GetGradient/lin-jac"); assert(B && Bt); @@ -123,10 +147,14 @@ Operator& SteadyNSOperator::GetGradient(const Vector &x) const if (direct_solve) { jac_hypre = new HypreParMatrix(MPI_COMM_SELF, sys_glob_size, sys_row_starts, mono_jac); + timer.Stop("GetGradient"); return *jac_hypre; } else + { + timer.Stop("GetGradient"); return *mono_jac; + } } /* @@ -249,27 +277,33 @@ void SteadyNSEQPROM::Mult(const Vector &x, Vector &y) const timer.Start("Mult"); y = 0.0; + timer.Start("Mult/lin"); linearOp->Mult(x, y); + timer.Stop("Mult/lin"); + timer.Start("Mult/nlin"); for (int m = 0; m < numSub; m++) { int midx = midxs[m]; x_comp.MakeRef(const_cast(x), block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); y_comp.MakeRef(y, block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); - timer.Start("Mult_nlin_domain"); hs[m]->AddMult(x_comp, y_comp); - timer.Stop("Mult_nlin_domain"); } + timer.Stop("Mult/nlin"); - if (!itf) return; + if (!itf) + { + timer.Stop("Mult"); + return; + } GetVel(x, x_u); y_u = 0.0; - timer.Start("Mult_nlin_itf"); + timer.Start("Mult/itf"); itf->InterfaceAddMult(x_u, y_u); - timer.Stop("Mult_nlin_itf"); + timer.Stop("Mult/itf"); AddVel(y_u, y); @@ -282,10 +316,13 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const delete jac_mono; delete jac_hypre; + timer.Start("GetGradient/lin"); jac_mono = new SparseMatrix(*linearOp); + timer.Stop("GetGradient/lin"); if (itf) { + timer.Start("GetGradient/itf"); BlockMatrix jac(block_offsets); jac.owns_blocks = true; @@ -301,9 +338,7 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const GetVel(x, x_u); - timer.Start("Grad_nlin_itf"); itf->InterfaceGetGradient(x_u, mats); - timer.Stop("Grad_nlin_itf"); MatrixBlocks u_mats(mats); @@ -313,8 +348,10 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const (*jac_mono) += *tmp; delete tmp; + timer.Stop("GetGradient/itf"); } + timer.Start("GetGradient/nlin"); DenseMatrix *jac_comp; for (int m = 0; m < numSub; m++) { @@ -322,11 +359,10 @@ Operator& SteadyNSEQPROM::GetGradient(const Vector &x) const x_comp.MakeRef(const_cast(x), block_offsets[midx], block_offsets[midx+1] - block_offsets[midx]); // NOTE(kevin): jac_comp is owned by hs[m]. No need of deleting it. - timer.Start("Grad_nlin_domain"); jac_comp = dynamic_cast(&hs[m]->GetGradient(x_comp)); - timer.Stop("Grad_nlin_domain"); jac_mono->AddSubMatrix(*block_idxs[m], *block_idxs[m], *jac_comp); } + timer.Stop("GetGradient/nlin"); jac_mono->Finalize(); timer.Stop("GetGradient"); From 53e2ce371c261b7563845b08cacab1dda631c91f Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 19 Jun 2025 11:23:25 -0700 Subject: [PATCH 18/20] reflecting dylans comments. --- bin/main.cpp | 2 +- include/etc.hpp | 2 +- include/unsteady_ns_solver.hpp | 3 + src/rom_handler.cpp | 3 +- src/unsteady_ns_solver.cpp | 160 ++++++++++----------------------- 5 files changed, 54 insertions(+), 116 deletions(-) diff --git a/bin/main.cpp b/bin/main.cpp index ced806aa..e84e19a3 100644 --- a/bin/main.cpp +++ b/bin/main.cpp @@ -33,7 +33,7 @@ int main(int argc, char *argv[]) else if (mode == "train_rom") TrainROM(MPI_COMM_WORLD); else if (mode == "aux_train_rom") AuxiliaryTrainROM(MPI_COMM_WORLD); else if (mode == "train_eqp") TrainEQP(MPI_COMM_WORLD); - else if (mode == "single_run") double dump = SingleRun(MPI_COMM_WORLD, output_file); + else if (mode == "single_run") SingleRun(MPI_COMM_WORLD, output_file); else if (mode == "print_eqp") PrintEQPCoords(MPI_COMM_WORLD); else { diff --git a/include/etc.hpp b/include/etc.hpp index 5b05d822..5b4ad977 100644 --- a/include/etc.hpp +++ b/include/etc.hpp @@ -53,7 +53,7 @@ class TimeProfiler public: TimeProfiler(MPI_Comm comm_=MPI_COMM_WORLD) - : comm(comm_), timers(0), starts(0) + : comm(comm_) { MPI_Comm_rank(MPI_COMM_WORLD, &rank); calls.SetSize(0); diff --git a/include/unsteady_ns_solver.hpp b/include/unsteady_ns_solver.hpp index d7a30ad1..c579b57a 100644 --- a/include/unsteady_ns_solver.hpp +++ b/include/unsteady_ns_solver.hpp @@ -16,6 +16,9 @@ class UnsteadyNSSolver : public SteadyNSSolver friend class ParameterizedProblem; friend class SteadyNSOperator; +private: + mutable TimeProfiler timer; + protected: // number of timesteps diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index 7700ad15..b231e42b 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -638,12 +638,13 @@ void MFEMROMHandler::NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec printf("Solve ROM.\n"); reduced_sol = new BlockVector(rom_block_offsets); bool use_restart = config.GetOption("solver/use_restart", false); + double amp = config.GetOption("solver/initial_random_perturbation", 1.0e-5); if (use_restart) ProjectGlobalToDomainBasis(U, reduced_sol); else { for (int k = 0; k < reduced_sol->Size(); k++) - (*reduced_sol)(k) = 1.0e-5 * UniformRandom(); + (*reduced_sol)(k) = amp * UniformRandom(); } int maxIter = config.GetOption("solver/max_iter", 100); diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp index d1e0982d..1c570ccf 100644 --- a/src/unsteady_ns_solver.cpp +++ b/src/unsteady_ns_solver.cpp @@ -13,7 +13,7 @@ using namespace mfem; */ UnsteadyNSSolver::UnsteadyNSSolver() - : SteadyNSSolver() + : SteadyNSSolver(), timer() { nonlinear_mode = true; @@ -88,10 +88,7 @@ bool UnsteadyNSSolver::Solve(SampleGenerator *sample_generator) std::string restart_file, file_fmt; file_fmt = "%s/%s_%08d.h5"; - StopWatch timers; - for (int k = 0; k < 10; k++) - times[k] = 0.0; - timers.Start(); + timer.Start("Solve/init"); SetupInitialCondition(initial_step, time); @@ -104,16 +101,14 @@ bool UnsteadyNSSolver::Solve(SampleGenerator *sample_generator) SaveVisualization(0, time); - timers.Stop(); - times[0] = timers.RealTime(); - timers.Clear(); + timer.Stop("Solve/init"); double cfl = 0.0; for (int step = initial_step; step < nt; step++) { Step(time, step); - timers.Start(); + timer.Start("Solve/save_sol"); cfl = ComputeCFL(dt); SanityCheck(step); @@ -137,27 +132,16 @@ bool UnsteadyNSSolver::Solve(SampleGenerator *sample_generator) ((step+1) > bootstrap) && (((step+1) % sample_interval) == 0)) SaveSnapshots(sample_generator); - timers.Stop(); - times[9] += timers.RealTime(); - timers.Clear(); + timer.Stop("Solve/save_sol"); } - timers.Start(); + timer.Start("Solve/final"); SortBySubdomains(*U_step, *U); - timers.Stop(); - times[9] += timers.RealTime(); - timers.Clear(); + timer.Stop("Solve/final"); - double total = 0.0; - for (int k = 0; k < 10; k++) - total += times[k]; - - printf("FOM Solve time: %.3e\n", total); - for (int k = 0; k < 10; k++) - printf("%.3e\t", times[k]); - printf("\n"); + timer.Print("UnsteadyNSSolver::Solve", true); return converged; } @@ -238,8 +222,7 @@ void UnsteadyNSSolver::SetupInitialCondition(int &initial_step, double &time) void UnsteadyNSSolver::Step(double &time, int step) { - StopWatch timers; - timers.Start(); + timer.Start("Step/copy_u"); /* set time for forcing/boundary. At this point, time remains at the previous timestep. */ SetTime(time); @@ -247,58 +230,44 @@ void UnsteadyNSSolver::Step(double &time, int step) /* copy velocity */ u1 = U_stepview->GetBlock(0); - timers.Stop(); - times[1] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("Step/copy_u"); + timer.Start("Step/advect_dom"); /* evaluate nonlinear advection at previous time step */ Hop->Mult(u1, Cu1); - timers.Stop(); - times[2] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("Step/advect_dom"); + timer.Start("Step/advect_itf"); nl_itf->InterfaceAddMult(u1, Cu1); - timers.Stop(); - times[3] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("Step/advect_itf"); + timer.Start("Step/rhs_sort"); /* Base right-hand side for boundary conditions and forcing */ SortByVariables(*RHS, *RHS_step); - timers.Stop(); - times[4] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("Step/rhs_sort"); + timer.Start("Step/add_advect"); /* Add nonlinear convection */ RHS_stepview->GetBlock(0).Add(-ab1, Cu1); - timers.Stop(); - times[5] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("Step/add_advect"); + timer.Start("Step/add_udot"); /* Add time derivative term */ // TODO: extend for high order bdf schemes massMat->AddMult(u1, RHS_stepview->GetBlock(0), -bd1 / dt); - timers.Stop(); - times[6] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("Step/add_udot"); + timer.Start("Step/solve_system"); /* Solve for the next step */ mumps->Mult(*RHS_step, *U_step); - timers.Stop(); - times[7] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("Step/solve_system"); + timer.Start("Step/pres_scale"); /* remove pressure scalar if all dirichlet bc */ if (!pres_dbc) @@ -308,9 +277,7 @@ void UnsteadyNSSolver::Step(double &time, int step) U_stepview->GetBlock(1) -= p_const; } - timers.Stop(); - times[8] += timers.RealTime(); - timers.Clear(); + timer.Stop("Step/pres_scale"); time += dt; } @@ -504,11 +471,7 @@ void UnsteadyNSSolver::SolveROM() { assert(rom_handler->GetOrdering() == ROMOrderBy::VARIABLE); - StopWatch timers; - for (int k = 0; k < 10; k++) - times[k] = 0.0; - - timers.Start(); + timer.Start("SolveROM/init"); int initial_step = 0; double time = 0.0; @@ -561,77 +524,59 @@ void UnsteadyNSSolver::SolveROM() } } - timers.Stop(); - times[0] = timers.RealTime(); - timers.Clear(); + timer.Stop("SolveROM/init"); for (int step = initial_step; step < nt; step++) { - timers.Start(); + timer.Start("SolveROM/set_time"); /* set time for forcing/boundary. At this point, time remains at the previous timestep. */ SetTime(time); - timers.Stop(); - times[0] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("SolveROM/set_time"); + timer.Start("SolveROM/copy_u"); /* copy velocity */ u1 = rsol_view->GetBlock(0); - timers.Stop(); - times[1] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("SolveROM/copy_u"); + timer.Start("SolveROM/advect_dom"); /* evaluate nonlinear advection at previous time step */ Hop->Mult(u1, Cu1); - timers.Stop(); - times[2] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("SolveROM/advect_dom"); + timer.Start("SolveROM/advect_itf"); itf_eqp->InterfaceAddMult(u1, Cu1); - timers.Stop(); - times[3] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("SolveROM/advect_itf"); + timer.Start("SolveROM/get_rhs"); /* Base right-hand side for boundary conditions and forcing */ reduced_rhs = *(rom_handler->GetReducedRHS()); - timers.Stop(); - times[4] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("SolveROM/get_rhs"); + timer.Start("SolveROM/add_advect"); /* Add nonlinear convection */ rrhs_view->GetBlock(0).Add(-ab1, Cu1); - timers.Stop(); - times[5] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("SolveROM/add_advect"); + timer.Start("SolveROM/add_udot"); /* Add time derivative term */ // TODO: extend for high order bdf schemes rom_mass->AddMult(u1, rrhs_view->GetBlock(0), -bd1 / dt); - timers.Stop(); - times[6] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("SolveROM/add_udot"); + timer.Start("SolveROM/solve_system"); /* Solve for the next step */ rom_handler->Solve(reduced_rhs, *reduced_sol); - timers.Stop(); - times[7] += timers.RealTime(); - timers.Clear(); - timers.Start(); + timer.Stop("SolveROM/solve_system"); + timer.Start("SolveROM/pres_scale"); /* remove pressure scalar if all dirichlet bc */ if (!pres_dbc) @@ -641,29 +586,18 @@ void UnsteadyNSSolver::SolveROM() rsol_view->GetBlock(1).Add(-p_const, rom_ones); } - timers.Stop(); - times[8] += timers.RealTime(); - timers.Clear(); + timer.Stop("SolveROM/pres_scale"); time += dt; } - timers.Start(); + timer.Start("SolveROM/liftup"); rom_handler->LiftUpGlobal(*reduced_sol, *U); - timers.Stop(); - times[9] += timers.RealTime(); - timers.Clear(); - - double total = 0.0; - for (int k = 0; k < 10; k++) - total += times[k]; + timer.Stop("SolveROM/liftup"); - printf("ROM Solve time: %.3e\n", total); - for (int k = 0; k < 10; k++) - printf("%.3e\t", times[k]); - printf("\n"); + timer.Print("UnsteadyNSSolver::SolveROM", true); delete rsol_view; delete reduced_sol; From 36775fdf2912026edd0ce8d1cfa12ceb559da1e1 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 19 Jun 2025 12:02:56 -0700 Subject: [PATCH 19/20] bug/warning fixes. --- include/etc.hpp | 6 +++--- src/rom_handler.cpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/include/etc.hpp b/include/etc.hpp index 5b4ad977..25ee5437 100644 --- a/include/etc.hpp +++ b/include/etc.hpp @@ -122,11 +122,11 @@ class TimeProfiler if (rank == 0) { - printf((title + "\n").c_str()); + printf("%s", (title + "\n").c_str()); std::string line = std::string(100, '='); line += "\n"; - printf(line.c_str()); + printf("%s", line.c_str()); printf("%20s\t%20s\t%20s\t%20s\n", "Function", "Total time", "Calls", "Time per call"); for (int k = 0; k < nfunc; k++) { @@ -134,7 +134,7 @@ class TimeProfiler } if (compute_sum) printf("%20s\t%20.5e\t%20d\t%20.5e\n", "Total time", total, 1, total); - printf(line.c_str()); + printf("%s", line.c_str()); } } }; diff --git a/src/rom_handler.cpp b/src/rom_handler.cpp index b231e42b..3edc51f3 100644 --- a/src/rom_handler.cpp +++ b/src/rom_handler.cpp @@ -638,7 +638,7 @@ void MFEMROMHandler::NonlinearSolve(Operator &oper, BlockVector* U, Solver *prec printf("Solve ROM.\n"); reduced_sol = new BlockVector(rom_block_offsets); bool use_restart = config.GetOption("solver/use_restart", false); - double amp = config.GetOption("solver/initial_random_perturbation", 1.0e-5); + double amp = config.GetOption("solver/initial_random_perturbation", 1.0e-5); if (use_restart) ProjectGlobalToDomainBasis(U, reduced_sol); else From 061dfb61607f0c7aa028259b3eb3fa519a0948f5 Mon Sep 17 00:00:00 2001 From: Kevin Chung Date: Thu, 19 Jun 2025 14:15:49 -0700 Subject: [PATCH 20/20] TimeProfile now possesses its own title. At its destruction, it prints with its title. --- include/etc.hpp | 17 +++++++++++++---- include/interface_form.hpp | 3 +-- include/interfaceinteg.hpp | 4 ++-- include/rom_interfaceform.hpp | 3 --- include/steady_ns_solver.hpp | 6 +----- src/interface_form.cpp | 6 ++---- src/interfaceinteg.cpp | 6 ------ src/rom_interfaceform.cpp | 7 +++---- src/rom_nonlinearform.cpp | 5 +---- src/steady_ns_solver.cpp | 9 ++++----- src/unsteady_ns_solver.cpp | 2 +- 11 files changed, 28 insertions(+), 40 deletions(-) diff --git a/include/etc.hpp b/include/etc.hpp index 25ee5437..9aad33e9 100644 --- a/include/etc.hpp +++ b/include/etc.hpp @@ -6,6 +6,7 @@ #define ETC_HPP #include "mfem.hpp" +#include "input_parser.hpp" using namespace mfem; using namespace std; @@ -41,6 +42,8 @@ std::string string_format( const std::string& format, Args ... args ) class TimeProfiler { private: + std::string title = ""; + Array timers; Array calls; @@ -52,8 +55,8 @@ class TimeProfiler int rank; public: - TimeProfiler(MPI_Comm comm_=MPI_COMM_WORLD) - : comm(comm_) + TimeProfiler(const std::string &title_, MPI_Comm comm_=MPI_COMM_WORLD) + : title(title_), comm(comm_) { MPI_Comm_rank(MPI_COMM_WORLD, &rank); calls.SetSize(0); @@ -63,9 +66,15 @@ class TimeProfiler virtual ~TimeProfiler() { + std::string option = "time_profile/" + title; + if (config.GetOption(option, false)) + Print(title); DeletePointers(timers); } + void SetTitle(const std::string &title_) + { title = title_; } + void Start(const std::string &name) { if (!indices.count(name)) @@ -106,7 +115,7 @@ class TimeProfiler starts[idx] = false; } - void Print(const std::string &title, const bool compute_sum=false) + void Print(const std::string &title_, const bool compute_sum=false) { int nfunc = timers.Size(); for (int k = 0; k < nfunc; k++) @@ -122,7 +131,7 @@ class TimeProfiler if (rank == 0) { - printf("%s", (title + "\n").c_str()); + printf("%s", (title_ + "\n").c_str()); std::string line = std::string(100, '='); line += "\n"; diff --git a/include/interface_form.hpp b/include/interface_form.hpp index 9a2185b9..24ffdf6a 100644 --- a/include/interface_form.hpp +++ b/include/interface_form.hpp @@ -15,10 +15,9 @@ namespace mfem // TODO(kevin): inherits or cherry-pick ROMNonlinearForm for hyper-reduction. class InterfaceForm { -private: +protected: mutable TimeProfiler timer; -protected: int numSub = -1; int skip_zeros = 1; diff --git a/include/interfaceinteg.hpp b/include/interfaceinteg.hpp index 6b8372b5..4f251372 100644 --- a/include/interfaceinteg.hpp +++ b/include/interfaceinteg.hpp @@ -303,9 +303,9 @@ class DGLaxFriedrichsFluxIntegrator : public InterfaceNonlinearFormIntegrator public: DGLaxFriedrichsFluxIntegrator(Coefficient &q, VectorCoefficient *ud = NULL, const IntegrationRule *ir = NULL) - : InterfaceNonlinearFormIntegrator(ir), Q(&q), UD(ud) {} + : InterfaceNonlinearFormIntegrator(ir), Q(&q), UD(ud), timer("DGLaxFriedrichsFluxIntegrator") {} - virtual ~DGLaxFriedrichsFluxIntegrator(); + virtual ~DGLaxFriedrichsFluxIntegrator() {}; void AssembleFaceVector(const FiniteElement &el1, const FiniteElement &el2, diff --git a/include/rom_interfaceform.hpp b/include/rom_interfaceform.hpp index 6c5d4fef..4193c71c 100644 --- a/include/rom_interfaceform.hpp +++ b/include/rom_interfaceform.hpp @@ -17,9 +17,6 @@ namespace mfem class ROMInterfaceForm : public InterfaceForm { -private: - mutable TimeProfiler timer; - protected: const int numPorts; const int numRefPorts; diff --git a/include/steady_ns_solver.hpp b/include/steady_ns_solver.hpp index 1fe70c98..c5bb02c8 100644 --- a/include/steady_ns_solver.hpp +++ b/include/steady_ns_solver.hpp @@ -118,11 +118,7 @@ class SteadyNSEQPROM : public SteadyNSROM SteadyNSEQPROM(ROMHandlerBase *rom_handler, Array &hs_, ROMInterfaceForm *itf_, const bool direct_solve_=true); - virtual ~SteadyNSEQPROM() - { - if (config.GetOption("time_profile/SteadyNSEQPROM", false)) - timer.Print("SteadyNSEQPROM"); - } + virtual ~SteadyNSEQPROM() {} virtual void Mult(const Vector &x, Vector &y) const; virtual Operator &GetGradient(const Vector &x) const; diff --git a/src/interface_form.cpp b/src/interface_form.cpp index 4bf1cc12..5eae5853 100644 --- a/src/interface_form.cpp +++ b/src/interface_form.cpp @@ -12,7 +12,8 @@ namespace mfem InterfaceForm::InterfaceForm( Array &meshes_, Array &fes_, TopologyHandler *topol_) - : meshes(meshes_), fes(fes_), topol_handler(topol_), numSub(meshes_.Size()), timer() + : meshes(meshes_), fes(fes_), topol_handler(topol_), + numSub(meshes_.Size()), timer("InterfaceForm") { assert(fes_.Size() == numSub); @@ -26,9 +27,6 @@ InterfaceForm::InterfaceForm( InterfaceForm::~InterfaceForm() { DeletePointers(fnfi); - - if (config.GetOption("time_profile/InterfaceForm", false)) - timer.Print("InterfaceForm::InterfaceAddMult"); } void InterfaceForm::AssembleInterfaceMatrices(Array2D &mats) const diff --git a/src/interfaceinteg.cpp b/src/interfaceinteg.cpp index dacd34d3..e791a07c 100644 --- a/src/interfaceinteg.cpp +++ b/src/interfaceinteg.cpp @@ -1308,12 +1308,6 @@ void InterfaceDGElasticityIntegrator::AssembleBlock( DGLaxFriedrichsFluxIntegrator */ -DGLaxFriedrichsFluxIntegrator::~DGLaxFriedrichsFluxIntegrator() -{ - if (config.GetOption("time_profile/DGLaxFriedrichsFluxIntegrator", false)) - timer.Print("DGLaxFriedrichsFluxIntegrator"); -} - void DGLaxFriedrichsFluxIntegrator::ComputeFluxDotN( const Vector &u1, const Vector &u2, const Vector &nor, const bool &eval2, Vector &flux) diff --git a/src/rom_interfaceform.cpp b/src/rom_interfaceform.cpp index 06359833..13e5dca7 100644 --- a/src/rom_interfaceform.cpp +++ b/src/rom_interfaceform.cpp @@ -14,7 +14,7 @@ namespace mfem ROMInterfaceForm::ROMInterfaceForm( Array &meshes_, Array &fes_, Array &comp_fes_, TopologyHandler *topol_) : InterfaceForm(meshes_, fes_, topol_), numPorts(topol_->GetNumPorts()), numRefPorts(topol_->GetNumRefPorts()), - num_comp(topol_->GetNumComponents()), comp_fes(comp_fes_), timer() + num_comp(topol_->GetNumComponents()), comp_fes(comp_fes_) { comp_basis.SetSize(num_comp); comp_basis = NULL; @@ -32,14 +32,13 @@ ROMInterfaceForm::ROMInterfaceForm( else if (nnls_str == "linf") nnls_criterion = CAROM::NNLS_termination::LINF; else mfem_error("ROMNonlinearForm: unknown NNLS criterion!\n"); + + timer.SetTitle("ROMInterfaceForm"); } ROMInterfaceForm::~ROMInterfaceForm() { DeletePointers(fnfi_ref_sample); - - if (config.GetOption("time_profile/ROMInterfaceForm", false)) - timer.Print("ROMInterfaceForm"); } void ROMInterfaceForm::SetBasisAtComponent(const int c, DenseMatrix &basis_, const int offset) diff --git a/src/rom_nonlinearform.cpp b/src/rom_nonlinearform.cpp index 75342753..9b7bf00d 100644 --- a/src/rom_nonlinearform.cpp +++ b/src/rom_nonlinearform.cpp @@ -12,7 +12,7 @@ namespace mfem { ROMNonlinearForm::ROMNonlinearForm(const int num_basis, FiniteElementSpace *f, const bool reference_) - : NonlinearForm(f), reference(reference_), timer() + : NonlinearForm(f), reference(reference_), timer("ROMNonlinearForm") { height = width = num_basis; @@ -41,9 +41,6 @@ ROMNonlinearForm::~ROMNonlinearForm() DeletePointers(fnfi_sample); DeletePointers(bfnfi_sample); } - - if (config.GetOption("time_profile/ROMNonlinearForm", false)) - timer.Print("ROMNonlinearForm"); } void ROMNonlinearForm::Mult(const Vector &x, Vector &y) const diff --git a/src/steady_ns_solver.cpp b/src/steady_ns_solver.cpp index 852332f7..7a3ae61a 100644 --- a/src/steady_ns_solver.cpp +++ b/src/steady_ns_solver.cpp @@ -24,7 +24,8 @@ SteadyNSOperator::SteadyNSOperator( Array &u_offsets_, const bool direct_solve_) : Operator(linearOp_->Height(), linearOp_->Width()), linearOp(linearOp_), hs(hs_), nl_itf(nl_itf_), u_offsets(u_offsets_), direct_solve(direct_solve_), - M(&(linearOp_->GetBlock(0, 0))), Bt(&(linearOp_->GetBlock(0, 1))), B(&(linearOp_->GetBlock(1, 0))) + M(&(linearOp_->GetBlock(0, 0))), Bt(&(linearOp_->GetBlock(0, 1))), B(&(linearOp_->GetBlock(1, 0))), + timer("SteadyNSOperator") { vblock_offsets.SetSize(3); vblock_offsets[0] = 0; @@ -57,9 +58,6 @@ SteadyNSOperator::~SteadyNSOperator() delete uu_mono; delete mono_jac; delete jac_hypre; - - if (config.GetOption("time_profile/SteadyNSOperator", false)) - timer.Print("SteadyNSOperator"); } void SteadyNSOperator::Mult(const Vector &x, Vector &y) const @@ -250,7 +248,8 @@ Operator& SteadyNSTensorROM::GetGradient(const Vector &x) const SteadyNSEQPROM::SteadyNSEQPROM( ROMHandlerBase *rom_handler, Array &hs_, ROMInterfaceForm *itf_, const bool direct_solve_) - : SteadyNSROM(hs_.Size(), rom_handler, direct_solve_), hs(hs_), itf(itf_), timer() + : SteadyNSROM(hs_.Size(), rom_handler, direct_solve_), hs(hs_), itf(itf_), + timer("SteadyNSEQPROM") { if (!separate_variable) { diff --git a/src/unsteady_ns_solver.cpp b/src/unsteady_ns_solver.cpp index 1c570ccf..bda721fc 100644 --- a/src/unsteady_ns_solver.cpp +++ b/src/unsteady_ns_solver.cpp @@ -13,7 +13,7 @@ using namespace mfem; */ UnsteadyNSSolver::UnsteadyNSSolver() - : SteadyNSSolver(), timer() + : SteadyNSSolver(), timer("UnsteadyNSSolver") { nonlinear_mode = true;