# How to approximate the error without exact solution?

Quite new to FEniCS so I’m sorry if the question is stupid.

I’m solving a fluid mechanics problem for non-Newtonian fluids, thus solving the p-Stokes equations (related to the Navier-Stokes, but with shear rate dependent viscosity). In the implementation there is a regularization factor \epsilon in the expression for the viscosity which is supposed to be very small in order to approximate the real model. I have no exact solution to compare with but I still want to calculate the evolution of the error with increasing mesh size, and to be able to compare the error for different values on \epsilon. I can’t use simply the residual in this case, since a too large value on \epsilon can give a small residual error but essentially solving the wrong problem.

How do I proceed with this? My first though was to solve the problem on a very fine grid and for a very small regularization factor, and then interpolate my numerical solution onto the function space of the more accurate solution, then simply take the error norm. This produces the following error:

“”"
Error: Unable to create mesh entity.
*** Reason: Mesh entity index -1 out of range [0, 1600] for entity of dimension 2.
*** Where: This error was encountered inside MeshEntity.cpp.
“”"

If you know the strong form of the PDE, you can use the method of manufactured solutions to turn whatever function you want into the exact solution. FEniCS UFL makes this very easy. For instance, for the Poisson problem:

from dolfin import *
N = 32
k = 1
mesh = UnitIntervalMesh(N)
V = FunctionSpace(mesh,"CG",k)
x = SpatialCoordinate(mesh)

# Choose some arbitrary exact solution, satisfying
# the BCs (and any other constraints):
u_exact = sin(pi*x[0])

# Generate the corresponding source term from
# the strong form of the PDE, which can be
# automated using UFL:

# Solve the weak problem with FEA, using that f:
u = TrialFunction(V)
v = TestFunction(V)
bc = DirichletBC(V,0,"on_boundary")
uh = Function(V)