I want to use Nitsche’s method to solve stationary NS eqaution.However,I encounter a wierd problem.
Here is my code.
from dolfin import *
from dolfin_adjoint import *
import moola
import matplotlib.pyplot as plt
...
import meshio
geometry = meshio.read("glacier_data/mesh.msh")
meshio.write("glacier_data/mesh.xdmf", meshio.Mesh(points=geometry.points[:,:2], cells={"triangle": geometry.cells["triangle"]}))
meshio.write("glacier_data/mf.xdmf", meshio.Mesh(points=geometry.points[:,:2], cells={"line": geometry.cells["line"]},
cell_data={"line": {"name_to_read": geometry.cell_data["line"]["gmsh:physical"]}}))
mesh = Mesh()
with XDMFFile("glacier_data/mesh.xdmf") as infile:
infile.read(mesh)
mvc = MeshValueCollection("size_t", mesh, 1)
with XDMFFile("glacier_data/mf.xdmf") as infile:
infile.read(mvc, "name_to_read")
boundary_markers =MeshFunction("size_t",mesh, mvc)
# Define function spaces (P2-P1)
V = VectorElement("CG", mesh.ufl_cell(), 2)
Q = FiniteElement("CG", mesh.ufl_cell(), 1)
TH = V * Q
W = FunctionSpace(mesh, TH)
# Define test and solution functions
v, q = TestFunctions(W)
s = Function(W,name='State')
u, p = split(s)
f = Function(W.sub(0).collapse(),name="control")
# Set parameter values
nu = Constant(1) #Viscosity coefficient
gamma=Constant(10) #Nitsche penalty parameter
n=FacetNormal(mesh)
h=CellVolume(mesh)
# Define boundary conditions
noslip = DirichletBC(W.sub(0), (0, 0), boundary_markers,2)
inflow = DirichletBC(W.sub(0), f, boundary_markers,0)
outflow = DirichletBC(W.sub(0), f, boundary_markers,1)
bcs=[inflow,outflow,noslip]
ds=Measure('ds',domain=mesh,subdomain_data=boundary_markers)
# Define the variational formulation of the Navier-Stokes equations
a=nu*inner(grad(u),grad(v))*dx+inner(grad(u)*u,v)*dx-nu*inner(grad(u)*n,v)*ds(2)+\
inner(p*n,v)*dx-p*div(v)*dx-q*div(u)*dx-nu*inner(grad(v)*n,u)*ds(0)-nu*inner(grad(v)*n,u)*ds(1)+\
gamma/h*nu*(inner(u,v)*ds(0)+inner(u,v)*ds(1))+\
inner(q*n,u)*ds(0)+inner(q*n,u)*ds(1)
L=-nu*inner(grad(v)*n,f)*ds(0)-nu*inner(grad(v)*n,f)*ds(1)+gamma/h*nu*inner(f,v)*ds(0)+gamma/h*nu*inner(f,v)*ds(1)+\
inner(q*n,f)*ds(0)+inner(q*n,f)*ds(1)
solve(a==L,s,bcs)
#moola optimization
...
Error
/usr/bin/python3.6 /home/rsrg/PycharmProjects/glacier_velocity/start_1.py
Traceback (most recent call last):
File "/home/rsrg/PycharmProjects/glacier_velocity/start_1.py", line 85, in <module>
inner(q*n,u)*ds(0)+inner(q*n,u)*ds(1)
TypeError: unsupported operand type(s) for *: 'Product' and 'Form'
Process finished with exit code 1