1 // Copyright (c) 2017-2026, Lawrence Livermore National Security, LLC and other CEED contributors.
2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
3 //
4 // SPDX-License-Identifier: BSD-2-Clause
5 //
6 // This file is part of CEED: http://github.com/ceed
7
8 // libCEED + MFEM Example: BP3
9 //
10 // This example illustrates a simple usage of libCEED with the MFEM (mfem.org) finite element library.
11 //
12 // The example reads a mesh from a file and solves a linear system with a diffusion stiffness matrix (with a prescribed analytic solution, provided by
13 // the function 'solution'). The diffusion matrix is expressed as a new class, CeedDiffusionOperator, derived from mfem::Operator. Internally,
14 // CeedDiffusionOperator uses a CeedOperator object constructed based on an mfem::FiniteElementSpace. All libCEED objects use a Ceed logical device
15 // object constructed based on a command line argument. (-ceed).
16 //
17 // The linear system is inverted using the conjugate gradients algorithm corresponding to CEED BP3, see http://ceed.exascaleproject.org/bps.
18 // Arbitrary mesh and solution orders in 1D, 2D and 3D are supported from the same code.
19 //
20 // Build with:
21 //
22 // make bp3 [MFEM_DIR=</path/to/mfem>] [CEED_DIR=</path/to/libceed>]
23 //
24 // Sample runs:
25 //
26 // ./bp3
27 // ./bp3 -ceed /cpu/self
28 // ./bp3 -ceed /gpu/cuda
29 // ./bp3 -m ../../../mfem/data/fichera.mesh -o 4
30 // ./bp3 -m ../../../mfem/data/square-disc-nurbs.mesh -o 6
31 // ./bp3 -m ../../../mfem/data/inline-segment.mesh -o 8
32
33 /// @file
34 /// MFEM diffusion operator based on libCEED
35
36 #include "bp3.hpp"
37
38 #include <ceed.h>
39
40 #include <mfem.hpp>
41
42 /// Exact solution
solution(const mfem::Vector & pt)43 double solution(const mfem::Vector &pt) {
44 static const double x[3] = {-0.32, 0.15, 0.24};
45 static const double k[3] = {1.21, 1.45, 1.37};
46 double val = sin(M_PI * (x[0] + k[0] * pt(0)));
47 for (int d = 1; d < pt.Size(); d++) val *= sin(M_PI * (x[d] + k[d] * pt(d)));
48 return val;
49 }
50
51 /// Right-hand side
rhs(const mfem::Vector & pt)52 double rhs(const mfem::Vector &pt) {
53 static const double x[3] = {-0.32, 0.15, 0.24};
54 static const double k[3] = {1.21, 1.45, 1.37};
55 double f[3], l[3], val, lap;
56 f[0] = sin(M_PI * (x[0] + k[0] * pt(0)));
57 l[0] = M_PI * M_PI * k[0] * k[0] * f[0];
58 val = f[0];
59 lap = l[0];
60 for (int d = 1; d < pt.Size(); d++) {
61 f[d] = sin(M_PI * (x[d] + k[d] * pt(d)));
62 l[d] = M_PI * M_PI * k[d] * k[d] * f[d];
63 lap = lap * f[d] + val * l[d];
64 val = val * f[d];
65 }
66 return lap;
67 }
68
69 //TESTARGS -ceed {ceed_resource} -t -no-vis --size 2000
main(int argc,char * argv[])70 int main(int argc, char *argv[]) {
71 // 1. Parse command-line options.
72 const char *ceed_spec = "/cpu/self";
73 #ifndef MFEM_DIR
74 const char *mesh_file = "../../../mfem/data/star.mesh";
75 #else
76 const char *mesh_file = MFEM_DIR "/data/star.mesh";
77 #endif
78 int order = 2;
79 bool visualization = true;
80 bool test = false;
81 double max_nnodes = 50000;
82
83 mfem::OptionsParser args(argc, argv);
84 args.AddOption(&ceed_spec, "-c", "-ceed", "Ceed specification.");
85 args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
86 args.AddOption(&order, "-o", "--order", "Finite element order (polynomial degree).");
87 args.AddOption(&max_nnodes, "-s", "--size", "Maximum size (number of DoFs)");
88 args.AddOption(&visualization, "-vis", "--visualization", "-no-vis", "--no-visualization", "Enable or disable GLVis visualization.");
89 args.AddOption(&test, "-t", "--test", "-no-test", "--no-test", "Enable or disable test mode.");
90 args.Parse();
91 if (!args.Good()) {
92 args.PrintUsage(std::cout);
93 return 1;
94 }
95 if (!test) {
96 args.PrintOptions(std::cout);
97 }
98
99 // 2. Initialize a Ceed device object using the given Ceed specification.
100 Ceed ceed;
101 CeedInit(ceed_spec, &ceed);
102
103 // 3. Read the mesh from the given mesh file.
104 mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
105 int dim = mesh->Dimension();
106
107 // 4. Refine the mesh to increase the resolution.
108 // In this example we do 'ref_levels' of uniform refinement.
109 // We choose 'ref_levels' to be the largest number that gives a final system with no more than 50,000 unknowns, approximately.
110 {
111 int ref_levels = (int)floor((log(max_nnodes / mesh->GetNE()) - dim * log(order)) / log(2.) / dim);
112 for (int l = 0; l < ref_levels; l++) {
113 mesh->UniformRefinement();
114 }
115 }
116 if (mesh->GetNodalFESpace() == NULL) {
117 mesh->SetCurvature(1, false, -1, mfem::Ordering::byNODES);
118 }
119 if (mesh->NURBSext) {
120 mesh->SetCurvature(order, false, -1, mfem::Ordering::byNODES);
121 }
122
123 // 5. Define a finite element space on the mesh.
124 // Here we use continuous Lagrange finite elements of the specified order.
125 MFEM_VERIFY(order > 0, "invalid order");
126 mfem::FiniteElementCollection *fec = new mfem::H1_FECollection(order, dim);
127 mfem::FiniteElementSpace *fespace = new mfem::FiniteElementSpace(mesh, fec);
128 if (!test) {
129 std::cout << "Number of finite element unknowns: " << fespace->GetTrueVSize() << std::endl;
130 }
131
132 mfem::FunctionCoefficient sol_coeff(solution);
133 mfem::Array<int> ess_tdof_list;
134 mfem::GridFunction sol(fespace);
135 if (mesh->bdr_attributes.Size()) {
136 mfem::Array<int> ess_bdr(mesh->bdr_attributes.Max());
137 ess_bdr = 1;
138 fespace->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
139 sol.ProjectBdrCoefficient(sol_coeff, ess_bdr);
140 }
141
142 // 6. Construct a rhs vector using the linear form f(v) = (rhs, v), where v is a test function.
143 mfem::LinearForm b(fespace);
144 mfem::FunctionCoefficient rhs_coeff(rhs);
145 b.AddDomainIntegrator(new mfem::DomainLFIntegrator(rhs_coeff));
146 b.Assemble();
147
148 // 7. Construct a CeedDiffusionOperator utilizing the 'ceed' device and using the 'fespace' object to extract data needed by the Ceed objects.
149 CeedDiffusionOperator diff(ceed, fespace);
150
151 mfem::Operator *D;
152 mfem::Vector X, B;
153 diff.FormLinearSystem(ess_tdof_list, sol, b, D, X, B);
154
155 // 8. Solve the discrete system using the conjugate gradients (CG) method.
156 mfem::CGSolver cg;
157 cg.SetRelTol(1e-6);
158 cg.SetMaxIter(1000);
159 if (test) {
160 cg.SetPrintLevel(0);
161 } else {
162 cg.SetPrintLevel(3);
163 }
164 cg.SetOperator(*D);
165
166 cg.Mult(B, X);
167
168 // 9. Compute and print the L2 norm of the error.
169 double err_l2 = sol.ComputeL2Error(sol_coeff);
170 if (!test) {
171 std::cout << "L2 projection error: " << err_l2 << std::endl;
172 } else {
173 if (fabs(sol.ComputeL2Error(sol_coeff)) > 2e-3) {
174 std::cout << "Error too large: " << err_l2 << std::endl;
175 }
176 }
177
178 // 10. Open a socket connection to GLVis and send the mesh and solution for visualization.
179 if (visualization) {
180 char vishost[] = "localhost";
181 int visport = 19916;
182 mfem::socketstream sol_sock(vishost, visport);
183 sol_sock.precision(8);
184 sol_sock << "solution\n" << *mesh << sol << std::flush;
185 }
186
187 // 11. Free memory and exit.
188 delete fespace;
189 delete fec;
190 delete mesh;
191 delete D;
192 CeedDestroy(&ceed);
193 return 0;
194 }
195