On 30 September 2013 11:01, Marco Morandini <[email protected]> wrote: > 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.): >
'Real' spaces are a problem and shouldn't be used for large problems. They lead to dense rows in the matrix. Sparsity pattern construction and sparse matrix insertion are designed/optimised for sparse rows. Garth > > #======================================== > 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
