From 83074c4c0a27c8553f388273239c2c4c0d3f75d9 Mon Sep 17 00:00:00 2001 From: Tobias Duswald Date: Wed, 25 May 2022 16:19:10 -0700 Subject: [PATCH 1/9] Generalize ex33 to alpha > 0 --- examples/ex33.cpp | 265 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 197 insertions(+), 68 deletions(-) diff --git a/examples/ex33.cpp b/examples/ex33.cpp index 1aec0063b2..cbcb13ad7d 100644 --- a/examples/ex33.cpp +++ b/examples/ex33.cpp @@ -8,29 +8,49 @@ // ex33 -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3 // ex33 -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4 // +// Verification runs: +// ex33 -m ../data/inline-quad.mesh -ver -alpha -o 2 -r 4 +// Note: the analytic solution to this problem is u(x) = sin(pi x) sin(pi y) +// for all alpha. +// // Description: // // In this example we solve the following fractional PDE with MFEM: // -// ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α < 1, +// ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α, // -// To solve this FPDE, we rely on a rational approximation [2] of the normal -// linear operator A^{-α}, where A = - Δ (with associated homogeneous -// boundary conditions). Namely, we first approximate the operator +// To solve this FPDE, we multiply with ( - Δ )^(-N) where the integer +// N is given by floor(α). We obtain // -// A^{-α} ≈ Σ_{i=0}^N c_i (A + d_i I)^{-1}, d_0 = 0, d_i > 0, +// ( - Δ )^(α-N) u = ( - Δ )^(-N) f in Ω, u = 0 on ∂Ω, 0 < α. +// +// We first compute the right hand side by solving the integer order PDE +// +// ( - Δ )^N g = f in Ω, g = 0 on ∂Ω, +// +// The remaining FPDE is then given by +// +// ( - Δ )^(α-N) u = g in Ω, u = 0 on ∂Ω. +// +// We rely on a rational approximation [2] of the normal linear operator +// A^{-α + N}, where A = - Δ (with associated homogeneous boundary conditions) +// and (a-N) in (0,1). We approximate the operator +// +// A^{-α+N} ≈ Σ_{i=0}^M c_i (A + d_i I)^{-1}, d_0 = 0, d_i > 0, // // where I is the L2-identity operator and the coefficients c_i and d_i // are generated offline to a prescribed accuracy in a pre-processing step. // We use the triple-A algorithm [1] to generate the rational approximation -// that this partial fractional expansion derives from. We then solve N+1 +// that this partial fractional expansion derives from. We then solve M+1 // independent integer-order PDEs, // -// A u_i + d_i u_i = c_i f in Ω, u_i = 0 on ∂Ω, i=0,...,N, +// A u_i + d_i u_i = c_i g in Ω, u_i = 0 on ∂Ω, i=0,...,M, // // using MFEM and sum u_i to arrive at an approximate solution of the FPDE // -// u ≈ Σ_{i=0}^N u_i. +// u ≈ Σ_{i=0}^M u_i. +// +// (If alpha is an integer, we stop after the first PDE was solved.) // // References: // @@ -47,6 +67,7 @@ #include "mfem.hpp" #include #include +#include #include "ex33.hpp" @@ -59,8 +80,9 @@ int main(int argc, char *argv[]) const char *mesh_file = "../data/star.mesh"; int order = 1; int num_refs = 3; - bool visualization = true; double alpha = 0.5; + bool visualization = true; + bool verification = false; OptionsParser args(argc, argv); args.AddOption(&mesh_file, "-m", "--mesh", @@ -75,6 +97,9 @@ int main(int argc, char *argv[]) args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", "--no-visualization", "Enable or disable GLVis visualization."); + args.AddOption(&verification, "-ver", "--verification", "-no-ver", + "--no-verification", + "Use sinusoidal function (rhs) for analytic comparison."); args.Parse(); if (!args.Good()) { @@ -84,9 +109,27 @@ int main(int argc, char *argv[]) args.PrintOptions(cout); Array coeffs, poles; + int progress_steps = 1; - // 2. Compute the coefficients that define the integer-order PDEs. - ComputePartialFractionApproximation(alpha,coeffs,poles); + // 2. Compute the rational expansion coefficients that define the + // integer-order PDEs. + int power_of_laplace = floor(alpha); + double exponent_to_approximate = alpha - power_of_laplace; + bool integer_order = false; + // Check if alpha is an integer or not. + if (alpha - power_of_laplace !=0) + { + mfem::out << "Approximating the fractional exponent " + << exponent_to_approximate + << endl; + ComputePartialFractionApproximation(exponent_to_approximate, coeffs, + poles); + } + else + { + integer_order = true; + mfem::out << "Treating integer order PDE." << endl; + } // 3. Read the mesh from the given mesh file. Mesh mesh(mesh_file, 1, 1); @@ -114,7 +157,11 @@ int main(int argc, char *argv[]) } // 7. Define diffusion coefficient, load, and solution GridFunction. - ConstantCoefficient f(1.0); + auto func = [&alpha](const Vector &x) + { + return pow(2*pow(M_PI,2), alpha) * sin(M_PI * x[0]) * sin(M_PI * x[1]); + }; + FunctionCoefficient f(func); ConstantCoefficient one(1.0); GridFunction u(&fespace); u = 0.; @@ -122,71 +169,153 @@ int main(int argc, char *argv[]) // 8. Prepare for visualization. char vishost[] = "localhost"; int visport = 19916; - socketstream xout, uout; - ostringstream oss_x, oss_u; - if (visualization) + + // 9 Compute Right Hand Side + LinearForm rhs(&fespace); + if (verification) { - xout.open(vishost, visport); - xout.precision(8); - uout.open(vishost, visport); - uout.precision(8); + // This statement is only relevant for the verification of the code. It + // uses a different RHS such that an analytic solution is known and easy + // to compare with the numerical one. The FPDE becomes: + // (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2 + // -> u(x,y) = sin(\pi x) sin(\pi y) + rhs.AddDomainIntegrator(new DomainLFIntegrator(f)); + } + else + { + rhs.AddDomainIntegrator(new DomainLFIntegrator(one)); + } + rhs.Assemble(); + + // ------------------------------------------------------------------------ + // 10. Solve the PDE -Δ ^ N g = f, i.e. compute g = (-Δ)^{-1}^N f. + // ------------------------------------------------------------------------ + + if (power_of_laplace > 0) + { + // 10.1 Compute Stiffnes Matrix + BilinearForm k(&fespace); + k.AddDomainIntegrator(new DiffusionIntegrator(one)); + k.Assemble(); + + // 10.2 Compute Mass Matrix + BilinearForm m(&fespace); + m.AddDomainIntegrator(new MassIntegrator(one)); + m.Assemble(); + SparseMatrix mass; + Array empty; + m.FormSystemMatrix(empty, mass); + + // 10.3 from the system of equations + Vector B, X; + GridFunction g(&fespace); + OperatorPtr Op; + k.FormLinearSystem(ess_tdof_list, g, rhs, Op, X, B); + GSSmoother M((SparseMatrix&)(*Op)); + + mfem::out << "\nComputing -Δ ^ -" << power_of_laplace + << " ( f ) " << endl; + for (int i = 0; i < power_of_laplace; i++) + { + // 10.4 Solve the linear system A X = B (N times). + PCG(*Op, M, B, X, 3, 200, 1e-12, 0.0); + // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step + if (visualization && i == power_of_laplace - 1) + { + socketstream fout; + ostringstream oss_f; + fout.open(vishost, visport); + fout.precision(8); + k.RecoverFEMSolution(X, rhs, g); + oss_f.str(""); oss_f.clear(); + oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ " + << power_of_laplace + << " g = f"; + fout << "solution\n" << mesh << g + << "window_title '" << oss_f.str() << "'" << flush; + } + mass.Mult(X, B); + X.SetSubVectorComplement(ess_tdof_list,0.0); + } + // 10.6 Extract solution for the next step. + rhs = B; } - for (int i = 0; i < coeffs.Size(); i++) + // ------------------------------------------------------------------------ + // 11. Solve the fractional PDE by solving M integer order PDEs and adding + // up the solutions. + // ------------------------------------------------------------------------ + if (!integer_order) { - // 9. Set up the linear form b(.) for integer-order PDE solve. - LinearForm b(&fespace); - ProductCoefficient cf(coeffs[i], f); - b.AddDomainIntegrator(new DomainLFIntegrator(cf)); - b.Assemble(); - - // 10. Define GridFunction for integer-order PDE solve. - GridFunction x(&fespace); - x = 0.0; - - // 11. Set up the bilinear form a(.,.) for integer-order PDE solve. - BilinearForm a(&fespace); - a.AddDomainIntegrator(new DiffusionIntegrator(one)); - ConstantCoefficient c2(-poles[i]); - a.AddDomainIntegrator(new MassIntegrator(c2)); - a.Assemble(); - - // 12. Assemble the bilinear form and the corresponding linear system. - OperatorPtr A; - Vector B, X; - a.FormLinearSystem(ess_tdof_list, x, b, A, X, B); - - // 13. Solve the linear system A X = B. - GSSmoother M((SparseMatrix&)(*A)); - - mfem::out << "\nSolving PDE -Δ u + " << -poles[i] - << " u = " << coeffs[i] << " f " << endl; - PCG(*A, M, B, X, 3, 200, 1e-12, 0.0); - - // 14. Recover the solution as a finite element grid function. - a.RecoverFEMSolution(X, b, x); - - // 15. Accumulate integer-order PDE solutions. - u+=x; - - // 16. Send the solutions by socket to a GLVis server. + // Setup visualization. + socketstream xout, uout; + ostringstream oss_x, oss_u; if (visualization) { - oss_x.str(""); oss_x.clear(); - oss_x << "Solution of PDE -Δ u + " << -poles[i] - << " u = " << coeffs[i] << " f"; - xout << "solution\n" << mesh << x - << "window_title '" << oss_x.str() << "'" << flush; + xout.open(vishost, visport); + xout.precision(8); + uout.open(vishost, visport); + uout.precision(8); + } + // Iterate over all expansion coefficient that contribute to the + // solution. + for (int i = 0; i < coeffs.Size(); i++) + { + mfem::out << "\nSolving PDE -Δ u + " << -poles[i] + << " u = " << coeffs[i] << " f " << endl; - oss_u.str(""); oss_u.clear(); - oss_u << "Solution of fractional PDE -Δ^" << alpha - << " u = f"; - uout << "solution\n" << mesh << u - << "window_title '" << oss_u.str() << "'" << flush; + // 11.1 Set up the linear form b(.) for integer-order PDE solve. + // (c_i * g) + Vector b (rhs); + b *= coeffs[i]; + + // 11.2 Define GridFunction for integer-order PDE solve. + GridFunction x(&fespace); + x = 0.0; + + // 11.3 Set up the bilinear form a(.,.) for integer-order PDE solve. + BilinearForm a(&fespace); + a.AddDomainIntegrator(new DiffusionIntegrator(one)); + ConstantCoefficient c2(-poles[i]); + a.AddDomainIntegrator(new MassIntegrator(c2)); + a.Assemble(); + + // 11.4 Assemble the bilinear form and the corresponding linear system. + OperatorPtr A; + Vector B, X; + a.FormLinearSystem(ess_tdof_list, x, b, A, X, B); + + // 11.5 Solve the linear system A X = B. + GSSmoother M((SparseMatrix&)(*A)); + + PCG(*A, M, B, X, 3, 200, 1e-12, 0.0); + + // 11.6 Recover the solution as a finite element grid function. + a.RecoverFEMSolution(X, b, x); + + // 11.7 Accumulate integer-order PDE solutions. + u+=x; + + // 11.8 Send the solutions by socket to a GLVis server. + if (visualization) + { + oss_x.str(""); oss_x.clear(); + oss_x << "Step " << progress_steps + << ": Solution of PDE -Δ u + " << -poles[i] + << " u = " << coeffs[i] << " f"; + xout << "solution\n" << mesh << x + << "window_title '" << oss_x.str() << "'" << flush; + + oss_u.str(""); oss_u.clear(); + oss_u << "Step " << progress_steps + 1 + << ": Solution of fractional PDE -Δ^" << alpha + << " u = f"; + uout << "solution\n" << mesh << u + << "window_title '" << oss_u.str() << "'" << flush; + } } } - - // 17. Free the used memory. + // 12. Free the used memory. delete fec; return 0; } From 043b338fee284dfe44b6b0dd342487a39050b2e9 Mon Sep 17 00:00:00 2001 From: Tobias Duswald Date: Wed, 25 May 2022 16:26:32 -0700 Subject: [PATCH 2/9] Add explanation for computation of coefficients --- examples/ex33.hpp | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/examples/ex33.hpp b/examples/ex33.hpp index 3da5a46895..aeb4c7594c 100644 --- a/examples/ex33.hpp +++ b/examples/ex33.hpp @@ -249,6 +249,13 @@ void PartialFractionExpansion(double scale, Array & poles, coeffs.SetSize(psize); coeffs = scale; + // Note: C p(z)/q(z) = Σ_i c_i / (z - p_i) results in an system of equations + // where the N unknowns are the coefficients c_i. After multiplying the + // system with q(z), the coefficients c_i can be computed analytically by + // choosing N values for z. Choosing z_j = = p_j diagonalizes the system and + // one can obtain an analytic form for the c_i coefficients. The result is + // implemented in the code block below. + for (int i=0; i Date: Wed, 25 May 2022 16:28:12 -0700 Subject: [PATCH 3/9] Astyle Format --- examples/ex33.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/examples/ex33.hpp b/examples/ex33.hpp index aeb4c7594c..b88744cb95 100644 --- a/examples/ex33.hpp +++ b/examples/ex33.hpp @@ -250,10 +250,10 @@ void PartialFractionExpansion(double scale, Array & poles, coeffs = scale; // Note: C p(z)/q(z) = Σ_i c_i / (z - p_i) results in an system of equations - // where the N unknowns are the coefficients c_i. After multiplying the - // system with q(z), the coefficients c_i can be computed analytically by - // choosing N values for z. Choosing z_j = = p_j diagonalizes the system and - // one can obtain an analytic form for the c_i coefficients. The result is + // where the N unknowns are the coefficients c_i. After multiplying the + // system with q(z), the coefficients c_i can be computed analytically by + // choosing N values for z. Choosing z_j = = p_j diagonalizes the system and + // one can obtain an analytic form for the c_i coefficients. The result is // implemented in the code block below. for (int i=0; i Date: Thu, 26 May 2022 17:35:40 -0700 Subject: [PATCH 4/9] Update ex33p to new logic --- examples/ex33p.cpp | 399 +++++++++++++++++++++++++++++---------------- 1 file changed, 258 insertions(+), 141 deletions(-) diff --git a/examples/ex33p.cpp b/examples/ex33p.cpp index f4723efc5e..3231e9cd0e 100644 --- a/examples/ex33p.cpp +++ b/examples/ex33p.cpp @@ -8,29 +8,49 @@ // mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3 // mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4 // +// Verification runs: +// mpirun -np 4ex33 -m ../data/inline-quad.mesh -ver -alpha -o 2 -r 4 +// Note: the analytic solution to this problem is u(x) = sin(pi x) sin(pi y) +// for all alpha. +// // Description: // // In this example we solve the following fractional PDE with MFEM: // -// ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α < 1, +// ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α, // -// To solve this FPDE, we rely on a rational approximation [2] of the normal -// linear operator A^{-α}, where A = - Δ (with associated homogeneous -// boundary conditions). Namely, we first approximate the operator +// To solve this FPDE, we multiply with ( - Δ )^(-N) where the integer +// N is given by floor(α). We obtain // -// A^{-α} ≈ Σ_{i=0}^N c_i (A + d_i I)^{-1}, d_0 = 0, d_i > 0, +// ( - Δ )^(α-N) u = ( - Δ )^(-N) f in Ω, u = 0 on ∂Ω, 0 < α. +// +// We first compute the right hand side by solving the integer order PDE +// +// ( - Δ )^N g = f in Ω, g = 0 on ∂Ω, +// +// The remaining FPDE is then given by +// +// ( - Δ )^(α-N) u = g in Ω, u = 0 on ∂Ω. +// +// We rely on a rational approximation [2] of the normal linear operator +// A^{-α + N}, where A = - Δ (with associated homogeneous boundary conditions) +// and (a-N) in (0,1). We approximate the operator +// +// A^{-α+N} ≈ Σ_{i=0}^M c_i (A + d_i I)^{-1}, d_0 = 0, d_i > 0, // // where I is the L2-identity operator and the coefficients c_i and d_i // are generated offline to a prescribed accuracy in a pre-processing step. // We use the triple-A algorithm [1] to generate the rational approximation -// that this partial fractional expansion derives from. We then solve N+1 +// that this partial fractional expansion derives from. We then solve M+1 // independent integer-order PDEs, // -// A u_i + d_i u_i = c_i f in Ω, u_i = 0 on ∂Ω, i=0,...,N, +// A u_i + d_i u_i = c_i g in Ω, u_i = 0 on ∂Ω, i=0,...,M, // // using MFEM and sum u_i to arrive at an approximate solution of the FPDE // -// u ≈ Σ_{i=0}^N u_i. +// u ≈ Σ_{i=0}^M u_i. +// +// (If alpha is an integer, we stop after the first PDE was solved.) // // References: // @@ -47,6 +67,8 @@ #include "mfem.hpp" #include #include +#include +#include #include "ex33.hpp" @@ -60,14 +82,18 @@ int main(int argc, char *argv[]) int num_procs = Mpi::WorldSize(); int myid = Mpi::WorldRank(); Hypre::Init(); + if (Mpi::Root()) + { + mfem::out << "\nTotal number of MPI ranks = " << num_procs << endl; + } // 1. Parse command-line options. const char *mesh_file = "../data/star.mesh"; int order = 1; int num_refs = 3; - bool visualization = true; - bool visualize_x = false; double alpha = 0.5; + bool visualization = true; + bool verification = false; OptionsParser args(argc, argv); args.AddOption(&mesh_file, "-m", "--mesh", @@ -79,12 +105,12 @@ int main(int argc, char *argv[]) "Number of uniform refinements"); args.AddOption(&alpha, "-alpha", "--alpha", "Fractional exponent"); - args.AddOption(&visualize_x, "-vis_x", "--visualize_x", "-no-vis_x", - "--no-visualization_x", - "Enable or disable GLVis visualization of each integer-order PDE solution."); args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", "--no-visualization", - "Enable or disable GLVis visualization of the fractional PDE solution."); + "Enable or disable GLVis visualization."); + args.AddOption(&verification, "-ver", "--verification", "-no-ver", + "--no-verification", + "Use sinusoidal function (f) for analytic comparison."); args.Parse(); if (!args.Good()) { @@ -97,61 +123,47 @@ int main(int argc, char *argv[]) } Array coeffs, poles; + int progress_steps = 1; - // 2. Compute the coefficients that define the integer-order PDEs. - ComputePartialFractionApproximation(alpha,coeffs,poles); - - int num_par_solves; - int max_par_solves = max(1,num_procs/2); - for (num_par_solves=max_par_solves; num_par_solves>0; num_par_solves--) + // 2. Compute the rational expansion coefficients that define the + // integer-order PDEs. + int power_of_laplace = floor(alpha); + double exponent_to_approximate = alpha - power_of_laplace; + bool integer_order = false; + // Check if alpha is an integer or not. + if (alpha - power_of_laplace !=0) { - if (num_procs%num_par_solves==0 && num_par_solves ess_tdof_list; if (pmesh.bdr_attributes.Size()) { @@ -169,120 +181,225 @@ int main(int argc, char *argv[]) fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list); } - // 8. Define diffusion coefficient, load, and solution GridFunction. - ConstantCoefficient f(1.0); + // 7. Define diffusion coefficient, load, and solution GridFunction. + auto func = [&alpha](const Vector &x) + { + return pow(2*pow(M_PI,2), alpha) * sin(M_PI * x[0]) * sin(M_PI * x[1]); + }; + FunctionCoefficient f(func); ConstantCoefficient one(1.0); ParGridFunction u(&fespace); ParGridFunction x(&fespace); + ParGridFunction g(&fespace); u = 0.0; + x = 0.0; + g = 0.0; + + // 8. Prepare for visualization. + char vishost[] = "localhost"; + int visport = 19916; // 9. Set up the linear form b(.) for integer-order PDE solves. ParLinearForm b(&fespace); - b.AddDomainIntegrator(new DomainLFIntegrator(f)); + if (verification) + { + // This statement is only relevant for the verification of the code. It + // uses a different f such that an analytic solution is known and easy + // to compare with the numerical one. The FPDE becomes: + // (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2 + // -> u(x,y) = sin(\pi x) sin(\pi y) + b.AddDomainIntegrator(new DomainLFIntegrator(f)); + } + else + { + b.AddDomainIntegrator(new DomainLFIntegrator(one)); + } b.Assemble(); - int my_coeff_size = max(coeffs.Size()/col_size,1); - int ibeg = col_rank*my_coeff_size; - if (ibeg + 2*my_coeff_size > coeffs.Size()) + // ------------------------------------------------------------------------ + // 10. Solve the PDE -Δ ^ N g = f, i.e. compute g = (-Δ)^{-1}^N f. + // ------------------------------------------------------------------------ + + if (power_of_laplace > 0) { - my_coeff_size = coeffs.Size()-col_rank*my_coeff_size; - } - else if (ibeg > coeffs.Size() - 1) - { - my_coeff_size = 0; - } + // 10.1 Compute Stiffnes Matrix + ParBilinearForm k(&fespace); + k.AddDomainIntegrator(new DiffusionIntegrator(one)); + k.Assemble(); - int iend = ibeg+my_coeff_size; + // 10.2 Compute Mass Matrix + ParBilinearForm m(&fespace); + m.AddDomainIntegrator(new MassIntegrator(one)); + m.Assemble(); + HypreParMatrix mass; + Array empty; + m.FormSystemMatrix(empty, mass); - - for (int i = ibeg; i < iend; i++) - { - // 10. Reset GridFunction for integer-order PDE solve. - x = 0.0; - - // 11. Set up the bilinear form a(.,.) for integer-order PDE solve. - ParBilinearForm a(&fespace); - a.AddDomainIntegrator(new DiffusionIntegrator(one)); - ConstantCoefficient d_i(-poles[i]); - a.AddDomainIntegrator(new MassIntegrator(d_i)); - a.Assemble(); - - // 12. Assemble the bilinear form and the corresponding linear system. - OperatorPtr A; + // 10.3 from the system of equations Vector B, X; - a.FormLinearSystem(ess_tdof_list, x, b, A, X, B); - - // 13. Solve the linear system A X = B. - HypreBoomerAMG * prec = new HypreBoomerAMG; - prec->SetPrintLevel(-1); - - int print_level = (col_rank==0) ? 3 : 0; - if (Mpi::Root()) - { - mfem::out << "\nMPI rank " << myid - << ": Solving PDE -Δ u + " << -poles[i] - << " u = " << coeffs[i] << " f " << endl; - } - CGSolver cg(row_comm); + OperatorPtr Op; + k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B); + HypreBoomerAMG prec; + prec.SetPrintLevel(-1); + int print_level = (myid==0) ? 3 : 0; + CGSolver cg(MPI_COMM_WORLD); cg.SetRelTol(1e-12); cg.SetMaxIter(2000); cg.SetPrintLevel(print_level); - cg.SetPreconditioner(*prec); - cg.SetOperator(*A); - cg.Mult(B, X); - delete prec; + cg.SetPreconditioner(prec); + cg.SetOperator(*Op); - // 14. Recover the solution as a finite element grid function. - a.RecoverFEMSolution(X, b, x); - - // 15. Accumulate integer-order PDE solutions. - x *= coeffs[i]; - u += x; - - // 16. Send integer-order PDE solutions to a GLVis server. - if (visualize_x) + if (Mpi::Root()) { - if (col_rank > 0 && i < iend-1) + mfem::out << "\nComputing -Δ ^ -" << power_of_laplace + << " ( f ) " << endl; + } + for (int i = 0; i < power_of_laplace; i++) + { + // 10.4 Solve the linear system A X = B (N times). + cg.Mult(B, X); + // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step + if (i == power_of_laplace - 1) { - MPI_Status status; - MPI_Recv(nullptr,0,MPI_INT, col_rank-1,0,col_comm,&status); + // Needed for visualization and solution verification. + k.RecoverFEMSolution(X, b, g); + if (integer_order && verification) + { + // For an integer order PDE, g is also our solution u. + u+=g; + } + if (visualization) + { + socketstream fout; + ostringstream oss_f; + fout.open(vishost, visport); + fout.precision(8); + oss_f.str(""); oss_f.clear(); + oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ " + << power_of_laplace + << " g = f"; + fout << "parallel " << num_procs << " " << myid << "\n" + << "solution\n" << pmesh << g + << "window_title '" << oss_f.str() << "'" << flush; + } } - char vishost[] = "localhost"; - int visport = 19916; - socketstream xout(vishost, visport); + + // 10.6 Prepare for next iteration (primal / dual space) + mass.Mult(X, B); + X.SetSubVectorComplement(ess_tdof_list,0.0); + } + + // 10.7 Extract solution for the next step. The b now corresponds to the + // function g in the PDE. + const SparseMatrix* rm = fespace.GetRestrictionMatrix(); + rm->MultTranspose(B, b); + } + + // ------------------------------------------------------------------------ + // 11. Solve the fractional PDE by solving M integer order PDEs and adding + // up the solutions. + // ------------------------------------------------------------------------ + if (!integer_order) + { + // Setup visualization. + socketstream xout, uout; + ostringstream oss_x, oss_u; + if (visualization) + { + xout.open(vishost, visport); xout.precision(8); - ostringstream oss; - oss << "Solution of PDE -Δ u + " << -poles[i] - << " u = " << coeffs[i] << " f" ; - xout << "parallel " << row_size << " " << row_rank << "\n"; - xout << "solution\n" << pmesh << x - << "window_title '" << oss.str() << "'" << flush; - if (col_rank < col_size-1) + uout.open(vishost, visport); + uout.precision(8); + } + // Iterate over all expansion coefficient that contribute to the + // solution. + for (int i = 0; i < coeffs.Size(); i++) + { + if (Mpi::Root()) { - MPI_Send(nullptr,0,MPI_INT,col_rank+1,0,col_comm); + mfem::out << "\nSolving PDE -Δ u + " << -poles[i] + << " u = " << coeffs[i] << " g " << endl; + } + + // 11.1 Reset GridFunction for integer-order PDE solve. + x = 0.0; + + // 11.2 Set up the bilinear form a(.,.) for integer-order PDE solve. + ParBilinearForm a(&fespace); + a.AddDomainIntegrator(new DiffusionIntegrator(one)); + ConstantCoefficient d_i(-poles[i]); + a.AddDomainIntegrator(new MassIntegrator(d_i)); + a.Assemble(); + + // 11.3 Assemble the bilinear form and the corresponding linear system. + OperatorPtr A; + Vector B, X; + a.FormLinearSystem(ess_tdof_list, x, b, A, X, B); + + // 11.4 Solve the linear system A X = B. + HypreBoomerAMG prec; + prec.SetPrintLevel(-1); + + int print_level = (myid==0) ? 3 : 0; + CGSolver cg(MPI_COMM_WORLD); + cg.SetRelTol(1e-12); + cg.SetMaxIter(2000); + cg.SetPrintLevel(print_level); + cg.SetPreconditioner(prec); + cg.SetOperator(*A); + cg.Mult(B, X); + + // 11.5 Recover the solution as a finite element grid function. + a.RecoverFEMSolution(X, b, x); + + // 11.6 Accumulate integer-order PDE solutions. + x *= coeffs[i]; + u += x; + + // 11.7 Send fractional PDE solution to a GLVis server. + if (visualization) + { + oss_x.str(""); oss_x.clear(); + oss_x << "Step " << progress_steps + << ": Solution of PDE -Δ u + " << -poles[i] + << " u = " << coeffs[i] << " g"; + xout << "parallel " << num_procs << " " << myid << "\n" + << "solution\n" << pmesh << x + << "window_title '" << oss_x.str() << "'" << flush; + + oss_u.str(""); oss_u.clear(); + oss_u << "Step " << progress_steps + 1 + << ": Solution of fractional PDE -Δ^" << alpha - floor(alpha) + << " u = g"; + uout << "parallel " << num_procs << " " << myid << "\n" + << "solution\n" << pmesh << u + << "window_title '" << oss_u.str() << "'" + << flush; } } } - // 17. Accumulate for the fractional PDE solution - MPI_Allreduce(MPI_IN_PLACE, u.GetData(), u.Size(), - MPI_DOUBLE, MPI_SUM,col_comm); - - // 18. Send fractional PDE solution to a GLVis server. - if (visualization) + // ------------------------------------------------------------------------ + // 12. (optional) Verify the solution. + // ------------------------------------------------------------------------ + if (verification) { - if (col_rank == 0) + auto solution = [] (const Vector &x) { - char vishost[] = "localhost"; - int visport = 19916; - socketstream uout(vishost, visport); - uout.precision(8); - ostringstream oss; - oss << "Solution of fractional PDE -Δ^" << alpha - << " u = f" ; - uout << "parallel " << row_size << " " << row_rank << "\n"; - uout << "solution\n" << pmesh << u - << "window_title '" << oss.str() << "'" << flush; + return sin(M_PI * x[0]) * sin(M_PI * x[1]); + }; + FunctionCoefficient sol(solution); + double l2_error = u.ComputeL2Error(sol); + + if (Mpi::Root()) + { + mfem::out << "\n" << string(80,'=') + << "\n\nSolution Verification\n\n" + << "Analytic solution : sin(pi x) sin(pi y)\n" + << "Expected mesh : inline_quad.mesh\n" + << "Your mesh : " << mesh_file << "\n" + << "L2 error : " << l2_error << "\n\n" + << string(80,'=') << endl; } } From 863638bb476f675cc033647f0cc25a8b158c7c88 Mon Sep 17 00:00:00 2001 From: Tobias Duswald Date: Thu, 26 May 2022 17:36:05 -0700 Subject: [PATCH 5/9] Make serial and parallel example identical --- examples/ex33.cpp | 131 +++++++++++++++++++++++++++++----------------- 1 file changed, 83 insertions(+), 48 deletions(-) diff --git a/examples/ex33.cpp b/examples/ex33.cpp index cbcb13ad7d..e79dff17c0 100644 --- a/examples/ex33.cpp +++ b/examples/ex33.cpp @@ -68,6 +68,7 @@ #include #include #include +#include #include "ex33.hpp" @@ -99,7 +100,7 @@ int main(int argc, char *argv[]) "Enable or disable GLVis visualization."); args.AddOption(&verification, "-ver", "--verification", "-no-ver", "--no-verification", - "Use sinusoidal function (rhs) for analytic comparison."); + "Use sinusoidal function (f) for analytic comparison."); args.Parse(); if (!args.Good()) { @@ -142,8 +143,8 @@ int main(int argc, char *argv[]) } // 5. Define a finite element space on the mesh. - FiniteElementCollection *fec = new H1_FECollection(order, dim); - FiniteElementSpace fespace(&mesh, fec); + H1_FECollection fec(order, dim); + FiniteElementSpace fespace(&mesh, &fec); cout << "Number of finite element unknowns: " << fespace.GetTrueVSize() << endl; @@ -164,28 +165,32 @@ int main(int argc, char *argv[]) FunctionCoefficient f(func); ConstantCoefficient one(1.0); GridFunction u(&fespace); - u = 0.; + GridFunction x(&fespace); + GridFunction g(&fespace); + u = 0.0; + x = 0.0; + g = 0.0; // 8. Prepare for visualization. char vishost[] = "localhost"; int visport = 19916; - // 9 Compute Right Hand Side - LinearForm rhs(&fespace); + // 9. Set up the linear form b(.) for integer-order PDE solves. + LinearForm b(&fespace); if (verification) { // This statement is only relevant for the verification of the code. It - // uses a different RHS such that an analytic solution is known and easy + // uses a different f such that an analytic solution is known and easy // to compare with the numerical one. The FPDE becomes: // (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2 // -> u(x,y) = sin(\pi x) sin(\pi y) - rhs.AddDomainIntegrator(new DomainLFIntegrator(f)); + b.AddDomainIntegrator(new DomainLFIntegrator(f)); } else { - rhs.AddDomainIntegrator(new DomainLFIntegrator(one)); + b.AddDomainIntegrator(new DomainLFIntegrator(one)); } - rhs.Assemble(); + b.Assemble(); // ------------------------------------------------------------------------ // 10. Solve the PDE -Δ ^ N g = f, i.e. compute g = (-Δ)^{-1}^N f. @@ -208,9 +213,8 @@ int main(int argc, char *argv[]) // 10.3 from the system of equations Vector B, X; - GridFunction g(&fespace); OperatorPtr Op; - k.FormLinearSystem(ess_tdof_list, g, rhs, Op, X, B); + k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B); GSSmoother M((SparseMatrix&)(*Op)); mfem::out << "\nComputing -Δ ^ -" << power_of_laplace @@ -219,26 +223,40 @@ int main(int argc, char *argv[]) { // 10.4 Solve the linear system A X = B (N times). PCG(*Op, M, B, X, 3, 200, 1e-12, 0.0); + // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step - if (visualization && i == power_of_laplace - 1) + if (i == power_of_laplace - 1) { - socketstream fout; - ostringstream oss_f; - fout.open(vishost, visport); - fout.precision(8); - k.RecoverFEMSolution(X, rhs, g); - oss_f.str(""); oss_f.clear(); - oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ " - << power_of_laplace - << " g = f"; - fout << "solution\n" << mesh << g - << "window_title '" << oss_f.str() << "'" << flush; + // Needed for visualization and solution verification. + k.RecoverFEMSolution(X, b, g); + if (integer_order && verification) + { + // For an integer order PDE, g is also our solution u. + u+=g; + } + if (visualization) + { + socketstream fout; + ostringstream oss_f; + fout.open(vishost, visport); + fout.precision(8); + oss_f.str(""); oss_f.clear(); + oss_f << "Step " << progress_steps++ << ": Solution of PDE -Δ ^ " + << power_of_laplace + << " g = f"; + fout << "solution\n" << mesh << g + << "window_title '" << oss_f.str() << "'" << flush; + } } + + // 10.6 Prepare for next iteration (primal / dual space) mass.Mult(X, B); X.SetSubVectorComplement(ess_tdof_list,0.0); } - // 10.6 Extract solution for the next step. - rhs = B; + + // 10.7 Extract solution for the next step. The b now corresponds to the + // function g in the PDE. + b = B; } // ------------------------------------------------------------------------ @@ -262,60 +280,77 @@ int main(int argc, char *argv[]) for (int i = 0; i < coeffs.Size(); i++) { mfem::out << "\nSolving PDE -Δ u + " << -poles[i] - << " u = " << coeffs[i] << " f " << endl; + << " u = " << coeffs[i] << " g " << endl; - // 11.1 Set up the linear form b(.) for integer-order PDE solve. - // (c_i * g) - Vector b (rhs); - b *= coeffs[i]; - // 11.2 Define GridFunction for integer-order PDE solve. - GridFunction x(&fespace); + // 11.1 Reset GridFunction for integer-order PDE solve. x = 0.0; - // 11.3 Set up the bilinear form a(.,.) for integer-order PDE solve. + // 11.2 Set up the bilinear form a(.,.) for integer-order PDE solve. BilinearForm a(&fespace); a.AddDomainIntegrator(new DiffusionIntegrator(one)); - ConstantCoefficient c2(-poles[i]); - a.AddDomainIntegrator(new MassIntegrator(c2)); + ConstantCoefficient d_i(-poles[i]); + a.AddDomainIntegrator(new MassIntegrator(d_i)); a.Assemble(); - // 11.4 Assemble the bilinear form and the corresponding linear system. + // 11.3 Assemble the bilinear form and the corresponding linear system. OperatorPtr A; Vector B, X; a.FormLinearSystem(ess_tdof_list, x, b, A, X, B); - // 11.5 Solve the linear system A X = B. + // 11.4 Solve the linear system A X = B. GSSmoother M((SparseMatrix&)(*A)); PCG(*A, M, B, X, 3, 200, 1e-12, 0.0); - // 11.6 Recover the solution as a finite element grid function. + // 11.5 Recover the solution as a finite element grid function. a.RecoverFEMSolution(X, b, x); - // 11.7 Accumulate integer-order PDE solutions. - u+=x; + // 11.6 Accumulate integer-order PDE solutions. + x *= coeffs[i]; + u += x; - // 11.8 Send the solutions by socket to a GLVis server. + // 11.7 Send fractional PDE solution to a GLVis server. if (visualization) { oss_x.str(""); oss_x.clear(); oss_x << "Step " << progress_steps << ": Solution of PDE -Δ u + " << -poles[i] - << " u = " << coeffs[i] << " f"; + << " u = " << coeffs[i] << " g"; xout << "solution\n" << mesh << x << "window_title '" << oss_x.str() << "'" << flush; oss_u.str(""); oss_u.clear(); oss_u << "Step " << progress_steps + 1 - << ": Solution of fractional PDE -Δ^" << alpha - << " u = f"; + << ": Solution of fractional PDE -Δ^" << alpha - floor(alpha) + << " u = g"; uout << "solution\n" << mesh << u - << "window_title '" << oss_u.str() << "'" << flush; + << "window_title '" << oss_u.str() << "'" + << flush; } } } - // 12. Free the used memory. - delete fec; + + // ------------------------------------------------------------------------ + // 12. (optional) Verify the solution. + // ------------------------------------------------------------------------ + if (verification) + { + auto solution = [] (const Vector &x) + { + return sin(M_PI * x[0]) * sin(M_PI * x[1]); + }; + FunctionCoefficient sol(solution); + double l2_error = u.ComputeL2Error(sol); + + mfem::out << "\n" << string(80,'=') + << "\n\nSolution Verification\n\n" + << "Analytic solution : sin(pi x) sin(pi y)\n" + << "Expected mesh : inline_quad.mesh\n" + << "Your mesh : " << mesh_file << "\n" + << "L2 error : " << l2_error << "\n\n" + << string(80,'=') << endl; + } + return 0; } From e91c582d032ea3c82c8f9d24f5d6e364bc2928e2 Mon Sep 17 00:00:00 2001 From: Tobias Duswald Date: Fri, 27 May 2022 18:04:57 -0700 Subject: [PATCH 6/9] Apply changes from review Co-authored-by: brendankeith Co-authored-by: psocratis --- examples/ex33.cpp | 34 ++++++++++++++++++++++------------ examples/ex33.hpp | 12 ++++++++---- examples/ex33p.cpp | 44 ++++++++++++++++++++++++-------------------- 3 files changed, 54 insertions(+), 36 deletions(-) diff --git a/examples/ex33.cpp b/examples/ex33.cpp index e79dff17c0..13bca1a5e2 100644 --- a/examples/ex33.cpp +++ b/examples/ex33.cpp @@ -3,15 +3,21 @@ // Compile with: make ex33 // // Sample runs: ex33 -m ../data/square-disc.mesh -alpha 0.33 -o 2 +// ex33 -m ../data/square-disc.mesh -alpha 4.5 -o 3 +// ex33 -m ../data/star.mesh -alpha 1.4 -o 3 // ex33 -m ../data/star.mesh -alpha 0.99 -o 3 // ex33 -m ../data/inline-quad.mesh -alpha 0.5 -o 3 // ex33 -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3 +// ex33 -m ../data/disc-nurbs.mesh -alpha 2.4 -o 3 -r 4 // ex33 -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4 +// ex33 -m ../data/l-shape.mesh -alpha 1.7 -o 3 -r 5 // // Verification runs: -// ex33 -m ../data/inline-quad.mesh -ver -alpha -o 2 -r 4 -// Note: the analytic solution to this problem is u(x) = sin(pi x) sin(pi y) -// for all alpha. +// ex33 -m ../data/inline-quad.mesh -ver -alpha 0.3 -o 2 -r 4 +// ex33 -m ../data/inline-quad.mesh -ver -alpha 1.4 -o 3 -r 5 +// ex33 -m ../data/inline-quad.mesh -ver -alpha 5.7 -o 2 -r 7 -no-vis +// Note: the analytic solution to this problem is u(x) = sin(pi x) sin(pi y) +// for all alpha. // // Description: // @@ -19,14 +25,14 @@ // // ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α, // -// To solve this FPDE, we multiply with ( - Δ )^(-N) where the integer -// N is given by floor(α). We obtain +// To solve this FPDE, we apply the operator ( - Δ )^(-N), where the integer +// N is given by floor(α). By doing so, we obtain // // ( - Δ )^(α-N) u = ( - Δ )^(-N) f in Ω, u = 0 on ∂Ω, 0 < α. // // We first compute the right hand side by solving the integer order PDE // -// ( - Δ )^N g = f in Ω, g = 0 on ∂Ω, +// ( - Δ )^N g = f in Ω, g = ( - Δ )^k g = 0 on ∂Ω, k = 1,..,N-1 // // The remaining FPDE is then given by // @@ -114,17 +120,21 @@ int main(int argc, char *argv[]) // 2. Compute the rational expansion coefficients that define the // integer-order PDEs. - int power_of_laplace = floor(alpha); + const int power_of_laplace = floor(alpha); double exponent_to_approximate = alpha - power_of_laplace; bool integer_order = false; // Check if alpha is an integer or not. - if (alpha - power_of_laplace !=0) + if (abs(exponent_to_approximate) > 1e-12) { mfem::out << "Approximating the fractional exponent " << exponent_to_approximate << endl; ComputePartialFractionApproximation(exponent_to_approximate, coeffs, poles); + + // If the example is build without LAPACK, the exponent_to_approximate + // might be modified by the function call above. + alpha = exponent_to_approximate + power_of_laplace; } else { @@ -193,7 +203,7 @@ int main(int argc, char *argv[]) b.Assemble(); // ------------------------------------------------------------------------ - // 10. Solve the PDE -Δ ^ N g = f, i.e. compute g = (-Δ)^{-1}^N f. + // 10. Solve the PDE (-Δ)^N g = f, i.e. compute g = (-Δ)^{-1}^N f. // ------------------------------------------------------------------------ if (power_of_laplace > 0) @@ -217,7 +227,7 @@ int main(int argc, char *argv[]) k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B); GSSmoother M((SparseMatrix&)(*Op)); - mfem::out << "\nComputing -Δ ^ -" << power_of_laplace + mfem::out << "\nComputing (-Δ) ^ -" << power_of_laplace << " ( f ) " << endl; for (int i = 0; i < power_of_laplace; i++) { @@ -322,8 +332,8 @@ int main(int argc, char *argv[]) oss_u.str(""); oss_u.clear(); oss_u << "Step " << progress_steps + 1 - << ": Solution of fractional PDE -Δ^" << alpha - floor(alpha) - << " u = g"; + << ": Solution of fractional PDE (-Δ)^" << alpha + << " u = f"; uout << "solution\n" << mesh << u << "window_title '" << oss_u.str() << "'" << flush; diff --git a/examples/ex33.hpp b/examples/ex33.hpp index b88744cb95..9313d67ca9 100644 --- a/examples/ex33.hpp +++ b/examples/ex33.hpp @@ -32,6 +32,7 @@ #include "mfem.hpp" #include #include +#include using namespace std; using namespace mfem; @@ -312,9 +313,12 @@ void ComputePartialFractionApproximation(double & alpha, if (print_warning) { mfem::out - << "\nMFEM is compiled without LAPACK.\nUsing precomputed values for PartialFractionApproximation. \n" - << "Only alpha = 0.33, 0.5, and 0.99 are available.\nThe default is alpha = 0.5." - << std::endl; + << "\n" << string(80, '=') + << "\nMFEM is compiled without LAPACK." + << "\nUsing precomputed values for PartialFractionApproximation." + << "\nOnly alpha = 0.33, 0.5, and 0.99 are available." + << "\nThe default is alpha = 0.5.\n" << string(80, '=') << "\n" + << endl; } const double eps = std::numeric_limits::epsilon(); @@ -358,7 +362,7 @@ void ComputePartialFractionApproximation(double & alpha, if (print_warning) { - mfem::out << "Using precomputed values for alpha = " + mfem::out << "=> Using precomputed values for alpha = " << alpha << "\n" << std::endl; } diff --git a/examples/ex33p.cpp b/examples/ex33p.cpp index 3231e9cd0e..0373cbda44 100644 --- a/examples/ex33p.cpp +++ b/examples/ex33p.cpp @@ -3,15 +3,21 @@ // Compile with: make ex33p // // Sample runs: mpirun -np 4 ex33p -m ../data/square-disc.mesh -alpha 0.33 -o 2 +// mpirun -np 4 ex33p -m ../data/square-disc.mesh -alpha 4.5 -o 3 +// mpirun -np 4 ex33p -m ../data/star.mesh -alpha 1.4 -o 3 // mpirun -np 4 ex33p -m ../data/star.mesh -alpha 0.99 -o 3 // mpirun -np 4 ex33p -m ../data/inline-quad.mesh -alpha 0.5 -o 3 // mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3 +// mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 2.4 -o 3 -r 4 // mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4 +// mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 1.7 -o 3 -r 5 // // Verification runs: -// mpirun -np 4ex33 -m ../data/inline-quad.mesh -ver -alpha -o 2 -r 4 -// Note: the analytic solution to this problem is u(x) = sin(pi x) sin(pi y) -// for all alpha. +// mpirun -np 4 ex33p -m ../data/inline-quad.mesh -ver -alpha 0.3 -o 2 -r 4 +// mpirun -np 4 ex33p -m ../data/inline-quad.mesh -ver -alpha 1.4 -o 3 -r 5 +// mpirun -np 4 ex33p -m ../data/inline-quad.mesh -ver -alpha 5.7 -o 2 -r 7 -no-vis +// Note: the analytic solution to this problem is u(x) = sin(pi x) sin(pi y) +// for all alpha. // // Description: // @@ -19,14 +25,14 @@ // // ( - Δ )^α u = f in Ω, u = 0 on ∂Ω, 0 < α, // -// To solve this FPDE, we multiply with ( - Δ )^(-N) where the integer -// N is given by floor(α). We obtain +// To solve this FPDE, we apply the operator ( - Δ )^(-N), where the integer +// N is given by floor(α). By doing so, we obtain // // ( - Δ )^(α-N) u = ( - Δ )^(-N) f in Ω, u = 0 on ∂Ω, 0 < α. // // We first compute the right hand side by solving the integer order PDE // -// ( - Δ )^N g = f in Ω, g = 0 on ∂Ω, +// ( - Δ )^N g = f in Ω, g = ( - Δ )^k g = 0 on ∂Ω, k = 1,..,N-1 // // The remaining FPDE is then given by // @@ -82,10 +88,6 @@ int main(int argc, char *argv[]) int num_procs = Mpi::WorldSize(); int myid = Mpi::WorldRank(); Hypre::Init(); - if (Mpi::Root()) - { - mfem::out << "\nTotal number of MPI ranks = " << num_procs << endl; - } // 1. Parse command-line options. const char *mesh_file = "../data/star.mesh"; @@ -127,11 +129,11 @@ int main(int argc, char *argv[]) // 2. Compute the rational expansion coefficients that define the // integer-order PDEs. - int power_of_laplace = floor(alpha); + const int power_of_laplace = floor(alpha); double exponent_to_approximate = alpha - power_of_laplace; bool integer_order = false; // Check if alpha is an integer or not. - if (alpha - power_of_laplace !=0) + if (abs(exponent_to_approximate) > 1e-12) { if (Mpi::Root()) { @@ -141,6 +143,10 @@ int main(int argc, char *argv[]) } ComputePartialFractionApproximation(exponent_to_approximate, coeffs, poles); + + // If the example is build without LAPACK, the exponent_to_approximate + // might be modified by the function call above. + alpha = exponent_to_approximate + power_of_laplace; } else { @@ -217,7 +223,7 @@ int main(int argc, char *argv[]) b.Assemble(); // ------------------------------------------------------------------------ - // 10. Solve the PDE -Δ ^ N g = f, i.e. compute g = (-Δ)^{-1}^N f. + // 10. Solve the PDE (-Δ)^N g = f, i.e. compute g = (-Δ)^{-1}^N f. // ------------------------------------------------------------------------ if (power_of_laplace > 0) @@ -241,17 +247,16 @@ int main(int argc, char *argv[]) k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B); HypreBoomerAMG prec; prec.SetPrintLevel(-1); - int print_level = (myid==0) ? 3 : 0; CGSolver cg(MPI_COMM_WORLD); cg.SetRelTol(1e-12); cg.SetMaxIter(2000); - cg.SetPrintLevel(print_level); + cg.SetPrintLevel(3); cg.SetPreconditioner(prec); cg.SetOperator(*Op); if (Mpi::Root()) { - mfem::out << "\nComputing -Δ ^ -" << power_of_laplace + mfem::out << "\nComputing (-Δ) ^ -" << power_of_laplace << " ( f ) " << endl; } for (int i = 0; i < power_of_laplace; i++) @@ -340,11 +345,10 @@ int main(int argc, char *argv[]) HypreBoomerAMG prec; prec.SetPrintLevel(-1); - int print_level = (myid==0) ? 3 : 0; CGSolver cg(MPI_COMM_WORLD); cg.SetRelTol(1e-12); cg.SetMaxIter(2000); - cg.SetPrintLevel(print_level); + cg.SetPrintLevel(3); cg.SetPreconditioner(prec); cg.SetOperator(*A); cg.Mult(B, X); @@ -369,8 +373,8 @@ int main(int argc, char *argv[]) oss_u.str(""); oss_u.clear(); oss_u << "Step " << progress_steps + 1 - << ": Solution of fractional PDE -Δ^" << alpha - floor(alpha) - << " u = g"; + << ": Solution of fractional PDE (-Δ)^" << alpha + << " u = f"; uout << "parallel " << num_procs << " " << myid << "\n" << "solution\n" << pmesh << u << "window_title '" << oss_u.str() << "'" From 38c60734b62494cfdc20686ef6f48b3de1f71e5e Mon Sep 17 00:00:00 2001 From: Socratis Petrides Date: Fri, 10 Jun 2022 10:37:38 -0700 Subject: [PATCH 7/9] fix typo, increase maxit in CG and remove slow sample run --- examples/ex33.cpp | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/examples/ex33.cpp b/examples/ex33.cpp index 13bca1a5e2..853a903cd3 100644 --- a/examples/ex33.cpp +++ b/examples/ex33.cpp @@ -15,7 +15,6 @@ // Verification runs: // ex33 -m ../data/inline-quad.mesh -ver -alpha 0.3 -o 2 -r 4 // ex33 -m ../data/inline-quad.mesh -ver -alpha 1.4 -o 3 -r 5 -// ex33 -m ../data/inline-quad.mesh -ver -alpha 5.7 -o 2 -r 7 -no-vis // Note: the analytic solution to this problem is u(x) = sin(pi x) sin(pi y) // for all alpha. // @@ -221,7 +220,7 @@ int main(int argc, char *argv[]) Array empty; m.FormSystemMatrix(empty, mass); - // 10.3 from the system of equations + // 10.3 Form the system of equations Vector B, X; OperatorPtr Op; k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B); @@ -231,8 +230,8 @@ int main(int argc, char *argv[]) << " ( f ) " << endl; for (int i = 0; i < power_of_laplace; i++) { - // 10.4 Solve the linear system A X = B (N times). - PCG(*Op, M, B, X, 3, 200, 1e-12, 0.0); + // 10.4 Solve the linear system Op X = B (N times). + PCG(*Op, M, B, X, 3, 300, 1e-12, 0.0); // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step if (i == power_of_laplace - 1) @@ -311,7 +310,7 @@ int main(int argc, char *argv[]) // 11.4 Solve the linear system A X = B. GSSmoother M((SparseMatrix&)(*A)); - PCG(*A, M, B, X, 3, 200, 1e-12, 0.0); + PCG(*A, M, B, X, 3, 300, 1e-12, 0.0); // 11.5 Recover the solution as a finite element grid function. a.RecoverFEMSolution(X, b, x); From 0917cebd108bedcc73423c229a12288da0863591 Mon Sep 17 00:00:00 2001 From: Socratis Petrides Date: Fri, 10 Jun 2022 11:16:39 -0700 Subject: [PATCH 8/9] adding support for 1D, 3D and AMR mesh --- examples/ex33.cpp | 56 +++++++++++++++++++++++++++++++++++++++-------- 1 file changed, 47 insertions(+), 9 deletions(-) diff --git a/examples/ex33.cpp b/examples/ex33.cpp index 853a903cd3..f362090162 100644 --- a/examples/ex33.cpp +++ b/examples/ex33.cpp @@ -7,15 +7,18 @@ // ex33 -m ../data/star.mesh -alpha 1.4 -o 3 // ex33 -m ../data/star.mesh -alpha 0.99 -o 3 // ex33 -m ../data/inline-quad.mesh -alpha 0.5 -o 3 +// ex33 -m ../data/amr-quad.mesh -alpha 1.5 -o 3 // ex33 -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3 // ex33 -m ../data/disc-nurbs.mesh -alpha 2.4 -o 3 -r 4 // ex33 -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4 // ex33 -m ../data/l-shape.mesh -alpha 1.7 -o 3 -r 5 // // Verification runs: -// ex33 -m ../data/inline-quad.mesh -ver -alpha 0.3 -o 2 -r 4 -// ex33 -m ../data/inline-quad.mesh -ver -alpha 1.4 -o 3 -r 5 -// Note: the analytic solution to this problem is u(x) = sin(pi x) sin(pi y) +// ex33 -m ../data/inline-segment.mesh -ver -alpha 1.7 -o 2 -r 2 +// ex33 -m ../data/inline-quad.mesh -ver -alpha 1.2 -o 2 -r 2 +// ex33 -m ../data/inline-hex.mesh -ver -alpha 0.3 -o 2 -r 1 +// +// Note: the analytic solution to this problem is u = ∏_{i=0}^{dim-1} sin(π x_i) // for all alpha. // // Description: @@ -169,7 +172,12 @@ int main(int argc, char *argv[]) // 7. Define diffusion coefficient, load, and solution GridFunction. auto func = [&alpha](const Vector &x) { - return pow(2*pow(M_PI,2), alpha) * sin(M_PI * x[0]) * sin(M_PI * x[1]); + double val = 1.0; + for (int i=0; iMultTranspose(B,b); + } + else + { + b = B; + } } // ------------------------------------------------------------------------ @@ -347,15 +363,37 @@ int main(int argc, char *argv[]) { auto solution = [] (const Vector &x) { - return sin(M_PI * x[0]) * sin(M_PI * x[1]); + double val = 1.0; + for (int i=0; i Date: Fri, 10 Jun 2022 11:27:07 -0700 Subject: [PATCH 9/9] parallel example 1D,3D support --- examples/ex33.cpp | 1 + examples/ex33p.cpp | 54 +++++++++++++++++++++++++++++++++++----------- 2 files changed, 43 insertions(+), 12 deletions(-) diff --git a/examples/ex33.cpp b/examples/ex33.cpp index f362090162..781437a2a9 100644 --- a/examples/ex33.cpp +++ b/examples/ex33.cpp @@ -16,6 +16,7 @@ // Verification runs: // ex33 -m ../data/inline-segment.mesh -ver -alpha 1.7 -o 2 -r 2 // ex33 -m ../data/inline-quad.mesh -ver -alpha 1.2 -o 2 -r 2 +// ex33 -m ../data/amr-quad.mesh -ver -alpha 2.6 -o 2 -r 2 // ex33 -m ../data/inline-hex.mesh -ver -alpha 0.3 -o 2 -r 1 // // Note: the analytic solution to this problem is u = ∏_{i=0}^{dim-1} sin(π x_i) diff --git a/examples/ex33p.cpp b/examples/ex33p.cpp index 0373cbda44..1a3850e9b4 100644 --- a/examples/ex33p.cpp +++ b/examples/ex33p.cpp @@ -7,16 +7,19 @@ // mpirun -np 4 ex33p -m ../data/star.mesh -alpha 1.4 -o 3 // mpirun -np 4 ex33p -m ../data/star.mesh -alpha 0.99 -o 3 // mpirun -np 4 ex33p -m ../data/inline-quad.mesh -alpha 0.5 -o 3 -// mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3 +// mpirun -np 4 ex33p -m ../data/amr-quad.mesh -alpha 1.5 -o 3 +// mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 0.33 -o 3 -r 2 // mpirun -np 4 ex33p -m ../data/disc-nurbs.mesh -alpha 2.4 -o 3 -r 4 // mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 0.33 -o 3 -r 4 // mpirun -np 4 ex33p -m ../data/l-shape.mesh -alpha 1.7 -o 3 -r 5 // // Verification runs: -// mpirun -np 4 ex33p -m ../data/inline-quad.mesh -ver -alpha 0.3 -o 2 -r 4 -// mpirun -np 4 ex33p -m ../data/inline-quad.mesh -ver -alpha 1.4 -o 3 -r 5 -// mpirun -np 4 ex33p -m ../data/inline-quad.mesh -ver -alpha 5.7 -o 2 -r 7 -no-vis -// Note: the analytic solution to this problem is u(x) = sin(pi x) sin(pi y) +// mpirun -np 4 ex33p -m ../data/inline-segment.mesh -ver -alpha 1.7 -o 2 -r 2 +// mpirun -np 4 ex33p -m ../data/inline-quad.mesh -ver -alpha 1.2 -o 2 -r 2 +// mpirun -np 4 ex33p -m ../data/amr-quad.mesh -ver -alpha 2.6 -o 2 -r 2 +// mpirun -np 4 ex33p -m ../data/inline-hex.mesh -ver -alpha 0.3 -o 2 -r 1 + +// Note: the analytic solution to this problem is u = ∏_{i=0}^{dim-1} sin(π x_i) // for all alpha. // // Description: @@ -190,7 +193,12 @@ int main(int argc, char *argv[]) // 7. Define diffusion coefficient, load, and solution GridFunction. auto func = [&alpha](const Vector &x) { - return pow(2*pow(M_PI,2), alpha) * sin(M_PI * x[0]) * sin(M_PI * x[1]); + double val = 1.0; + for (int i=0; i empty; m.FormSystemMatrix(empty, mass); - // 10.3 from the system of equations + // 10.3 Form the system of equations Vector B, X; OperatorPtr Op; k.FormLinearSystem(ess_tdof_list, g, b, Op, X, B); @@ -261,7 +269,7 @@ int main(int argc, char *argv[]) } for (int i = 0; i < power_of_laplace; i++) { - // 10.4 Solve the linear system A X = B (N times). + // 10.4 Solve the linear system Op X = B (N times). cg.Mult(B, X); // 10.5 Visualize the solution g of -Δ ^ N g = f in the last step if (i == power_of_laplace - 1) @@ -390,17 +398,39 @@ int main(int argc, char *argv[]) { auto solution = [] (const Vector &x) { - return sin(M_PI * x[0]) * sin(M_PI * x[1]); + double val = 1.0; + for (int i=0; i