cbc.block team mailing list archive

[Claas Abert <question220333@xxxxxxxxxxxxxxxxxxxxx>]

Question #220333 on CBC.Block changed:
https://answers.launchpad.net/cbc.block/+question/220333

    Status: Answered => Open

Claas Abert is still having a problem:
Thanks Joachim Haga and Kent-Andre Mardal for the quick response.
I'm using the dolfin-dev rev 7250 (from the beginning of 2013) and the most recent dev-version of cbc.block (rev 121).
Here is a (not very minimal) example showing the problem.

The problem itself is a saddle point problem and I know that I should
probably deal with this block structure properly using CBC.Block.
However this is just the starting point for a bigger block system and I
wanted to keep things simple in the beginning...

from dolfin import *
from block  import *
from block.iterative import *
from block.algebraic.trilinos import *

# Set backend to Epetra
parameters['linear_algebra_backend'] = 'Epetra'

# Mesh and Function Spaces
mesh = BoxMesh(-50, -50, -5, 50, 50, 5, 10, 10, 1)
VV   = VectorFunctionSpace(mesh, "CG", 1)
VS   = FunctionSpace(mesh, "CG", 1)
V    = VV * VS

# Test and Trial Functions
(v, sigma) = TrialFunctions(V)
(w2, w3)   = TestFunctions(V)

# Define some constants and coeffecients
arg    = "sqrt((3.141592*(x[0]/1e1))*(3.141592*(x[0]/1e1)))"
m_expr = Expression(("cos(%s)" % arg, "sin(%s)" % arg, "0.0"))
m      = interpolate(m_expr, VV)

f_ex  = -5.71823816178e+12
alpha = 0.1
dt    = 1e-12

# Define weak forms
a  = alpha * dot(v, w2) * dx + dot(cross(m, v), w2) * dx
a += - 0.5 * dt * f_ex * Dx(v[i],j) * Dx(w2[i],j) * dx
a += sigma * inner(m, w2) * dx
a += inner(m, v) * w3 * dx

L = f_ex * Dx(m[i],j) * Dx(w2[i],j) * dx

# Assemble
A = assemble(a)
b = assemble(L)

# Solve with Trilinos solvers
x = Function(V)
solve(A, x.vector(), b, "gmres", "ilu")
solve(A, x.vector(), b, "bicgstab", "ilu")

# Solve with CBC.Block solvers
Ap   = ILU(A)
Ainv = LGMRES(A, precond=Ap)
#Ainv = BiCGStab(A, precond=Ap)
x    = Ainv * b

-- 
You received this question notification because you are a member of
cbc.block maintainers, which is an answer contact for CBC.Block.