[cig-commits] commit 2414 by heister to /var/svn/dealii/aspect

[dealii.demon at gmail.com]

Revision 2414

implement angular momentum in 2d (thanks Ian)

U   trunk/aspect/include/aspect/simulator.h
U   trunk/aspect/source/simulator/nullspace.cc
U   trunk/aspect/source/simulator/solver.cc


http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=2414&peg=2414

Diff:
Modified: trunk/aspect/include/aspect/simulator.h
===================================================================
--- trunk/aspect/include/aspect/simulator.h	2014-04-07 00:50:22 UTC (rev 2413)
+++ trunk/aspect/include/aspect/simulator.h	2014-04-07 00:55:17 UTC (rev 2414)
@@ -930,11 +930,26 @@
       void setup_nullspace_removal();
 
       /**
-       * Eliminate the nullspace of the velocity in the given vector @p vector.
+       * Eliminate the nullspace of the velocity in the given vector. Both vectors
+       * are expected to contain the up to date data.
+       *
+       * @param relevant_dst locally relevant vector for the whole FE, will be filled at the end.
+       * @param tmp_distributed_stokes only contains velocity and pressure.
        */
-      void remove_nullspace(LinearAlgebra::BlockVector &vector);
+      void remove_nullspace(LinearAlgebra::BlockVector &relevant_dst,
+          LinearAlgebra::BlockVector &tmp_distributed_stokes);
 
       /**
+       * Remove the angular momentum of the given vector
+       */
+      void remove_net_angular_momentum( LinearAlgebra::BlockVector &relevant_dst, LinearAlgebra::BlockVector &tmp_distributed_stokes);
+
+      /**
+       * Remove the linear momentum of the given vector
+       */
+      void remove_net_linear_momentum( LinearAlgebra::BlockVector &relevant_dst, LinearAlgebra::BlockVector &tmp_distributed_stokes);
+
+      /**
        * Compute the maximal velocity throughout the domain. This is needed
        * to compute the size of the time step.
        *
@@ -1147,7 +1162,6 @@
       /**
       * @}
       */
-
       /**
        * @name Variables that describe the time discretization
        * @{

Modified: trunk/aspect/source/simulator/nullspace.cc
===================================================================
--- trunk/aspect/source/simulator/nullspace.cc	2014-04-07 00:50:22 UTC (rev 2413)
+++ trunk/aspect/source/simulator/nullspace.cc	2014-04-07 00:55:17 UTC (rev 2414)
@@ -35,6 +35,9 @@
 #include <deal.II/lac/pointer_matrix.h>
 #include <deal.II/base/tensor_function.h>
 
+#include <deal.II/base/quadrature_lib.h>
+#include <deal.II/fe/fe_values.h>
+
 namespace aspect
 {
   namespace internal
@@ -68,8 +71,6 @@
   template <int dim>
   void Simulator<dim>::setup_nullspace_removal()
   {
-    if (parameters.nullspace_removal & NullspaceRemoval::angular_momentum)
-        AssertThrow(false, ExcNotImplemented());
     if (parameters.nullspace_removal & NullspaceRemoval::translational_momentum)
         AssertThrow(false, ExcNotImplemented());
 
@@ -110,7 +111,7 @@
   }
 
   template <int dim>
-  void Simulator<dim>::remove_nullspace(LinearAlgebra::BlockVector &vector)
+  void Simulator<dim>::remove_nullspace(LinearAlgebra::BlockVector &relevant_dst, LinearAlgebra::BlockVector &tmp_distributed_stokes)
   {
     if (parameters.nullspace_removal & NullspaceRemoval::net_rotation ||
         parameters.nullspace_removal & NullspaceRemoval::net_translation)
@@ -118,15 +119,104 @@
         for(unsigned int i=0; i<net_rotations_translations.size(); ++i)
         {
            double power = net_rotations_translations[i]
-                          * vector.block(introspection.block_indices.velocities);
-           vector.block(introspection.block_indices.velocities).sadd(1.0,
+                          * tmp_distributed_stokes.block(introspection.block_indices.velocities);
+           tmp_distributed_stokes.block(introspection.block_indices.velocities).sadd(1.0,
                      -1.0*power,
                      net_rotations_translations[i]);
-           pcout << "removing NULLSPACE " << i << " power: " << power << std::endl;
         }
+        relevant_dst.block(0) = tmp_distributed_stokes.block(0);
       }
+    if (parameters.nullspace_removal & NullspaceRemoval::angular_momentum)
+      {
+        remove_net_angular_momentum( relevant_dst, tmp_distributed_stokes);
+      }
   }
 
+  template <>
+  void
+  Simulator<3>::remove_net_angular_momentum( LinearAlgebra::BlockVector &relevant_dst, LinearAlgebra::BlockVector &tmp_distributed_stokes )
+  {
+    AssertThrow(false, ExcNotImplemented());
+  }
+
+  template <>
+  void
+  Simulator<2>::remove_net_angular_momentum( LinearAlgebra::BlockVector &relevant_dst, LinearAlgebra::BlockVector &tmp_distributed_stokes )
+  {
+
+    // compute and remove angular momentum from velocity field, by computing
+    // int rho V