[cig-commits] [commit] master: Included the pressure in the calculation for cases that are not completely incompressible (67296e2)

[cig_noreply at geodynamics.org]

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

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/1cea9a9309fea18bec922b55d3e1aed04c0372d2...0801e5573990bb6b2572872a2c084577ffc9fa4a

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

commit 67296e2bdb345f8d268f1dd072ed62ed1ba83eae
Author: anne-glerum <a.c.glerum at uu.nl>
Date:   Thu May 22 18:46:50 2014 +0200

    Included the pressure in the calculation for cases that are not completely incompressible


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

67296e2bdb345f8d268f1dd072ed62ed1ba83eae
 .../postprocess/viscous_dissipation_statistics.cc  | 25 ++++++++++++++++------
 1 file changed, 19 insertions(+), 6 deletions(-)

diff --git a/source/postprocess/viscous_dissipation_statistics.cc b/source/postprocess/viscous_dissipation_statistics.cc
index ad1627e..8c42df6 100644
--- a/source/postprocess/viscous_dissipation_statistics.cc
+++ b/source/postprocess/viscous_dissipation_statistics.cc
@@ -40,7 +40,9 @@ namespace aspect
     ViscousDissipationStatistics<dim>::execute (TableHandler &statistics)
     {
       const QGauss<dim> quadrature_formula (this->get_fe()
-                                            .base_element(0).degree+1);
+                                            .base_element(this->introspection().base_elements.velocities)
+                                            .degree+1);
+
       const unsigned int n_q_points = quadrature_formula.size();
 
       FEValues<dim> fe_values (this->get_mapping(),
@@ -93,11 +95,17 @@ namespace aspect
             // get the viscosity from the material model
             this->get_material_model().evaluate(in, out);
 
-            // integrate the local viscous dissipation $2\mu(\dot{\epsilon}::\dot{\epsilon})$ over the cell
+            // calculate the local viscous dissipation integral
             for (unsigned int q = 0; q < n_q_points; ++q)
               {
-                local_dissipation_integral += out.viscosities[q] * 2.0 * in.strain_rate[q].norm() *
-                                              in.strain_rate[q].norm() * fe_values.JxW(q);
+/*                local_dissipation_integral += (2.0 * out.viscosities[q] * in.strain_rate[q].norm() *
+                                              in.strain_rate[q].norm() 
+                                              - in.pressure[q] * trace(in.strain_rate[q]))
+                                              * fe_values.JxW(q);
+*/                local_dissipation_integral += (2.0 * out.viscosities[q] * in.strain_rate[q].norm() *
+                                                in.strain_rate[q].norm() 
+                                              - in.pressure[q] * trace(in.strain_rate[q]))
+                                              * fe_values.JxW(q);
               }
           }
 
@@ -109,7 +117,8 @@ namespace aspect
         {
           // make sure that the columns filled by the this object
           // all show up with sufficient accuracy and in scientific notation
-          const char *columns[] = { "Total viscous dissipation (J/yr/m)", //TODO does it make sense to output this?
+          //TODO does it make sense to output this?
+          const char *columns[] = { "Total viscous dissipation (J/yr/m)",
                                     "Total viscous dissipation (kJ/yr/m)"
                                   };
 
@@ -168,6 +177,10 @@ namespace aspect
     ASPECT_REGISTER_POSTPROCESSOR(ViscousDissipationStatistics,
                                   "viscous dissipation statistics",
                                   "A postprocessor that computes the viscous dissipation"
-                                  "for the whole domain.")
+                                  "for the whole domain as: "
+                                  "$\\frac{1}{2} \\int_{V} \\mathbf{\\sigma} : \\mathbr{\\dot{\\epsilon}}dV "
+                                  "= $\\frac{1}{2} \\int_{V} (-p\\delta_{ij}+2\\mu\\dot{\\epsilon}_{ij}) "
+                                  "(\\dot{\\epsilon_{ij}})dV$ ."
+                                 )
   }
 }