[cig-commits] commit 1999 by buerg to /var/svn/dealii/aspect

[dealii.demon at gmail.com]

Revision 1999

Make compressible case work.

U   trunk/aspire/include/aspect/material_model/interface.h
U   trunk/aspire/source/material_model/interface.cc
U   trunk/aspire/source/simulator/assembly.cc


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

Diff:
Modified: trunk/aspire/include/aspect/material_model/interface.h
===================================================================
--- trunk/aspire/include/aspect/material_model/interface.h	2013-10-30 16:12:59 UTC (rev 1998)
+++ trunk/aspire/include/aspect/material_model/interface.h	2013-10-31 15:29:09 UTC (rev 1999)
@@ -570,6 +570,11 @@
              * Compressibility at the given positions.
              */
             std::vector<double> compressibilities;
+            /**
+             * Changes in density with respect to temperature and
+             * composition at the given points.
+             */
+            std::vector<std::vector<double> > density_changes;
           };
 
         /**

Modified: trunk/aspire/source/material_model/interface.cc
===================================================================
--- trunk/aspire/source/material_model/interface.cc	2013-10-30 16:12:59 UTC (rev 1998)
+++ trunk/aspire/source/material_model/interface.cc	2013-10-31 15:29:09 UTC (rev 1999)
@@ -338,9 +338,13 @@
       specific_heat.resize (n_points);
       thermal_conductivities.resize (n_points);
       compositional_sources.resize (n_points);
+      density_changes.resize (n_points);
       
       for (unsigned int q = 0; q < n_points; ++q)
+      {
         compositional_sources[q].resize (n_comp);
+        density_changes[q].resize (n_comp + 1);
+      }
       
       thermal_sources.resize (n_points);
       compressibilities.resize (n_points);

Modified: trunk/aspire/source/simulator/assembly.cc
===================================================================
--- trunk/aspire/source/simulator/assembly.cc	2013-10-30 16:12:59 UTC (rev 1998)
+++ trunk/aspire/source/simulator/assembly.cc	2013-10-31 15:29:09 UTC (rev 1999)
@@ -972,13 +972,12 @@
 
         const double eta = scratch.material_model_outputs.viscosities[q];
         const double density = scratch.material_model_outputs.densities[q];
-        const Tensor<1,dim> gravity = gravity_model->gravity_vector (scratch.fe_values_velocity.quadrature_point (q));
-        double right_hand_side = (scratch.velocity_values[q] * scratch.temperature_gradients[q]) / scratch.temperature_values[q];
+        double right_hand_side = scratch.material_model_outputs.density_changes[q][0] * (scratch.velocity_values[q] * scratch.temperature_gradients[q]);
 
-//        for (unsigned int n = 0; n < parameters.n_compositional_fields; ++n)
-//          right_hand_side += (scratch.velocity_values[q] * scratch.composition_gradients[n][q]) / scratch.composition_values[n][q];
+        for (unsigned int n = 0; n < parameters.n_compositional_fields; ++n)
+          right_hand_side += scratch.material_model_outputs.density_changes[q][n + 1] * (scratch.velocity_values[q] * scratch.composition_gradients[n][q]);
         
-        right_hand_side *= scratch.fe_values_velocity.JxW (q) * (axisymmetric_terms ? scratch.fe_values_velocity.quadrature_point (q) (0) : 1);
+        right_hand_side *= -1.0 * scratch.fe_values_velocity.JxW (q) * (axisymmetric_terms ? scratch.fe_values_velocity.quadrature_point (q) (0) : 1) / density;
         
         for (unsigned int i=0; i<dofs_per_cell; ++i)
           {
@@ -1169,7 +1168,7 @@
         const double c_P = scratch.material_model_outputs.specific_heat[q];
         const double lambda = scratch.material_model_outputs.thermal_conductivities[q];
         const double old_value = scratch.old_temperature_values[q];
-        const double rho = 0.0*scratch.material_model_outputs.densities[q];
+        const double rho = scratch.material_model_outputs.densities[q];
         const Tensor<1, dim> current_u = scratch.current_velocity_values[q];
 
         for (unsigned int i = 0; i < dofs_per_cell; ++i)
@@ -1237,7 +1236,7 @@
             scratch.phi_field[i]      = scratch.fe_values_advection.shape_value (i, q);
           }
         
-        const double rho = 0.0*scratch.material_model_outputs.densities[q];
+        const double rho = scratch.material_model_outputs.densities[q];
         const Tensor<1, dim> current_u = scratch.current_velocity_values[q];
         
         for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)