VMS stabilization of 2D steady advection-diffusion fails
11 weeks ago by
As I was implementing VMS stabilization of a related problem I ran into some issues. I reduced it to the following MWE. It concerns a simple steady 2D advection-diffusion problem on a square mesh with a diagonally directed advection vector. The problem occurs when I use higher order basis functions (p=2).
from dolfin import * # Parameters nu = Constant(0.01) #diffusion coefficient a_mag = 0.8 #advection magnitude a = Constant((a_mag*0.5**0.5,a_mag*0.5**0.5)) #advection vector VMS = False #toggle to include or exclude VMS stabilization # Create mesh and define function space mesh = RectangleMesh(Point(-0.5,-0.5),Point(0.5,0.5),10,10) V = FunctionSpace(mesh, "CG", 2) # Boundary conditions def bdy_outer(x,on_boundary): return on_boundary bc_outer = Expression("x+x", degree=1,domain=mesh) bcs = [DirichletBC(V, bc_outer, bdy_outer)] # Stability parameter h = CellSize(mesh) tau = pow( (pow(4.*nu/pow(h,2),2)+pow(2*a_mag/h,2)) ,-0.5) # Define unknown and test function v = TestFunction(V) u = TrialFunction(V) # Define variational forms F = -inner(grad(v),a)*u*dx+ nu*inner(grad(u), grad(v))*dx if VMS: F += -tau*( -inner(grad(v),a)-nu*div(grad(v)) )*( inner(grad(u),a)-nu*div(grad(u)))*dx # Solving the problem u = Function(V) problem = LinearVariationalProblem(lhs(F),rhs(F), u, bcs) solver = LinearVariationalSolver(problem) solver.solve() # Projecting on a refined mesh for visualization meshRefine = RectangleMesh(Point(-0.5,-0.5),Point(0.5,0.5),100,100) V2 = FunctionSpace(meshRefine, "CG", 1) u2 = project(u,V2) u2.rename("u", "u") # Create files for storing results file = File("solution.pvd") file << u2
The solution completely blows up when the VMS stabilization is included. This can be fixed by switching to SUPG or GLS, but im actually interested in the VMS part. Still, this illustrates that the (sign of the) Laplace operator on the test function is the root of the problem.
- Am I making a mistake with my FEniCS implementation?
- Am I making a mistake with my VMS formulation?
- Is VMS unstuited for the current problem?
Even though I've checked and rechecked the first two points, I can't believe that it would be the third. This MWE should be the simplest example of true VMS (in the sense that p=1 reduces to SUPG and 1D could be reconstructed analytically). I would find it very odd if VMS doesn't work for this case.
I would appreciate any help or advice,
Community: FEniCS Project
11 weeks ago by
# Define variational forms F = -inner(grad(v),a)*u*dx+ nu*inner(grad(u), grad(v))*dx if VMS: #F += -tau*( -inner(grad(v),a)-nu*div(grad(v)) )*( inner(grad(u),a)-nu*div(grad(u)))*dx def L(u): return inner(grad(u),a)-nu*div(grad(u)) def Ladj_formal(u): return -inner(grad(u),a)-nu*div(grad(u)) def badj(u,n): return nu*inner(grad(u),n) res = L(u) # - f uPrime = -tau*res # Equivalent to what was there originally: F += Ladj_formal(v)*uPrime*dx # Boundary terms from deriving the adjoint operator # in the non-smooth FEM setting: n = FacetNormal(mesh) F += badj(v('+'),n('+'))*uPrime('+')*dS \ + badj(v('-'),n('-'))*uPrime('-')*dS \ + badj(v,n)*uPrime*ds
Please login to add an answer/comment or follow this question.