As a follow-up of 'Broken PETSc wrappers?' thread on this list, can
anyone reproduce incorrect (orders of magnitude) norm using
superlu_dist on following example? Both in serial and parallel. Thanks,
Jan
=============================================================
from dolfin import *
# Create mesh
mesh = UnitSquareMesh(32, 32)
# Define function spaces and mixed (product) space
BDM = FunctionSpace(mesh, "BDM", 1)
DG = FunctionSpace(mesh, "DG", 0)
W = BDM * DG
# Define trial and test functions
(sigma, u) = TrialFunctions(W)
(tau, v) = TestFunctions(W)
# Define source function
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
# Define variational form
a = (dot(sigma, tau) + div(tau)*u + div(sigma)*v)*dx
L = - f*v*dx
# Define function G such that G \cdot n = g
class BoundarySource(Expression):
def __init__(self, mesh):
self.mesh = mesh
def eval_cell(self, values, x, ufc_cell):
cell = Cell(self.mesh, ufc_cell.index)
n = cell.normal(ufc_cell.local_facet)
g = sin(5*x[0])
values[0] = g*n[0]
values[1] = g*n[1]
def value_shape(self):
return (2,)
G = BoundarySource(mesh)
# Define essential boundary
def boundary(x):
return x[1] < DOLFIN_EPS or x[1] > 1.0 - DOLFIN_EPS
bc = DirichletBC(W.sub(0), G, boundary)
# Compute solution
w = Function(W)
solve(a == L, w, bc, solver_parameters={'linear_solver': 'mumps'})
norm_mumps = norm(w)
w.vector().zero()
solve(a == L, w, bc, solver_parameters={'linear_solver': 'superlu_dist'})
norm_superlu_dist = norm(w)
print 'L2 norm mumps ', norm_mumps
print 'L2 norm superlu_dist', norm_superlu_dist
=============================================================
_______________________________________________
fenics-support mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics-support
