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)
a = inner(grad(u), grad(v))*dx
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;
}

Any help or advice appreciated!

Community: FEniCS Project

1 Answer


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.

Please login to add an answer/comment or follow this question.

Similar posts:
Search »