diff --git a/examples/prom/maxwell_global_rom.cpp b/examples/prom/maxwell_global_rom.cpp index 10b71fa3b..8019185ad 100644 --- a/examples/prom/maxwell_global_rom.cpp +++ b/examples/prom/maxwell_global_rom.cpp @@ -64,81 +64,80 @@ int dim; int main(int argc, char *argv[]) { - // 1. Initialize MPI and HYPRE - Mpi::Init(argc, argv); - int num_procs = Mpi::WorldSize(); - int myid = Mpi::WorldRank(); - Hypre::Init(); + // 1. Initialize MPI and HYPRE + Mpi::Init(argc, argv); + int num_procs = Mpi::WorldSize(); + int myid = Mpi::WorldRank(); - // 2. Parse command-line options. + // 2. Parse command-line options. // const char *mesh_file = "../data/beam-tet.mesh"; - const char *mesh_file = "../data/star.mesh"; - int order = 1; - bool static_cond = false; - bool pa = false; - const char *device_config = "cpu"; - bool visualization = true; - bool visit = true; - bool fom = false; - bool offline = false; - bool merge = false; - bool online = false; - int precision = 8; - int id = 0; - int nsets = 0; - double coef = 1.0; + const char *mesh_file = "../data/star.mesh"; + int order = 1; + bool static_cond = false; + bool pa = false; + const char *device_config = "cpu"; + bool visualization = true; + bool visit = true; + bool fom = false; + bool offline = false; + bool merge = false; + bool online = false; + int precision = 8; + int id = 0; + int nsets = 0; + double coef = 1.0; #ifdef MFEM_USE_AMGX - bool useAmgX = false; + bool useAmgX = false; #endif - - - OptionsParser args(argc, argv); - args.AddOption(&mesh_file, "-m", "--mesh", - "Mesh file to use."); - args.AddOption(&order, "-o", "--order", - "Finite element order (polynomial degree)."); - args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact" - " solution."); - args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc", - "--no-static-condensation", "Enable static condensation."); - args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa", - "--no-partial-assembly", "Enable Partial Assembly."); - args.AddOption(&device_config, "-d", "--device", - "Device configuration string, see Device::Configure()."); - args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", - "--no-visualization", - "Enable or disable GLVis visualization."); + + + OptionsParser args(argc, argv); + args.AddOption(&mesh_file, "-m", "--mesh", + "Mesh file to use."); + args.AddOption(&order, "-o", "--order", + "Finite element order (polynomial degree)."); + args.AddOption(&freq, "-f", "--frequency", "Set the frequency for the exact" + " solution."); + args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc", + "--no-static-condensation", "Enable static condensation."); + args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa", + "--no-partial-assembly", "Enable Partial Assembly."); + args.AddOption(&device_config, "-d", "--device", + "Device configuration string, see Device::Configure()."); + args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", + "--no-visualization", + "Enable or disable GLVis visualization."); #ifdef MFEM_USE_AMGX - args.AddOption(&useAmgX, "-amgx", "--useAmgX", "-no-amgx", - "--no-useAmgX", - "Enable or disable AmgX in MatrixFreeAMS."); + args.AddOption(&useAmgX, "-amgx", "--useAmgX", "-no-amgx", + "--no-useAmgX", + "Enable or disable AmgX in MatrixFreeAMS."); #endif - args.AddOption(&fom, "-fom", "--fom", "-no-fom", "--no-fom", - "Enable or disable the fom phase."); - args.AddOption(&offline, "-offline", "--offline", "-no-offline","--no-offline", - "Enable or disable the offline phase."); - args.AddOption(&online, "-online", "--online", "-no-online", "--no-online", - "Enable or disable the online phase."); - args.AddOption(&merge, "-merge", "--merge", "-no-merge", "--no-merge", - "Enable or disable the merge phase."); - args.AddOption(&id, "-id", "--id", "Parametric id"); - args.AddOption(&nsets, "-ns", "--nset", "Number of parametric snapshot sets"); - - args.Parse(); - if (!args.Good()) - { - if (myid == 0) - { - args.PrintUsage(cout); - } - return 1; - } - if (myid == 0) - { - args.PrintOptions(cout); - } - kappa = freq * M_PI; - + args.AddOption(&fom, "-fom", "--fom", "-no-fom", "--no-fom", + "Enable or disable the fom phase."); + args.AddOption(&offline, "-offline", "--offline", "-no-offline","--no-offline", + "Enable or disable the offline phase."); + args.AddOption(&online, "-online", "--online", "-no-online", "--no-online", + "Enable or disable the online phase."); + args.AddOption(&merge, "-merge", "--merge", "-no-merge", "--no-merge", + "Enable or disable the merge phase."); + args.AddOption(&id, "-id", "--id", "Parametric id"); + args.AddOption(&nsets, "-ns", "--nset", "Number of parametric snapshot sets"); + + args.Parse(); + if (!args.Good()) + { + if (myid == 0) + { + args.PrintUsage(cout); + } + return 1; + } + if (myid == 0) + { + args.PrintOptions(cout); + } + kappa = freq * M_PI; + if (fom) { MFEM_VERIFY(fom && !offline && !online @@ -151,65 +150,67 @@ int main(int argc, char *argv[]) MFEM_VERIFY(check, "only one of offline, merge, or online must be true!"); } - // 3. Enable hardware devices such as GPUs, and programming models such as - // CUDA, OCCA, RAJA and OpenMP based on command line options. - Device device(device_config); - if (myid == 0) { device.Print(); } - - // 4. Read the (serial) mesh from the given mesh file on all processors. We - // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface - // and volume meshes with the same code. - Mesh *mesh = new Mesh(mesh_file, 1, 1); - dim = mesh->Dimension(); - int sdim = mesh->SpaceDimension(); - - // 5. Refine the serial mesh on all processors to increase the resolution. In - // this example we do 'ref_levels' of uniform refinement. We choose - // 'ref_levels' to be the largest number that gives a final mesh with no - // more than 1,000 elements. - { - int ref_levels = (int)floor(log(1000./mesh->GetNE())/log(2.)/dim); - for (int l = 0; l < ref_levels; l++) - { - mesh->UniformRefinement(); - } - } - - // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine - // this mesh further in parallel to increase the resolution. Once the - // parallel mesh is defined, the serial mesh can be deleted. - ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh); - delete mesh; - { - int par_ref_levels = 2; - for (int l = 0; l < par_ref_levels; l++) - { - pmesh->UniformRefinement(); - } - } - - // 7. Define a parallel finite element space on the parallel mesh. Here we - // use the Nedelec finite elements of the specified order. - FiniteElementCollection *fec = new ND_FECollection(order, dim); - ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec); - HYPRE_BigInt size = fespace->GlobalTrueVSize(); - if (myid == 0) - { - cout << "Number of finite element unknowns: " << size << endl; - } - - // 8. Determine the list of true (i.e. parallel conforming) essential - // boundary dofs. In this example, the boundary conditions are defined - // by marking all the boundary attributes from the mesh as essential - // (Dirichlet) and converting them to a list of true dofs. - Array ess_tdof_list; - Array ess_bdr; - if (pmesh->bdr_attributes.Size()) - { - ess_bdr.SetSize(pmesh->bdr_attributes.Max()); - ess_bdr = 1; - fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list); - } + // 3. Enable hardware devices such as GPUs, and programming models such as + // CUDA, OCCA, RAJA and OpenMP based on command line options. + Device device(device_config); + if (myid == 0) { + device.Print(); + } + + // 4. Read the (serial) mesh from the given mesh file on all processors. We + // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface + // and volume meshes with the same code. + Mesh *mesh = new Mesh(mesh_file, 1, 1); + dim = mesh->Dimension(); + int sdim = mesh->SpaceDimension(); + + // 5. Refine the serial mesh on all processors to increase the resolution. In + // this example we do 'ref_levels' of uniform refinement. We choose + // 'ref_levels' to be the largest number that gives a final mesh with no + // more than 1,000 elements. + { + int ref_levels = (int)floor(log(1000./mesh->GetNE())/log(2.)/dim); + for (int l = 0; l < ref_levels; l++) + { + mesh->UniformRefinement(); + } + } + + // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine + // this mesh further in parallel to increase the resolution. Once the + // parallel mesh is defined, the serial mesh can be deleted. + ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh); + delete mesh; + { + int par_ref_levels = 2; + for (int l = 0; l < par_ref_levels; l++) + { + pmesh->UniformRefinement(); + } + } + + // 7. Define a parallel finite element space on the parallel mesh. Here we + // use the Nedelec finite elements of the specified order. + FiniteElementCollection *fec = new ND_FECollection(order, dim); + ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec); + HYPRE_BigInt size = fespace->GlobalTrueVSize(); + if (myid == 0) + { + cout << "Number of finite element unknowns: " << size << endl; + } + + // 8. Determine the list of true (i.e. parallel conforming) essential + // boundary dofs. In this example, the boundary conditions are defined + // by marking all the boundary attributes from the mesh as essential + // (Dirichlet) and converting them to a list of true dofs. + Array ess_tdof_list; + Array ess_bdr; + if (pmesh->bdr_attributes.Size()) + { + ess_bdr.SetSize(pmesh->bdr_attributes.Max()); + ess_bdr = 1; + fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list); + } // 9. Initiate ROM related variables int max_num_snapshots = 100; @@ -222,7 +223,7 @@ int main(int argc, char *argv[]) CAROM::BasisGenerator *generator; int numRowRB, numColumnRB; StopWatch solveTimer, assembleTimer, mergeTimer; - + // 10. Set BasisGenerator if offline if (offline) { @@ -230,7 +231,7 @@ int main(int argc, char *argv[]) update_right_SV); generator = new CAROM::BasisGenerator(*options, isIncremental, basisFileName); } - + // 11. The merge phase if (merge) { @@ -258,281 +259,290 @@ int main(int argc, char *argv[]) return 0; } - // 12. Set up the parallel linear form b(.) which corresponds to the - // right-hand side of the FEM linear system, which in this case is - // (f,phi_i) where f is given by the function f_exact and phi_i are the - // basis functions in the finite element fespace. - VectorFunctionCoefficient f(sdim, f_exact); - ParLinearForm *b = new ParLinearForm(fespace); - b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f)); - b->Assemble(); - - // 13. Define the solution vector x as a parallel finite element grid function - // corresponding to fespace. Initialize x by projecting the exact - // solution. Note that only values from the boundary edges will be used - // when eliminating the non-homogeneous boundary condition to modify the - // r.h.s. vector b. - ParGridFunction x(fespace); - VectorFunctionCoefficient E(sdim, E_exact); - x.ProjectCoefficient(E); - - // 14. Set up the parallel bilinear form corresponding to the EM diffusion - // operator curl muinv curl + sigma I, by adding the curl-curl and the - // mass domain integrators. - Coefficient *muinv = new ConstantCoefficient(1.0); - Coefficient *sigma = new ConstantCoefficient(1.0); - ParBilinearForm *a = new ParBilinearForm(fespace); - if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); } - a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv)); - a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma)); - - // 15. Assemble the parallel bilinear form and the corresponding linear - // system, applying any necessary transformations such as: parallel - // assembly, eliminating boundary conditions, applying conforming - // constraints for non-conforming AMR, static condensation, etc. - if (static_cond) { a->EnableStaticCondensation(); } - a->Assemble(); - - OperatorPtr A; - Vector B, X; - a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B); - assembleTimer.Stop(); - - // 16. Offline Phase: Solve the system AX=B using PCG with an AMS preconditioner. - if(fom || offline) - { - // 17. Solve the full order linear system A X = B - if (pa) - { + // 12. Set up the parallel linear form b(.) which corresponds to the + // right-hand side of the FEM linear system, which in this case is + // (f,phi_i) where f is given by the function f_exact and phi_i are the + // basis functions in the finite element fespace. + VectorFunctionCoefficient f(sdim, f_exact); + ParLinearForm *b = new ParLinearForm(fespace); + b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f)); + b->Assemble(); + + // 13. Define the solution vector x as a parallel finite element grid function + // corresponding to fespace. Initialize x by projecting the exact + // solution. Note that only values from the boundary edges will be used + // when eliminating the non-homogeneous boundary condition to modify the + // r.h.s. vector b. + ParGridFunction x(fespace); + VectorFunctionCoefficient E(sdim, E_exact); + x.ProjectCoefficient(E); + + // 14. Set up the parallel bilinear form corresponding to the EM diffusion + // operator curl muinv curl + sigma I, by adding the curl-curl and the + // mass domain integrators. + Coefficient *muinv = new ConstantCoefficient(1.0); + Coefficient *sigma = new ConstantCoefficient(1.0); + ParBilinearForm *a = new ParBilinearForm(fespace); + if (pa) { + a->SetAssemblyLevel(AssemblyLevel::PARTIAL); + } + a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv)); + a->AddDomainIntegrator(new VectorFEMassIntegrator(*sigma)); + + // 15. Assemble the parallel bilinear form and the corresponding linear + // system, applying any necessary transformations such as: parallel + // assembly, eliminating boundary conditions, applying conforming + // constraints for non-conforming AMR, static condensation, etc. + if (static_cond) { + a->EnableStaticCondensation(); + } + a->Assemble(); + + OperatorPtr A; + Vector B, X; + a->FormLinearSystem(ess_tdof_list, x, *b, A, X, B); + assembleTimer.Stop(); + + // 16. Offline Phase: Solve the system AX=B using PCG with an AMS preconditioner. + if(fom || offline) + { + // 17. Solve the full order linear system A X = B + if (pa) + { #ifdef MFEM_USE_AMGX - MatrixFreeAMS ams(*a, *A, *fespace, muinv, sigma, NULL, ess_bdr, useAmgX); + MatrixFreeAMS ams(*a, *A, *fespace, muinv, sigma, NULL, ess_bdr, useAmgX); #else - MatrixFreeAMS ams(*a, *A, *fespace, muinv, sigma, NULL, ess_bdr); + MatrixFreeAMS ams(*a, *A, *fespace, muinv, sigma, NULL, ess_bdr); #endif - CGSolver cg(MPI_COMM_WORLD); - cg.SetRelTol(1e-12); - cg.SetMaxIter(1000); - cg.SetPrintLevel(1); - cg.SetOperator(*A); - cg.SetPreconditioner(ams); - solveTimer.Start(); - cg.Mult(B, X); - solveTimer.Stop(); - } - else - { - if (myid == 0) - { - cout << "Size of linear system: " - << A.As()->GetGlobalNumRows() << endl; - } - - ParFiniteElementSpace *prec_fespace = - (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace); - HypreAMS ams(*A.As(), prec_fespace); - HyprePCG pcg(*A.As()); - pcg.SetTol(1e-12); - pcg.SetMaxIter(500); - pcg.SetPrintLevel(2); - pcg.SetPreconditioner(ams); - solveTimer.Start(); - pcg.Mult(B, X); - solveTimer.Stop(); - } - // 18. take and write snapshot for ROM - if (offline) - { - bool addSample = generator->takeSample(X.GetData(), 0.0, 0.01); - generator->writeSnapshot(); - delete generator; - delete options; - } - } + CGSolver cg(MPI_COMM_WORLD); + cg.SetRelTol(1e-12); + cg.SetMaxIter(1000); + cg.SetPrintLevel(1); + cg.SetOperator(*A); + cg.SetPreconditioner(ams); + solveTimer.Start(); + cg.Mult(B, X); + solveTimer.Stop(); + } + else + { + if (myid == 0) + { + cout << "Size of linear system: " + << A.As()->GetGlobalNumRows() << endl; + } + + ParFiniteElementSpace *prec_fespace = + (a->StaticCondensationIsEnabled() ? a->SCParFESpace() : fespace); + HypreAMS ams(*A.As(), prec_fespace); + HyprePCG pcg(*A.As()); + pcg.SetTol(1e-12); + pcg.SetMaxIter(500); + pcg.SetPrintLevel(2); + pcg.SetPreconditioner(ams); + solveTimer.Start(); + pcg.Mult(B, X); + solveTimer.Stop(); + } + // 18. take and write snapshot for ROM + if (offline) + { + bool addSample = generator->takeSample(X.GetData(), 0.0, 0.01); + generator->writeSnapshot(); + delete generator; + delete options; + } + } // 19. The online phase - if (online) { - // 20. read the reduced basis - assembleTimer.Start(); - CAROM::BasisReader reader(basisName); - spatialbasis = reader.getSpatialBasis(0.0); - numRowRB = spatialbasis->numRows(); - numColumnRB = spatialbasis->numColumns(); - if (myid == 0) printf("spatial basis dimension is %d x %d\n", numRowRB, - numColumnRB); - - // libROM stores the matrix row-wise, so wrapping as a DenseMatrix in MFEM means it is transposed. - DenseMatrix *reducedBasisT = new DenseMatrix(spatialbasis->getData(), - numColumnRB, numRowRB); - - // 21. form inverse ROM operator - CAROM::Matrix invReducedA(numColumnRB, numColumnRB, false); - ComputeCtAB( *A.As(), *spatialbasis, *spatialbasis, invReducedA); - invReducedA.inverse(); - - CAROM::Vector B_carom(B.GetData(), B.Size(), true, false); - CAROM::Vector X_carom(X.GetData(), X.Size(), true, false); - CAROM::Vector *reducedRHS = spatialbasis->transposeMult(&B_carom); - CAROM::Vector reducedSol(numColumnRB, false); - assembleTimer.Stop(); - - // 22. solve ROM - solveTimer.Start(); - invReducedA.mult(*reducedRHS, reducedSol); - solveTimer.Stop(); - - // 23. reconstruct FOM state - spatialbasis->mult(reducedSol, X_carom); - delete spatialbasis; - delete reducedRHS; - } - - // 24. Recover the parallel grid function corresponding to X. This is the - // local finite element solution on each processor. - a->RecoverFEMSolution(X, *b, x); - - // 26. Calculate the relative error of the ROM prediction compared to FOM - ostringstream sol_dofs_name, sol_dofs_name_fom; - if (fom || offline) - { - sol_dofs_name << "sol_dofs_fom." << setfill('0') << setw(6) << myid; - } - if (online) - { - sol_dofs_name << "sol_dofs." << setfill('0') << setw(6) << myid; - sol_dofs_name_fom << "sol_dofs_fom." << setfill('0') << setw(6) << myid; - } - - if (online) - { - // Initialize FOM solution - Vector x_fom(x.Size()); - - ifstream fom_file; - - // Open and load file - fom_file.open(sol_dofs_name_fom.str().c_str()); - - x_fom.Load(fom_file, x_fom.Size()); - - fom_file.close(); - - Vector diff_x(x.Size()); - - subtract(x, x_fom, diff_x); - - // Get norms - double tot_diff_norm = sqrt(InnerProduct(MPI_COMM_WORLD, diff_x, diff_x)); - - double tot_fom_norm = sqrt(InnerProduct(MPI_COMM_WORLD, + if (online) { + // 20. read the reduced basis + assembleTimer.Start(); + CAROM::BasisReader reader(basisName); + spatialbasis = reader.getSpatialBasis(0.0); + numRowRB = spatialbasis->numRows(); + numColumnRB = spatialbasis->numColumns(); + if (myid == 0) printf("spatial basis dimension is %d x %d\n", numRowRB, + numColumnRB); + + // libROM stores the matrix row-wise, so wrapping as a DenseMatrix in MFEM means it is transposed. + DenseMatrix *reducedBasisT = new DenseMatrix(spatialbasis->getData(), + numColumnRB, numRowRB); + + // 21. form inverse ROM operator + CAROM::Matrix invReducedA(numColumnRB, numColumnRB, false); + ComputeCtAB( *A.As(), *spatialbasis, *spatialbasis, + invReducedA); + invReducedA.inverse(); + + CAROM::Vector B_carom(B.GetData(), B.Size(), true, false); + CAROM::Vector X_carom(X.GetData(), X.Size(), true, false); + CAROM::Vector *reducedRHS = spatialbasis->transposeMult(&B_carom); + CAROM::Vector reducedSol(numColumnRB, false); + assembleTimer.Stop(); + + // 22. solve ROM + solveTimer.Start(); + invReducedA.mult(*reducedRHS, reducedSol); + solveTimer.Stop(); + + // 23. reconstruct FOM state + spatialbasis->mult(reducedSol, X_carom); + delete spatialbasis; + delete reducedRHS; + } + + // 24. Recover the parallel grid function corresponding to X. This is the + // local finite element solution on each processor. + a->RecoverFEMSolution(X, *b, x); + + // 26. Calculate the relative error of the ROM prediction compared to FOM + ostringstream sol_dofs_name, sol_dofs_name_fom; + if (fom || offline) + { + sol_dofs_name << "sol_dofs_fom." << setfill('0') << setw(6) << myid; + } + if (online) + { + sol_dofs_name << "sol_dofs." << setfill('0') << setw(6) << myid; + sol_dofs_name_fom << "sol_dofs_fom." << setfill('0') << setw(6) << myid; + } + + if (online) + { + // Initialize FOM solution + Vector x_fom(x.Size()); + + ifstream fom_file; + + // Open and load file + fom_file.open(sol_dofs_name_fom.str().c_str()); + + x_fom.Load(fom_file, x_fom.Size()); + + fom_file.close(); + + Vector diff_x(x.Size()); + + subtract(x, x_fom, diff_x); + + // Get norms + double tot_diff_norm = sqrt(InnerProduct(MPI_COMM_WORLD, diff_x, diff_x)); + + double tot_fom_norm = sqrt(InnerProduct(MPI_COMM_WORLD, x_fom, x_fom)); - if (myid == 0) - { - cout << "Relative error of ROM solution = " - << tot_diff_norm / tot_fom_norm << endl; - } - } - - // 27. Save the refined mesh and the solution in parallel. This output can - // be viewed later using GLVis: "glvis -np -m mesh -g sol". - { - ostringstream mesh_name, sol_name; - mesh_name << "mesh." << setfill('0') << setw(6) << myid; - sol_name << "sol." << setfill('0') << setw(6) << myid; - - ofstream mesh_ofs(mesh_name.str().c_str()); - mesh_ofs.precision(8); - pmesh->Print(mesh_ofs); - - ofstream sol_ofs(sol_name.str().c_str()); - sol_ofs.precision(8); - x.Save(sol_ofs); - } - - // 28. Save data in the VisIt format. - DataCollection *dc = NULL; - if (visit) - { - if (offline) dc = new VisItDataCollection("Example3", pmesh); - else if (fom) dc = new VisItDataCollection("Example3_fom", pmesh); - else if (online) dc = new VisItDataCollection("Example3_rom", pmesh); - dc->SetPrecision(precision); - dc->RegisterField("solution", &x); - dc->Save(); - delete dc; - } - - // 29. Send the solution by socket to a GLVis server. - if (visualization) - { - char vishost[] = "localhost"; - int visport = 19916; - socketstream sol_sock(vishost, visport); - sol_sock << "parallel " << num_procs << " " << myid << "\n"; - sol_sock.precision(8); - sol_sock << "solution\n" << *pmesh << x << flush; - } - - // 30. print timing info - if (myid == 0) - { - if(fom || offline) - { - printf("Elapsed time for assembling FOM: %e second\n", - assembleTimer.RealTime()); - printf("Elapsed time for solving FOM: %e second\n", solveTimer.RealTime()); - } - if(online) - { - printf("Elapsed time for assembling ROM: %e second\n", + if (myid == 0) + { + cout << "Relative error of ROM solution = " + << tot_diff_norm / tot_fom_norm << endl; + } + } + + // 27. Save the refined mesh and the solution in parallel. This output can + // be viewed later using GLVis: "glvis -np -m mesh -g sol". + { + ostringstream mesh_name, sol_name; + mesh_name << "mesh." << setfill('0') << setw(6) << myid; + sol_name << "sol." << setfill('0') << setw(6) << myid; + + ofstream mesh_ofs(mesh_name.str().c_str()); + mesh_ofs.precision(8); + pmesh->Print(mesh_ofs); + + ofstream sol_ofs(sol_name.str().c_str()); + sol_ofs.precision(8); + x.Save(sol_ofs); + } + + // 28. Save data in the VisIt format. + DataCollection *dc = NULL; + if (visit) + { + if (offline) dc = new VisItDataCollection("Example3", pmesh); + else if (fom) dc = new VisItDataCollection("Example3_fom", pmesh); + else if (online) dc = new VisItDataCollection("Example3_rom", pmesh); + dc->SetPrecision(precision); + dc->RegisterField("solution", &x); + dc->Save(); + delete dc; + } + + // 29. Send the solution by socket to a GLVis server. + if (visualization) + { + char vishost[] = "localhost"; + int visport = 19916; + socketstream sol_sock(vishost, visport); + sol_sock << "parallel " << num_procs << " " << myid << "\n"; + sol_sock.precision(8); + sol_sock << "solution\n" << *pmesh << x << flush; + } + + // 30. print timing info + if (myid == 0) + { + if(fom || offline) + { + printf("Elapsed time for assembling FOM: %e second\n", assembleTimer.RealTime()); - printf("Elapsed time for solving ROM: %e second\n", solveTimer.RealTime()); - } - } + printf("Elapsed time for solving FOM: %e second\n", solveTimer.RealTime()); + } + if(online) + { + printf("Elapsed time for assembling ROM: %e second\n", + assembleTimer.RealTime()); + printf("Elapsed time for solving ROM: %e second\n", solveTimer.RealTime()); + } + } - // 31. Free the used memory. - delete a; - delete sigma; - delete muinv; - delete b; - delete fespace; - delete fec; - delete pmesh; + // 31. Free the used memory. + delete a; + delete sigma; + delete muinv; + delete b; + delete fespace; + delete fec; + delete pmesh; - return 0; + return 0; } // 32. define Exact E void E_exact(const Vector &x, Vector &E) { - if (dim == 3) - { - E(0) = sin(kappa * x(1)); - E(1) = sin(kappa * x(2)); - E(2) = sin(kappa * x(0)); - } - else - { - E(0) = sin(kappa * x(1)); - E(1) = sin(kappa * x(0)); - if (x.Size() == 3) { E(2) = 0.0; } - } + if (dim == 3) + { + E(0) = sin(kappa * x(1)); + E(1) = sin(kappa * x(2)); + E(2) = sin(kappa * x(0)); + } + else + { + E(0) = sin(kappa * x(1)); + E(1) = sin(kappa * x(0)); + if (x.Size() == 3) { + E(2) = 0.0; + } + } } // 33. define spatially varying righthand side function void f_exact(const Vector &x, Vector &f) { - if (dim == 3) - { - f(0) = (1. + kappa * kappa) * sin(kappa * x(1)); - f(1) = (1. + kappa * kappa) * sin(kappa * x(2)); - f(2) = (1. + kappa * kappa) * sin(kappa * x(0)); - } - else - { - f(0) = (1. + kappa * kappa) * sin(kappa * x(1)); - f(1) = (1. + kappa * kappa) * sin(kappa * x(0)); - if (x.Size() == 3) { f(2) = 0.0; } - } + if (dim == 3) + { + f(0) = (1. + kappa * kappa) * sin(kappa * x(1)); + f(1) = (1. + kappa * kappa) * sin(kappa * x(2)); + f(2) = (1. + kappa * kappa) * sin(kappa * x(0)); + } + else + { + f(0) = (1. + kappa * kappa) * sin(kappa * x(1)); + f(1) = (1. + kappa * kappa) * sin(kappa * x(0)); + if (x.Size() == 3) { + f(2) = 0.0; + } + } }