On 09/06/2013 04:17 PM, Jan Blechta wrote:
On Fri, 6 Sep 2013 15:59:19 +0200
"Garth N. Wells" <[email protected]> wrote:
On 6 September 2013 12:45, Jan Blechta <[email protected]>
wrote:
It seems that for P1*R function space DofMapBuilder::reorder_local
scales quadratically with number of DOFs. Is it inevitable?
Yes. The 'Real' dof is coupled to all dofs, destroying graph sparsity.
The solutions are (a) extend DofMapBuilder to re-order just sub-dofmap
or (b) do not use Real for anything other than small problems (for
nullspace, they can be handled by an iterative solve). It doesn't play
well with linear solvers, it's bad for assembly and its bad in
parallel.
Ok, I see. Is a result of the following code somewhat similar to
approach (a)?
=========================================
# Build reordered subspaces
parameters['reorder_dofs_serial'] = True
P1 = FunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh, 'R', 0)
# Build mixed space, do not reorder
parameters['reorder_dofs_serial'] = False
V = P1*R
=========================================
This solves the quadratic scaling of DofMapBuilder::reorder_local;
however there must be something else that is quadratic during the
computation of the matrix sparsity pattern (see example belows with
poisson eq.):
#========================================
from dolfin import *
import time
import numpy as np
import matplotlib.pyplot as plt
set_log_level(DEBUG)
parameters['reorder_dofs_serial'] = False
for n in range(25, 33):
mesh = UnitCubeMesh(n, n, n)
#mesh = UnitSquareMesh(10*n, 10*n)
P1 = FunctionSpace(mesh, 'CG', 1)
R = FunctionSpace(mesh, 'R', 0)
# Measure time for creation of P1*R space
t = time.time()
V = P1*R
t = time.time() - t
PP = P1*P1
(uu, ul) = TrialFunctions(V)
(vu, vl) = TestFunctions(V)
(pu, pp) = TrialFunctions(PP)
(pv, pq) = TestFunctions(PP)
u = TrialFunction(P1)
v = TestFunction(P1)
am = inner(grad(uu), grad(vu))*dx
a = inner(grad(pu), grad(pv))*dx
af = inner(grad(u), grad(v))*dx
#mixed cg-real
t1 = time.time()
AM = assemble(am)
t1 = time.time() - t1
#mixed cg-cg
t2 = time.time()
A = assemble(a)
t2 = time.time() - t2
#cg
t3 = time.time()
AF = assemble(af)
t3 = time.time() - t3
# Report
print 'DOFs:', V.dim(),\
'T/D:', t/V.dim(),\
'T1/D:', t1/P1.dim(),\
'T2/D:', t2/P1.dim(),\
'T3/D:', t3/P1.dim()
plt.loglog(V.dim(), t1, 'o')
plt.loglog(V.dim(), t2, 'x')
plt.loglog(V.dim(), t3, '*')
#
# Plot some quadratic function for comparison
xl = plt.xlim()
yl = plt.ylim()
x = np.linspace(*xl)
y = map(lambda x:yl[0]/(xl[0]*xl[0])*x*x, x)
plt.loglog(x, y)
#
# # Show plot
plt.show()
#=============================================
Marco
_______________________________________________
fenics mailing list
[email protected]
http://fenicsproject.org/mailman/listinfo/fenics
