[cig-commits] commit 1876 by dannberg to /var/svn/dealii/aspect

[dealii.demon at gmail.com]

Revision 1876

fix some bugs

U   branches/j-dannberg/include/aspect/simulator_access.h
U   branches/j-dannberg/source/material_model/composition_reaction.cc
U   branches/j-dannberg/source/material_model/latent_heat.cc
U   branches/j-dannberg/source/simulator/simulator_access.cc


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

Diff:
Modified: branches/j-dannberg/include/aspect/simulator_access.h
===================================================================
--- branches/j-dannberg/include/aspect/simulator_access.h	2013-09-11 14:21:08 UTC (rev 1875)
+++ branches/j-dannberg/include/aspect/simulator_access.h	2013-09-11 14:28:30 UTC (rev 1876)
@@ -174,6 +174,12 @@
       include_adiabatic_heating () const;
 
       /**
+       * Return whether we use the latent heat term.
+       */
+      bool
+      include_latent_heat () const;
+
+      /**
        * Return the adiabatic surface temperature.
        */
       double

Modified: branches/j-dannberg/source/material_model/composition_reaction.cc
===================================================================
--- branches/j-dannberg/source/material_model/composition_reaction.cc	2013-09-11 14:21:08 UTC (rev 1875)
+++ branches/j-dannberg/source/material_model/composition_reaction.cc	2013-09-11 14:28:30 UTC (rev 1876)
@@ -56,7 +56,13 @@
 
         // last, calculate the percentage of material that has undergone the transition
         // (also in dependence of the phase transition width - this is an input parameter)
-        return 0.5*(1.0 + std::tanh(pressure_deviation / pressure_width));
+    	double phase_func;
+    	// use delta function for width = 0
+    	if(transition_widths[phase]==0)
+    	  (pressure_deviation > 0) ? phase_func = 1 : phase_func = 0;
+    	else
+    	  phase_func = 0.5*(1.0 + std::tanh(pressure_deviation / pressure_width));
+        return phase_func;
       }
       // if we do not have the adiabatic conditions, we have to use the depth instead
       // this is less precise, because we do not have the exact pressure gradient, instead we use pressure/depth
@@ -64,9 +70,21 @@
       else
       {
     	double depth = this->get_geometry_model().depth(position);
-    	double depth_deviation = depth - transition_depths[phase]
-    	                         - transition_slopes[phase] * (depth / pressure) * (temperature - transition_temperatures[phase]);
-    	return 0.5*(1.0 + std::tanh(depth_deviation / transition_widths[phase]));
+    	double depth_deviation = (pressure > 0
+    			                  ?
+    			                  depth - transition_depths[phase]
+    	                          - transition_slopes[phase] * (depth / pressure) * (temperature - transition_temperatures[phase])
+    	                          :
+    			                  depth - transition_depths[phase]
+    	                          - transition_slopes[phase] / (this->get_gravity_model().gravity_vector(position).norm() * reference_rho)
+    	                          * (temperature - transition_temperatures[phase]));
+    	double phase_func;
+    	// use delta function for width = 0
+    	if(transition_widths[phase]==0)
+    	  (depth_deviation > 0) ? phase_func = 1 : phase_func = 0;
+    	else
+    	  phase_func = 0.5*(1.0 + std::tanh(depth_deviation / transition_widths[phase]));
+    	return phase_func;
       }
   	}
 
@@ -96,8 +114,11 @@
                                 - transition_slopes[phase] * (temperature - transition_temperatures[phase]);
 
       // last, calculate the analytical derivative of the phase function
-      return 0.5 / pressure_width * (1.0 - std::tanh(pressure_deviation / pressure_width)
-    				                     * std::tanh(pressure_deviation / pressure_width));
+  	  if(transition_widths[phase]==0)
+  	    return 0;
+  	  else
+        return 0.5 / pressure_width * (1.0 - std::tanh(pressure_deviation / pressure_width)
+    				                       * std::tanh(pressure_deviation / pressure_width));
   	}
 
     template <int dim>
