### Calculate L2 error norm in C++

392
views
0
9 months ago by

I recently switched over from Python to C++, and am trying to calculate the L2 error norm:

L2error = sqrt(assemble(inner(u_h - u_exact, u_h - u_exact)*dx))

I don't directly see any C++ examples that explicitly show how to do the above. Suppose I am solving the poisson problem in 3D. I have the following analytical solution and source term:

class ExactSolution : public Expression
{
void eval(Array<double>& values, const Array<double>& x) const
{
values[0] = sin(2*M_PI*x[0])*sin(2*M_PI*x[1])*sin(2*M_PI*x[2]);
}
};

class Source : public Expression
{
void eval(Array<double>& values, const Array<double>& x) const
{
values[0] = 12*M_PI*M_PI*sin(2*M_PI*x[0])*sin(2*M_PI*x[1])*sin(2*M_PI*x[2]);
}
};

I also have the following UFL file:

element = FiniteElement("Lagrange", tetrahedron, 1)
u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)
L = f*v*dx

and the main cpp code:

class DirichletBoundary : public SubDomain
{
bool inside(const Array<double>& x, bool on_boundary) const
{
return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS or x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS or x[2] < DOLFIN_EPS or x[2] > 1.0 - DOLFIN_EPS;
}
};

int main(int argc, char* argv[])
{
// Initialize stuff
double tstart,tend;
int rank, size;
int seed;
seed = atoi(argv[1]);

// Create mesh and function space
auto mesh = std::make_shared<UnitCubeMesh>(seed,seed,seed);
auto V = std::make_shared<CG1::FunctionSpace>(mesh);
MPI_Comm_size(MPI_COMM_WORLD,&size);
MPI_Comm_rank(MPI_COMM_WORLD,&rank);

// Define boundary condition
auto u0 = std::make_shared<Constant>(0.0);
auto boundary = std::make_shared<DirichletBoundary>();
auto bc = std::make_shared<DirichletBC>(V, u0, boundary);

// Define variational forms
CG1::BilinearForm a(V, V);
CG1::LinearForm L(V);
auto f = std::make_shared<Source>();
L.f = f;

// Define linear algebra
auto A = std::make_shared<Matrix>();
Vector b;
KrylovSolver solver("gmres", "petsc_amg");
solver.parameters["relative_tolerance"] = 1.0e-7;
PETScOptions::set("ksp_converged_reason");
Function u(V);

// Show problem size
auto usize = u.vector();
if (!rank) {
PetscPrintf(MPI_COMM_WORLD,"\tDegrees-of-freedom: %d\n", usize->size());
}

// Compute solution
tstart = MPI_Wtime();
assemble_system(*A, b, a, L, {bc});
solver.set_operators(A,A);
solver.solve(*u.vector(),b);
tend = MPI_Wtime();

// Show time-to-solution
if (!rank) {
PetscPrintf(MPI_COMM_WORLD,"\tWall-clock time: %1.3e seconds\n", tend-tstart);
}

// Save solution in VTK format
//File file("poisson.pvd");
//file << u;

return 0;
}

Community: FEniCS Project

2
9 months ago by

Okay I found one solution, I needed this additional UFL:

element = FiniteElement("Lagrange", tetrahedron, 1)
element2 = FiniteElement("Lagrange", tetrahedron, 4)
ComputedSolution = Coefficient(element)
ExactSolution = Coefficient(element2)
a = (ComputedSolution - ExactSolution)*(ComputedSolution - ExactSolution)*dx

And the modified main.cpp file:

int main(int argc, char* argv[])
{
// Initialize stuff
double tstart,tend;
int rank, size;
int seed;
seed = atoi(argv[1]);

// Create mesh and function space
auto mesh = std::make_shared<UnitCubeMesh>(seed,seed,seed);
auto V = std::make_shared<CG1::FunctionSpace>(mesh);
MPI_Comm_size(MPI_COMM_WORLD,&size);
MPI_Comm_rank(MPI_COMM_WORLD,&rank);

// Define boundary condition
auto u0 = std::make_shared<ExactSolution>();
auto boundary = std::make_shared<DirichletBoundary>();
auto bc = std::make_shared<DirichletBC>(V, u0, boundary);

// Define variational forms
CG1::BilinearForm a(V, V);
CG1::LinearForm L(V);
auto f = std::make_shared<Source>();
L.f = f;

// Define linear algebra
auto A = std::make_shared<Matrix>();
Vector b;
KrylovSolver solver("gmres", "petsc_amg");
solver.parameters["relative_tolerance"] = 1.0e-7;
PETScOptions::set("ksp_converged_reason");
//Function u(V);
auto u = std::make_shared<Function>(V);

// Show problem size
auto usize = u->vector();
if (!rank) {
PetscPrintf(MPI_COMM_WORLD,"\tDegrees-of-freedom: %d\n", usize->size());
}

// Compute solution
tstart = MPI_Wtime();
assemble_system(*A, b, a, L, {bc});
solver.set_operators(A,A);
solver.solve(*u->vector(),b);
tend = MPI_Wtime();

// Compute L2 error
L2error::Functional L2_errorform(mesh);
L2_errorform.ExactSolution = u0;
L2_errorform.ComputedSolution = u;
double l2errornorm = sqrt(assemble(L2_errorform));

// Show time-to-solution
if (!rank) {
PetscPrintf(MPI_COMM_WORLD,"\tWall-clock time: %1.3e seconds\n", tend-tstart);
PetscPrintf(MPI_COMM_WORLD,"\tL2 error norm: %1.3e\n", l2errornorm);
}

// Save solution in VTK format
//File file("poisson.pvd");
//file << u;

return 0;
}

Not sure this is the best or most efficient way of calculating the error norm, but if anyone has any alternative, I am all ears.