[cig-commits] [commit] master: Only remove the nullspace in a direction if the user specifies it (60908b9)

[cig_noreply at geodynamics.org]

Repository : https://github.com/geodynamics/aspect

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/2df080f4905a6be43fd1729ea0d6a7a956cce070...a7135c1f7697d39efff2f47a79ca1e1395cff504

>---------------------------------------------------------------

commit 60908b9de6d00e688df90711bc2abd2b6e0cda21
Author: ian-r-rose <ian.r.rose at gmail.com>
Date:   Wed May 21 12:33:16 2014 -0500

    Only remove the nullspace in a direction if the user specifies it


>---------------------------------------------------------------

60908b9de6d00e688df90711bc2abd2b6e0cda21
 include/aspect/simulator.h     | 12 ++++++++----
 source/simulator/nullspace.cc  | 25 ++++++++++++++++++-------
 source/simulator/parameters.cc | 26 ++++++++++++++++++++------
 3 files changed, 46 insertions(+), 17 deletions(-)

diff --git a/include/aspect/simulator.h b/include/aspect/simulator.h
index 82ea21a..f44e9cd 100644
--- a/include/aspect/simulator.h
+++ b/include/aspect/simulator.h
@@ -113,10 +113,14 @@ namespace aspect
         enum Kind
         {
           none = 0,
-          net_rotation = 0x1,
-          net_translation = 0x2,
-          angular_momentum = 0x4,
-          translational_momentum = 0x8
+          net_translation_x = 0x2,
+          net_translation_y = 0x4,
+          net_translation_z = 0x8,
+          linear_momentum_x = 0x10,
+          linear_momentum_y = 0x20,
+          linear_momentum_z = 0x40,
+          net_rotation      = 0x80,
+          angular_momentum  = 0x100
         };
       };
 
diff --git a/source/simulator/nullspace.cc b/source/simulator/nullspace.cc
index 2091c35..440293a 100644
--- a/source/simulator/nullspace.cc
+++ b/source/simulator/nullspace.cc
@@ -100,7 +100,9 @@ namespace aspect
   template <int dim>
   void Simulator<dim>::setup_nullspace_removal()
   {
-    if (parameters.nullspace_removal & NullspaceRemoval::translational_momentum)
+    if (parameters.nullspace_removal & NullspaceRemoval::linear_momentum_x || 
+        parameters.nullspace_removal & NullspaceRemoval::linear_momentum_y || 
+        parameters.nullspace_removal & NullspaceRemoval::linear_momentum_z )
       AssertThrow(false, ExcNotImplemented());
 
     std::vector<std_cxx1x::shared_ptr<TensorFunction<1,dim> > > funcs;
@@ -114,11 +116,18 @@ namespace aspect
             funcs.push_back(std_cxx1x::shared_ptr<TensorFunction<1,dim> >(new internal::Rotation<dim>(a)));
       }
 
-    if (parameters.nullspace_removal & NullspaceRemoval::net_translation)
-      {
-        for (unsigned int a=0; a<dim; ++a)
-          funcs.push_back(std_cxx1x::shared_ptr<TensorFunction<1,dim> >(new internal::Translation<dim>(a)));
-      }
+    if (parameters.nullspace_removal & NullspaceRemoval::net_translation_x)
+      funcs.push_back(std_cxx1x::shared_ptr<TensorFunction<1,dim> >(new internal::Translation<dim>(0))); //x dir
+
+    if (parameters.nullspace_removal & NullspaceRemoval::net_translation_y)
+      funcs.push_back(std_cxx1x::shared_ptr<TensorFunction<1,dim> >(new internal::Translation<dim>(1))); //y dir
+
+    if (parameters.nullspace_removal & NullspaceRemoval::net_translation_z)
+    { 
+      //Only do z direction if dim == 3
+      AssertThrow( dim == 3, ExcMessage("Can't remove z translational mode in 2 dimensions"));
+      funcs.push_back(std_cxx1x::shared_ptr<TensorFunction<1,dim> >(new internal::Translation<dim>(2))); //z dir
+    }
 
     if (funcs.size()>0)
       {
@@ -149,7 +158,9 @@ namespace aspect
                                         LinearAlgebra::BlockVector &tmp_distributed_stokes)
   {
     if (parameters.nullspace_removal & NullspaceRemoval::net_rotation ||
-        parameters.nullspace_removal & NullspaceRemoval::net_translation)
+        parameters.nullspace_removal & NullspaceRemoval::net_translation_x || 
+        parameters.nullspace_removal & NullspaceRemoval::net_translation_y ||
+        parameters.nullspace_removal & NullspaceRemoval::net_translation_z )
       {
         Assert(introspection.block_indices.velocities != introspection.block_indices.pressure,
             ExcNotImplemented());
diff --git a/source/simulator/parameters.cc b/source/simulator/parameters.cc
index 2ec6db5..c03681a 100644
--- a/source/simulator/parameters.cc
+++ b/source/simulator/parameters.cc
@@ -362,7 +362,9 @@ namespace aspect
                          "current parameter section.");
 
       prm.declare_entry ("Remove nullspace", "",
-                         Patterns::MultipleSelection("net rotation|net translation|angular momentum|translational momentum"),
+                         Patterns::MultipleSelection("net rotation|angular momentum|"
+                                                     "net x translation|net y translation|net z translation|"
+                                                     "linear x momentum|linear y momentum|linear z momnetum"),
                          "A selection of operations to remove certain parts of the nullspace from "
                          "the velocity after solving. For some geometries and certain boundary conditions "
                          "the velocity field is not uniquely determined but contains free translations "
@@ -758,15 +760,27 @@ namespace aspect
             if (nullspace_names[i]=="net rotation")
               nullspace_removal = typename NullspaceRemoval::Kind(
                                     nullspace_removal | NullspaceRemoval::net_rotation);
-            else if (nullspace_names[i]=="net translation")
-              nullspace_removal = typename NullspaceRemoval::Kind(
-                                    nullspace_removal | NullspaceRemoval::net_translation);
             else if (nullspace_names[i]=="angular momentum")
               nullspace_removal = typename NullspaceRemoval::Kind(
                                     nullspace_removal | NullspaceRemoval::angular_momentum);
-            else if (nullspace_names[i]=="translational momentum")
+            else if (nullspace_names[i]=="net x translation")
+              nullspace_removal = typename NullspaceRemoval::Kind(
+                                    nullspace_removal | NullspaceRemoval::net_translation_x);
+            else if (nullspace_names[i]=="net y translation")
+              nullspace_removal = typename NullspaceRemoval::Kind(
+                                    nullspace_removal | NullspaceRemoval::net_translation_y);
+            else if (nullspace_names[i]=="net z translation")
+              nullspace_removal = typename NullspaceRemoval::Kind(
+                                    nullspace_removal | NullspaceRemoval::net_translation_z);
+            else if (nullspace_names[i]=="linear x momentum")
+              nullspace_removal = typename       NullspaceRemoval::Kind(
+                                    nullspace_removal | NullspaceRemoval::linear_momentum_x);
+            else if (nullspace_names[i]=="linear y momentum")
+              nullspace_removal = typename       NullspaceRemoval::Kind(
+                                    nullspace_removal | NullspaceRemoval::linear_momentum_y);
+            else if (nullspace_names[i]=="linear z momentum")
               nullspace_removal = typename       NullspaceRemoval::Kind(
-                                    nullspace_removal | NullspaceRemoval::translational_momentum);
+                                    nullspace_removal | NullspaceRemoval::linear_momentum_z);
             else
               AssertThrow(false, ExcInternalError());
           }