Error with GenericLinearOperator and argument type


284
views
-1
6 months ago by
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)


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
6 months ago by
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).
Please login to add an answer/comment or follow this question.

Similar posts:
Search »