Skip to content

Commit 988aaf7

Browse files
skang67skang67
committed
Working code for element_least_squares.
1 parent cf64636 commit 988aaf7

File tree

1 file changed

+88
-88
lines changed

1 file changed

+88
-88
lines changed

‎miniapps/tools/reconstruction.cpp‎

Lines changed: 88 additions & 88 deletions
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ int main(int argc, char* argv[])
2222
Mpi::Init(argc, argv);
2323

2424
// Default command-line options
25-
std::string lor_method = "optimization"; // "optimization" or "l2_projection"
25+
std::string lor_method = "element_least_squares"; // "element_least_squares" or "l2_projection"
2626
int refinement_levels = 0;
2727
int order_lo = 0;
2828
int order_ho = 3;
@@ -73,7 +73,7 @@ int main(int argc, char* argv[])
7373

7474
// Parse options
7575
OptionsParser args(argc, argv);
76-
args.AddOption(&lor_method, "-m", "--method", "LOR method.");
76+
args.AddOption(&lor_method, "-m", "--method", "LOR method: \"element_least_squares\" or \"l2_projection\".");
7777
args.AddOption(&refinement_levels, "-r", "--refine",
7878
"Number of times to refine the mesh uniformly.");
7979
args.AddOption(&order_ho, "-ho", "--order_ho",
@@ -105,9 +105,9 @@ int main(int argc, char* argv[])
105105
const int num_x = 8;
106106
const int num_y = 8;
107107

108-
// std::cout << "LOR method: " << lor_method << ", (lor_method == l2_projection) = "<< (lor_method == "l2_projection")<< std::endl;
108+
std::cout << "LOR method: " << lor_method << std::endl;
109109

110-
if (lor_method == "optimization") {
110+
// if (lor_method == "element_least_squares") {
111111

112112
Mesh serial_mesh = Mesh::MakeCartesian2D(num_x, num_y, Element::QUADRILATERAL);
113113
for (int i = 0; i < refinement_levels; i++) {
@@ -117,23 +117,23 @@ int main(int argc, char* argv[])
117117
mesh = ParMesh(MPI_COMM_WORLD, serial_mesh);
118118
serial_mesh.Clear();
119119

120-
} else if (lor_method == "l2_projection") {
121-
122-
order_im = order_ho;
123-
lref = order_im + 1;
124-
MFEM_VERIFY((num_x % lref) == 0 && (num_y % lref) == 0,
125-
"For l2_projection, lref = order_im (=order_ho) + 1 must divide both num_x and num_y.");
126-
int num_x_im = num_x / lref;
127-
int num_y_im = num_y / lref;
128-
Mesh mesh_im_serial = Mesh::MakeCartesian2D(num_x_im, num_y_im, Element::QUADRILATERAL);
129-
for (int i = 0; i < refinement_levels; i++) {
130-
mesh_im_serial.UniformRefinement();
131-
}
132-
mesh_im_serial.EnsureNCMesh();
133-
mesh_im = ParMesh(MPI_COMM_WORLD, mesh_im_serial);
134-
Mesh mesh_refined = Mesh::MakeRefined(mesh_im_serial, lref, BasisType::ClosedUniform);
135-
mesh = ParMesh(MPI_COMM_WORLD, mesh_refined); // GaussLobatto, ClosedUniform
136-
}
120+
// } else if (lor_method == "l2_projection") {
121+
122+
// order_im = order_ho;
123+
// lref = order_im + 1;
124+
// MFEM_VERIFY((num_x % lref) == 0 && (num_y % lref) == 0,
125+
// "For l2_projection, lref = order_im (=order_ho) + 1 must divide both num_x and num_y.");
126+
// int num_x_im = num_x / lref;
127+
// int num_y_im = num_y / lref;
128+
// Mesh mesh_im_serial = Mesh::MakeCartesian2D(num_x_im, num_y_im, Element::QUADRILATERAL);
129+
// for (int i = 0; i < refinement_levels; i++) {
130+
// mesh_im_serial.UniformRefinement();
131+
// }
132+
// mesh_im_serial.EnsureNCMesh();
133+
// mesh_im = ParMesh(MPI_COMM_WORLD, mesh_im_serial);
134+
// Mesh mesh_refined = Mesh::MakeRefined(mesh_im_serial, lref, BasisType::ClosedUniform);
135+
// mesh = ParMesh(MPI_COMM_WORLD, mesh_refined); // GaussLobatto, ClosedUniform
136+
// }
137137
int dim = mesh.Dimension();
138138

139139
// create FEM things
@@ -150,8 +150,8 @@ int main(int argc, char* argv[])
150150
ParGridFunction u_lo(&fespace_lo);
151151
ParGridFunction u_hi(&fespace_hi);
152152

153-
if (lor_method == "optimization")
154-
{
153+
// if (lor_method == "element_least_squares")
154+
// {
155155
FiniteElementCollection *fec_exact;
156156
fec_exact = new L2_FECollection(order_ho, dim);
157157
ParFiniteElementSpace fespace_exact(&mesh, fec_exact);
@@ -164,71 +164,71 @@ int main(int argc, char* argv[])
164164

165165
// compute reconstruction
166166
L2Reconstruction(u_lo, u_hi);
167-
} else if (lor_method == "l2_projection") {
168-
FiniteElementCollection *fec_im;
169-
fec_im = new H1_FECollection(order_im, dim); // Both L2 and H1 give the same convergence.
170-
ParFiniteElementSpace fespace_im(&mesh_im, fec_im);
171-
ParGridFunction u_im(&fespace_im);
172-
173-
BilinearForm M_lo(&fespace_lo);
174-
M_lo.AddDomainIntegrator(new MassIntegrator);
175-
M_lo.Assemble();
176-
M_lo.Finalize();
177-
178-
BilinearForm M_im(&fespace_im);
179-
M_im.AddDomainIntegrator(new MassIntegrator);
180-
M_im.Assemble();
181-
M_im.Finalize();
182-
183-
BilinearForm M_hi(&fespace_hi);
184-
M_hi.AddDomainIntegrator(new MassIntegrator);
185-
M_hi.Assemble();
186-
M_hi.Finalize();
187-
188-
// Set up the right-hand side vector for the exact solution
189-
LinearForm b_lo(&fespace_lo);
190-
// b_lo.AddDomainIntegrator(new DomainLFIntegrator(RHO));
191-
DomainLFIntegrator *lf_integ = new DomainLFIntegrator(u_function_exact);
192-
const IntegrationRule &ir_rhs = IntRules.Get(fespace_lo.GetFE(0)->GetGeomType(),
193-
order_ho+1);
194-
lf_integ->SetIntRule(&ir_rhs);
195-
b_lo.AddDomainIntegrator(lf_integ);
196-
b_lo.Assemble();
197-
198-
GridTransfer *gt1 = nullptr;
199-
GridTransfer *gt2 = nullptr;
200-
gt1 = new L2ProjectionGridTransfer(fespace_im, fespace_lo); // (dom_fes_, ran_fes_)
201-
gt2 = new L2ProjectionGridTransfer(fespace_im, fespace_hi);
202-
203-
// Configure element assembly for device acceleration
204-
gt1->UseEA(use_ea);
205-
gt2->UseEA(use_ea);
206-
207-
const Operator &P1 = gt1->BackwardOperator(); // Prolongation 1 (LO->IM)
208-
const Operator &P2 = gt2->ForwardOperator(); // Prolongation 2 (IM->HO)
209-
210-
const Operator &R1 = gt1->ForwardOperator(); // Restriction 1 (IM->LO)
211-
const Operator &R2 = gt2->BackwardOperator(); // Restriction 2 (HO->IM)
212-
213-
// STEP1: L2 projection of RHO onto u_lo
214-
SparseMatrix &M_mat_lo = M_lo.SpMat();
215-
CGSolver cg;
216-
cg.SetOperator(M_mat_lo);
217-
cg.SetRelTol(1e-16);
218-
cg.SetMaxIter(1000);
219-
cg.SetPrintLevel(0);
220-
u_lo = 0.0;
221-
cg.Mult(b_lo, u_lo); // Solve: M * u_lo = b_lo
222-
u_lo.SetTrueVector();
223-
u_lo.SetFromTrueVector();
224-
225-
// STEP2: Prolongation 1 (LO->IM)
226-
P1.Mult(u_lo, u_im); // u_im = P1 * u_lo
227-
228-
// STEP3: Prolongation 2 (IM->HO)
229-
P2.Mult(u_im, u_hi); // u_hi = P2 * u_im
230-
231-
}
167+
// } else if (lor_method == "l2_projection") {
168+
// FiniteElementCollection *fec_im;
169+
// fec_im = new H1_FECollection(order_im, dim); // Both L2 and H1 give the same convergence.
170+
// ParFiniteElementSpace fespace_im(&mesh_im, fec_im);
171+
// ParGridFunction u_im(&fespace_im);
172+
173+
// BilinearForm M_lo(&fespace_lo);
174+
// M_lo.AddDomainIntegrator(new MassIntegrator);
175+
// M_lo.Assemble();
176+
// M_lo.Finalize();
177+
178+
// BilinearForm M_im(&fespace_im);
179+
// M_im.AddDomainIntegrator(new MassIntegrator);
180+
// M_im.Assemble();
181+
// M_im.Finalize();
182+
183+
// BilinearForm M_hi(&fespace_hi);
184+
// M_hi.AddDomainIntegrator(new MassIntegrator);
185+
// M_hi.Assemble();
186+
// M_hi.Finalize();
187+
188+
// // Set up the right-hand side vector for the exact solution
189+
// LinearForm b_lo(&fespace_lo);
190+
// // b_lo.AddDomainIntegrator(new DomainLFIntegrator(RHO));
191+
// DomainLFIntegrator *lf_integ = new DomainLFIntegrator(u_function_exact);
192+
// const IntegrationRule &ir_rhs = IntRules.Get(fespace_lo.GetFE(0)->GetGeomType(),
193+
// order_ho+1);
194+
// lf_integ->SetIntRule(&ir_rhs);
195+
// b_lo.AddDomainIntegrator(lf_integ);
196+
// b_lo.Assemble();
197+
198+
// GridTransfer *gt1 = nullptr;
199+
// GridTransfer *gt2 = nullptr;
200+
// gt1 = new L2ProjectionGridTransfer(fespace_im, fespace_lo); // (dom_fes_, ran_fes_)
201+
// gt2 = new L2ProjectionGridTransfer(fespace_im, fespace_hi);
202+
203+
// // Configure element assembly for device acceleration
204+
// gt1->UseEA(use_ea);
205+
// gt2->UseEA(use_ea);
206+
207+
// const Operator &P1 = gt1->BackwardOperator(); // Prolongation 1 (LO->IM)
208+
// const Operator &P2 = gt2->ForwardOperator(); // Prolongation 2 (IM->HO)
209+
210+
// const Operator &R1 = gt1->ForwardOperator(); // Restriction 1 (IM->LO)
211+
// const Operator &R2 = gt2->BackwardOperator(); // Restriction 2 (HO->IM)
212+
213+
// // STEP1: L2 projection of RHO onto u_lo
214+
// SparseMatrix &M_mat_lo = M_lo.SpMat();
215+
// CGSolver cg;
216+
// cg.SetOperator(M_mat_lo);
217+
// cg.SetRelTol(1e-16);
218+
// cg.SetMaxIter(1000);
219+
// cg.SetPrintLevel(0);
220+
// u_lo = 0.0;
221+
// cg.Mult(b_lo, u_lo); // Solve: M * u_lo = b_lo
222+
// u_lo.SetTrueVector();
223+
// u_lo.SetFromTrueVector();
224+
225+
// // STEP2: Prolongation 1 (LO->IM)
226+
// P1.Mult(u_lo, u_im); // u_im = P1 * u_lo
227+
228+
// // STEP3: Prolongation 2 (IM->HO)
229+
// P2.Mult(u_im, u_hi); // u_hi = P2 * u_im
230+
231+
// }
232232

233233
// evaluate reconstruction
234234
char vishost[] = "localhost";

0 commit comments

Comments
 (0)