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:
f = -div(grad(u_exact))
# Solve the weak problem with FEA, using that f:
u = TrialFunction(V)
v = TestFunction(V)
bc = DirichletBC(V,0,"on_boundary")
uh = Function(V)
solve(dot(grad(u),grad(v))*dx==f*v*dx,uh,bc)
# Check error:
e = uh-u_exact
print(sqrt(assemble(e*e*dx)))
This can be iterated to more complicated problems, although you may need to think a little harder to satisfy constraints. E.g., if you need a solenoidal exact velocity solution for the “p-Stokes” problem, you can take the curl of a vector potential.
Hi, thanks for your response!
Exactly how would you proceed to find a manufactured solution to an equation which I do not know if it even has analytical solutions? Or do you mean that it only has to fulfill some boundary condition + additional constraint?
If you have the strong form of the problem, whatever solution you manufacture the corresponding source term for will, by construction, solve the PDE system as long as it satisfies the boundary conditions and constraints. There may not be a solution for every source term, but there is a source term for every solution. For complicated nonlinear problems, you might not know whether the manufactured solution is stable or unique, which may or may not become an issue in practice. (For a concrete example, I was using manufactured solutions to test a numerical method for geometrically-nonlinear dynamic fluid–membrane interaction problem a few years ago, and I found that I needed to consider solutions with a tensile pre-stress added to the membrane, to prevent buckling instabilities.)