Hi, I have to solve Stokes problem on a cube domain with a sphere obstacle in the middle of the cube. When I don't use mesh adaptivity, the execution of code on a coarse mesh ends succesfully. The code is:
#Coarse problem
from dolfin import *
import numpy as np
parameters["allow_extrapolation"]= True
#Set test variables
eta = 0.001;
#Inflow pressure
p1 = 1;
# Load mesh
mesh = Mesh("./Mesh/Sphere.xml")
boundaries = MeshFunction("size_t", mesh, "./Mesh/Sphere_facet_region.xml")
# Define function spaces
V = VectorFunctionSpace (mesh, "CG", 2)
Q = FunctionSpace (mesh, "CG", 1)
W = MixedFunctionSpace([V,Q])
n = FacetNormal(mesh)
ds = Measure("ds", subdomain_data=boundaries)
# Create functions for boundary conditions
zero = Constant(0.0)
collector_bc = Constant((0,0,0))
# Boundary condition for velocity on extern boundaries
bc3 = DirichletBC(W.sub(0).sub(2), zero, boundaries, 3) #front
bc4 = DirichletBC(W.sub(0).sub(2), zero, boundaries, 4) #back
bc5 = DirichletBC(W.sub(0).sub(1), zero, boundaries, 5) #upper
bc6 = DirichletBC(W.sub(0).sub(1), zero, boundaries, 6) #lower
# Boundary condition for velocity on collector boundary
bc7 = DirichletBC(W.sub(0), collector_bc, boundaries, 7) #collector
# Collect boundary conditions
bcs = [bc3,bc4,bc5,bc6,bc7]
# Define variational problem
v, q = TestFunctions(W)
u, p = TrialFunctions(W)
a = eta*inner(grad(u), grad(v))*dx - div(v)*p*dx + q*div(u)*dx
L = - p1*inner(v,n)*ds(1)
# Compute solution
w = Function(W)
solve(a == L, w, bcs, solver_parameters={'linear_solver':'mumps'})
Now I want use adaptivity to improve coarse mesh, but when I run the code below I have this error
Traceback (most recent call last):
File "adaptive_problem.py", line 108, in
solve(F == 0, w, bcs, tol=1.0e-1, M=M)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 265, in solve
_solve_varproblem_adaptive(*args, **kwargs)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/solving.py", line 361, in _solve_varproblem_adaptive
solver.solve(tol)
File "/usr/lib/python2.7/dist-packages/dolfin/fem/adaptivesolving.py", line 123, in solve
cpp.AdaptiveNonlinearVariationalSolver.solve(self, tol)
RuntimeError:
*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
*** fenics@fenicsproject.org
*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.
*** -------------------------------------------------------------------------
*** Error: Unable to extract subsystem of finite element.
*** Reason: Requested subsystem (2) out of range [0, 2).
*** Where: This error was encountered inside FiniteElement.cpp.
*** Process: unknown
*** DOLFIN version: 1.5.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------
Code of adaptive problem:
#Adaptive Problem
from dolfin import *
import numpy as np
parameters["allow_extrapolation"]= True
parameters["refinement_algorithm"] = "plaza"
#Set test variables
eta = 0.001;
b = 10.0;
radius = 3.0;
#Inflow pressure
p1 = 1;
# Load mesh
mesh = Mesh("./Mesh/Sphere.xml")
class Inflow(SubDomain):
def inside(self,x,on_boundary):
return near(x[0], -b) and on_boundary
class Outflow(SubDomain):
def inside(self,x,on_boundary):
return near(x[0], b) and on_boundary
class UpperLower(SubDomain):
def inside(self,x,on_boundary):
return near(x[1], -b) or near(x[1], b) and on_boundary
class FrontBack(SubDomain):
def inside(self,x,on_boundary):
return near(x[2], -b) or near(x[2], b) and on_boundary
class Collector(SubDomain):
def inside(self, x, on_boundary):
r = sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2])
return r < radius + DOLFIN_EPS and on_boundary
# Define function spaces
V = VectorFunctionSpace (mesh, "CG", 2)
Q = FunctionSpace (mesh, "CG", 1)
W = V * Q
n = FacetNormal(mesh)
boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
boundaries.set_all(0)
inflow = Inflow()
inflow.mark(boundaries, 1)
outflow = Outflow()
outflow.mark(boundaries, 2)
upperlower = UpperLower()
upperlower.mark(boundaries, 3)
frontback = FrontBack()
frontback.mark(boundaries, 4)
collector = Collector()
collector.mark(boundaries, 5)
ds = Measure("ds")[boundaries]
# Create functions for boundary conditions
zero = Constant(0.0)
collector_bc = Constant((0,0,0))
# Boundary condition for velocity on extern boundaries
bc3 = DirichletBC(W.sub(0).sub(1), zero, boundaries, 3) #upperlower
bc4 = DirichletBC(W.sub(0).sub(2), zero, boundaries, 4) #frontback
# Boundary condition for velocity on collector boundary
bc5 = DirichletBC(W.sub(0), collector_bc, boundaries, 5) #collector
# Collect boundary conditions
bcs = [bc3, bc4, bc5]
# Define variational problem
v, q = TestFunctions(W)
w = Function(W)
u, p = split(w)
F = (eta*inner(grad(u), grad(v)) - div(v)*p + q*div(u))*dx + \
p1*inner(n,v)*ds(1)
M = inner(u,n)*ds(2)
# Compute solution
solve(F == 0, w, bcs, tol=1.0e-1, M=M)
Why I have that error?
I also try to use
problem = LinearVariationalProblem (a, L, w, bcs)
solver = AdaptiveLinearVariationalSolver (problem, M)
solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "mumps"
tol = 1.0e-1
solver.solve(tol)
but without success.
Another question, can I use boundary definitions as in coarse problem instead a manual definition with adaptivity?
Here I link code files and meshes.
Thank you in advance.