@@ -198,11 +219,14 @@
              const Point<dim> &position) const
     {
       // first, calculate temperature dependence of density
-	  double temperature_dependence = 1.0;
-	  if (this->include_adiabatic_heating ())
-        // temperature dependence is 1 - alpha * (T - T(adiabatic))
-        temperature_dependence -= (temperature - this->get_adiabatic_conditions().temperature(position))
-	    	      * thermal_expansion_coefficient(temperature, pressure, compositional_fields, position);
+  	  double temperature_dependence = 1.0;
+  	  if (this->include_adiabatic_heating ())
+  	  {
+          // temperature dependence is 1 - alpha * (T - T(adiabatic))
+  	    if (&this->get_adiabatic_conditions())
+            temperature_dependence -= (temperature - this->get_adiabatic_conditions().temperature(position))
+  	    	        * thermal_expansion_coefficient(temperature, pressure, compositional_fields, position);
+  	  }
 	  else
 		temperature_dependence -= temperature * thermal_expansion_coefficient(temperature, pressure, compositional_fields, position);
 
@@ -295,7 +319,7 @@
                         const NonlinearDependence::Dependence dependence) const
     {
       double entropy_gradient = 0.0;
-      if (&this->get_adiabatic_conditions())
+      if (&this->get_adiabatic_conditions() && this->include_latent_heat())
         for (unsigned int phase=0;phase<transition_depths.size();++phase)
         {
     	  //calculate derivative of the phase function

Modified: branches/j-dannberg/source/material_model/latent_heat.cc
===================================================================
--- branches/j-dannberg/source/material_model/latent_heat.cc	2013-09-11 14:21:08 UTC (rev 1875)
+++ branches/j-dannberg/source/material_model/latent_heat.cc	2013-09-11 14:28:30 UTC (rev 1876)
@@ -56,7 +56,13 @@
 
         // last, calculate the percentage of material that has undergone the transition
         // (also in dependence of the phase transition width - this is an input parameter)
-        return 0.5*(1.0 + std::tanh(pressure_deviation / pressure_width));
+    	double phase_func;
+    	// use delta function for width = 0
+    	if(transition_widths[phase]==0)
+    	  (pressure_deviation > 0) ? phase_func = 1 : phase_func = 0;
+    	else
+    	  phase_func = 0.5*(1.0 + std::tanh(pressure_deviation / pressure_width));
+        return phase_func;
       }
       // if we do not have the adiabatic conditions, we have to use the depth instead
       // this is less precise, because we do not have the exact pressure gradient, instead we use pressure/depth
@@ -64,9 +70,21 @@
       else
       {
     	double depth = this->get_geometry_model().depth(position);
-    	double depth_deviation = depth - transition_depths[phase]
-    	                         - transition_slopes[phase] * (depth / pressure) * (temperature - transition_temperatures[phase]);
-    	return 0.5*(1.0 + std::tanh(depth_deviation / transition_widths[phase]));
+    	double depth_deviation = (pressure > 0
+    			                  ?
+    			                  depth - transition_depths[phase]
+    	                          - transition_slopes[phase] * (depth / pressure) * (temperature - transition_temperatures[phase])
+    	                          :
+    			                  depth - transition_depths[phase]
+    	                          - transition_slopes[phase] / (this->get_gravity_model().gravity_vector(position).norm() * reference_rho)
+    	                          * (temperature - transition_temperatures[phase]));
+    	double phase_func;
+    	// use delta function for width = 0
+    	if(transition_widths[phase]==0)
+    	  (depth_deviation > 0) ? phase_func = 1 : phase_func = 0;
+    	else
+    	  phase_func = 0.5*(1.0 + std::tanh(depth_deviation / transition_widths[phase]));
+    	return phase_func;
       }
   	}
 
@@ -96,8 +114,11 @@
                                 - transition_slopes[phase] * (temperature - transition_temperatures[phase]);
 
       // last, calculate the analytical derivative of the phase function
-      return 0.5 / pressure_width * (1.0 - std::tanh(pressure_deviation / pressure_width)
-    				                     * std::tanh(pressure_deviation / pressure_width));
+  	  if(transition_widths[phase]==0)
+  	    return 0;
+  	  else
+        return 0.5 / pressure_width * (1.0 - std::tanh(pressure_deviation / pressure_width)
+    				                       * std::tanh(pressure_deviation / pressure_width));
   	}
 
     template <int dim>
@@ -200,9 +221,12 @@
       // first, calculate temperature dependence of density
 	  double temperature_dependence = 1.0;
 	  if (this->include_adiabatic_heating ())
+	  {
         // temperature dependence is 1 - alpha * (T - T(adiabatic))
-        temperature_dependence -= (temperature - this->get_adiabatic_conditions().temperature(position))
-	    	      * thermal_expansion_coefficient(temperature, pressure, compositional_fields, position);
+	    if (&this->get_adiabatic_conditions())
+          temperature_dependence -= (temperature - this->get_adiabatic_conditions().temperature(position))
+	    	        * thermal_expansion_coefficient(temperature, pressure, compositional_fields, position);
+	  }
 	  else
 		temperature_dependence -= temperature * thermal_expansion_coefficient(temperature, pressure, compositional_fields, position);
 
@@ -295,7 +319,7 @@
                         const NonlinearDependence::Dependence dependence) const
     {
       double entropy_gradient = 0.0;
-      if (&this->get_adiabatic_conditions())
+      if (&this->get_adiabatic_conditions() && this->include_latent_heat())
         for (unsigned int phase=0;phase<transition_depths.size();++phase)
         {
     	  //calculate derivative of the phase function

Modified: branches/j-dannberg/source/simulator/simulator_access.cc
===================================================================
--- branches/j-dannberg/source/simulator/simulator_access.cc	2013-09-11 14:21:08 UTC (rev 1875)
+++ branches/j-dannberg/source/simulator/simulator_access.cc	2013-09-11 14:28:30 UTC (rev 1876)
@@ -145,6 +145,14 @@
   }
 
   template <int dim>
+  bool
+  SimulatorAccess<dim>::include_latent_heat () const
+  {
+    return simulator->parameters.include_latent_heat;
+  }
+
+
+  template <int dim>
   double
   SimulatorAccess<dim>::get_adiabatic_surface_temperature () const
   {