Merge lp:~pefarrell/dolfin/periodic-malloc into lp:~fenics-core/dolfin/trunk

[Patrick Farrell]

Implement Anders' temporary fix for nonlinear periodic boundary conditions crashing, as discussed in
https://answers.launchpad.net/dolfin/+question/212407 .

This should only live until periodic boundary conditions are "properly fixed" so that the correct sparsity is generated. In the meantime, this works (but is inefficient).

I also added a snippet to the PeriodicBC.py unit test which fails on the current trunk, but which passes on this branch.

[Garth Wells (garth-wells)]

Looks like this breaks with the Epetra backed again. It should be fixed properly.

[Patrick Farrell (pefarrell)]

Are you referring to https://bugs.launchpad.net/dolfin/+bug/1056106 ?

[Anders Logg (logg)]

On Thu, Nov 01, 2012 at 03:39:22PM -0000, Garth Wells wrote:
> Looks like this breaks with the Epetra backed again. It should be fixed properly.

How will this patch affect the Epetra backend? The patch only
introduces a fix that makes PeriodicBC run again with PETSc.

--
Anders

[Garth Wells (garth-wells)]

On Thu, Nov 1, 2012 at 5:52 PM, Anders Logg <email address hidden> wrote:
> On Thu, Nov 01, 2012 at 03:39:22PM -0000, Garth Wells wrote:
>> Looks like this breaks with the Epetra backed again. It should be fixed properly.
>
> How will this patch affect the Epetra backend? The patch only
> introduces a fix that makes PeriodicBC run again with PETSc.
>

That's not the case. The added test in this change in

   test/unit/fem/python/PeriodicBC.py

does not work with Epetra for the same reason the code did not work
with PETSc, and which is why a proper fix is needed.

Garth

> --
> Anders
>
> --
> https://code.launchpad.net/~pefarrell/dolfin/periodic-malloc/+merge/132257
> Your team DOLFIN Core Team is subscribed to branch lp:dolfin.

[Patrick Farrell (pefarrell)]

Garth: can you sketch what a "proper fix" might entail?

[Garth Wells (garth-wells)]

On Thu, Nov 1, 2012 at 9:55 PM, Patrick Farrell
<email address hidden> wrote:
> Garth: can you sketch what a "proper fix" might entail?

Creating a periodic FunctionSpace. This would involve passing the
periodic boundaries though the FunctionSpace constructor and to the
dof map builder. The dof map would just have to make sure that the
'slave' dofs map to the master dofs (i.e. there will be no 'master'
and 'slave' in practice). With this in place, the sparsity pattern and
all the rest will work without needing to be aware of the periodic
boundaries, any symmetry will be preserved and preconditioners won't
get screwed up.

Garth

> --
> https://code.launchpad.net/~pefarrell/dolfin/periodic-malloc/+merge/132257
> Your team DOLFIN Core Team is subscribed to branch lp:dolfin.

[Anders Logg (logg)]

On Thu, Nov 01, 2012 at 06:22:18PM -0000, Garth Wells wrote:
> On Thu, Nov 1, 2012 at 5:52 PM, Anders Logg <email address hidden> wrote:
> > On Thu, Nov 01, 2012 at 03:39:22PM -0000, Garth Wells wrote:
> >> Looks like this breaks with the Epetra backed again. It should be fixed properly.
> >
> > How will this patch affect the Epetra backend? The patch only
> > introduces a fix that makes PeriodicBC run again with PETSc.
> >
>
> That's not the case. The added test in this change in
>
> test/unit/fem/python/PeriodicBC.py
>
> does not work with Epetra for the same reason the code did not work
> with PETSc, and which is why a proper fix is needed.

Fixed now.

--
Anders

=== modified file 'dolfin/fem/PeriodicBC.cpp'
--- dolfin/fem/PeriodicBC.cpp	2012-08-20 07:18:33 +0000
+++ dolfin/fem/PeriodicBC.cpp	2012-10-31 08:38:23 +0000
@@ -32,6 +32,7 @@
 #include <dolfin/function/FunctionSpace.h>
 #include <dolfin/log/log.h>
 #include <dolfin/la/GenericMatrix.h>
+#include <dolfin/la/PETScMatrix.h>
 #include <dolfin/la/GenericVector.h>
 #include <dolfin/mesh/Cell.h>
 #include <dolfin/mesh/Facet.h>
@@ -214,6 +215,20 @@
 void PeriodicBC::apply(GenericMatrix* A, GenericVector* b,
                        const GenericVector* x) const
 {
+
+  // The current handling of periodic boundary conditions adds entries in the matrix that
+  // were not accounted for in the original sparsity pattern. This causes PETSc to crash,
+  // since by default it will not let you do this (as it requires extra memory allocation
+  // and is very inefficient). However, until that's fixed, this hack should stay in, so
+  // that we can run these problems at all.
+#ifdef HAS_PETSC
+  if (A && has_type<PETScMatrix>(*A))
+  {
+    PETScMatrix& petsc_A = A->down_cast<PETScMatrix>();
+    MatSetOption(*petsc_A.mat(), MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
+  }
+#endif
+
   if (MPI::num_processes() > 1)
   {
     parallel_apply(A, b, x);

=== modified file 'test/unit/fem/python/PeriodicBC.py'
--- test/unit/fem/python/PeriodicBC.py	2012-08-18 21:37:40 +0000
+++ test/unit/fem/python/PeriodicBC.py	2012-10-31 08:38:23 +0000
@@ -96,7 +96,7 @@
         bc1 = PeriodicBC(V, PeriodicBoundary2())
         bcs = [bc0, bc1]
 
-        # Define variational problem
+        # Define variational problem, linear formulation
         u, v = TrialFunction(V), TestFunction(V)
         f = Expression("sin(x[0])", degree=2)
         a = dot(grad(u), grad(v))*dx
@@ -108,6 +108,18 @@
 
         self.assertAlmostEqual(u.vector().norm("l2"), 0.3567245204026249, 10)
 
+        # Define variational problem, nonlinear formulation
+        u, v = Function(V), TestFunction(V)
+        f = Expression("sin(x[0])", degree=2)
+        a = dot(grad(u), grad(v))*dx
+        L = f*v*dx
+        F = a - L
+
+        # Compute solution
+        solve(F == 0, u, bcs)
+
+        self.assertAlmostEqual(u.vector().norm("l2"), 0.3567245204026249, 10)
+
 if __name__ == "__main__":
     print ""
     print "Testing Dirichlet boundary conditions"