### Error with GenericLinearOperator and argument type

284

views

-1

Hi,

I'm getting this error in my script, for dolfin::GenericLinearOperator. I've tried using different solvers but I haven't figured it out. I attached a shorter version of my code that gives the same error, I hope somebody can help! :)

File "NCAR_NT1.py", line 269, in <module>

solve(A1, u_.vector(), b1)

File "/scratch/anaconda2/envs/fenicsproject/lib/python2.7/site-packages/dolfin/fem/solving.py", line 310, in solve

return cpp.la_solve(*args)

File "/scratch/anaconda2/envs/fenicsproject/lib/python2.7/site-packages/dolfin/cpp/la.py", line 5274, in la_solve

return _la.la_solve(*args)

TypeError: in method 'la_solve', argument 1 of type 'dolfin::GenericLinearOperator const &'

Aborted (core dumped)

I'm getting this error in my script, for dolfin::GenericLinearOperator. I've tried using different solvers but I haven't figured it out. I attached a shorter version of my code that gives the same error, I hope somebody can help! :)

File "NCAR_NT1.py", line 269, in <module>

solve(A1, u_.vector(), b1)

File "/scratch/anaconda2/envs/fenicsproject/lib/python2.7/site-packages/dolfin/fem/solving.py", line 310, in solve

return cpp.la_solve(*args)

File "/scratch/anaconda2/envs/fenicsproject/lib/python2.7/site-packages/dolfin/cpp/la.py", line 5274, in la_solve

return _la.la_solve(*args)

TypeError: in method 'la_solve', argument 1 of type 'dolfin::GenericLinearOperator const &'

Aborted (core dumped)

```
from __future__ import print_function
from fenics import *
from mshr import *
import numpy as np
## TIME SPECS ##
T = 23.0
dt = 7.5 # For progress bar
k = Constant(dt) # For variational form
# DOMAIN SPECS
Lx = 320
Ly = 320
Lz = 96
rsl = 100
# # PARAMETERS
Coriolis = Constant(0.0000729) # Coriolis magnitude
rho0 = Constant(1000.0) # reference density [kg/m^3] not for buoyancy
g= Constant(9.81) #for buoyancy implementation
B_T = Constant(2.0 * 10 **(-4)) # H2O thermal expansion coefficient - buoyancy
theta0 = Constant(290.16) # reference temperature [K]
u_z0 = 0.022 # Surface wind stress boundary condition
### MESH ###
# Domain from Hamlington paper, 100^3 resolution
mesh = BoxMesh(Point(0,0,0), Point(Lx,Ly,-96), 100,100,100)
x = SpatialCoordinate(mesh)
# FUNCTION SPACE DEFINITIONS #
V = VectorFunctionSpace(mesh,'CG', 1, dim=3) # velocity
Q = FunctionSpace(mesh,'CG', 1) # pressure
R = FunctionSpace(mesh,'CG', 1) # temperature
#### FUNCTION DEFINITION FOR VARIABLES ####
# VELOCITY #
u = TrialFunction(V)
u = Function(V)
v = TestFunction(V)
# PRESSURE #
p = TrialFunction(Q)
p = Function(Q)
q = TestFunction(Q)
# TEMPERATURE #
theta = TrialFunction(R)
theta = Function(R)
r = TestFunction(R)
u_n = Function(V) # previous velocity solution
u_ = Function(V) # current velocity (for velocity correction in Poisson)
p_n = Function(Q) # previous pressure solution
theta_n = Function(R) # previous temperature solution
#### BOUNDARY CONDITIONS ####
# Surfaces, for boundary conditions
top = 'near(x[2],0)'
bottom = 'near(x[2],-96)'
# Boundary conditions on the domain
bc_top = DirichletBC(V.sub(2) , Constant(0) , top) # No vertical velocity @ top
bc_bottom = DirichletBC(V.sub(2) , Constant(0), bottom) # NVV @ bottom
bc_wind = DirichletBC(V.sub(0) , Constant(u_z0), top) # Wind sirface forcing
bcu = [bc_top, bc_bottom, bc_wind] # gather up BCs for function space constrained_domain
# (f x u) + b(z)) term in Craik Leibovich equations
def GeosB(u,theta):
return as_vector([Coriolis*u[1]] + [-Coriolis*u[0]] + [g*(1+B_T*theta0-B_T*theta)] )
# Initialize velocity field
u0 = Constant((DOLFIN_EPS,DOLFIN_EPS,DOLFIN_EPS))
u_n = interpolate(u0,V)
# # VARIATIONAL FORM DEFINITION
# Velocity equation (without grad(P) because of Poisson)
a1 = dot(u_,v)*dx
L1 = dot(u_n,v)*dx + k*( - dot(dot(u_n, nabla_grad(u_n)),v)*dx - dot(GeosB(u_n,theta_n),v)*dx) #- inner(tij(u_n),nabla_grad(v))*dx)
# From Poisson pressure correction
a2 = inner(grad(p), grad(q))*dx
L2 = - (1/k)*div(u_)*q*dx
a3 = inner(u, v)*dx
L3 = inner(u_, v)*dx - k*inner(grad(p_n), v)*dx
# Potential Temperature/Buoyancy
a4 = theta*r*dx
L4 = theta_n*r*dx - k*(dot(u_n , grad(theta_n))*r*dx)
# Assemble forms
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
A4 = assemble(a4)
# Time stepping
t=0
while t < T + DOLFIN_EPS:
# Step 1: Tentative velocity step
begin("Computing tentative velocity")
b1 = assemble(L1)
[bc.apply(b1) for bc in bcu]
solve(A1, u_.vector(), b1)
end()
# Pressure correction
begin("Computing pressure correction")
solve(a2==L2, p, solver_parameters={"linear_solver":"cg"},form_compiler_parameters={"optimize":True})
end()
#Velocity correction
begin("Computing velocity correction")
b3 = assemble(L3)
[bc.apply(b3) for bc in bcu]
solve(A3, u.vector(), b3, "bicgstab","default")
end()
# Temperature
begin("Computing temperature")
solve(a4 == L4, theta, solver_parameters={"linear_solver":"cg"},form_compiler_parameters={"optimize":True})
end()
# Update previous solution
u_n.assign(u)
p_n.assign(p)
theta_n.assign(theta)
t += dt
```

Community: FEniCS Project

### 1 Answer

0

Hi,

on the Left Hand Side of your equations, should only exist trial functions. You have declared u_ which is not a trial function. See step 3 and make the appropriate changes. Also, use another name for your solution vectors (see lines 38 and 41 should be removed).

on the Left Hand Side of your equations, should only exist trial functions. You have declared u_ which is not a trial function. See step 3 and make the appropriate changes. Also, use another name for your solution vectors (see lines 38 and 41 should be removed).

Please login to add an answer/comment or follow this question.