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

dealii.demon at gmail.com dealii.demon at gmail.com
Fri Nov 22 14:07:12 PST 2013


Revision 2040

Merge from mainline.

_U  branches/j-dannberg/
A   branches/j-dannberg/CTestConfig.cmake
U   branches/j-dannberg/Makefile
U   branches/j-dannberg/doc/modules/changes.h
A   branches/j-dannberg/include/aspect/postprocess/visualization/melt_fraction.h
U   branches/j-dannberg/include/aspect/simulator.h
U   branches/j-dannberg/include/aspect/termination_criteria/end_time.h
U   branches/j-dannberg/include/aspect/termination_criteria/interface.h
U   branches/j-dannberg/source/geometry_model/box.cc
A   branches/j-dannberg/source/postprocess/visualization/melt_fraction.cc
U   branches/j-dannberg/source/postprocess/visualization.cc
U   branches/j-dannberg/source/simulator/assembly.cc
U   branches/j-dannberg/source/simulator/core.cc
U   branches/j-dannberg/source/simulator/helper_functions.cc
U   branches/j-dannberg/source/termination_criteria/end_time.cc
U   branches/j-dannberg/source/termination_criteria/interface.cc
U   branches/j-dannberg/tests/CMakeLists.txt
U   branches/j-dannberg/tests/Makefile
U   branches/j-dannberg/tests/adiabatic_conditions/screen-output
U   branches/j-dannberg/tests/adiabatic_initial_conditions/screen-output
U   branches/j-dannberg/tests/adiabatic_initial_conditions_subadiabaticity/screen-output
U   branches/j-dannberg/tests/box-end_time_1e7-terminate/screen-output
U   branches/j-dannberg/tests/box-end_time_1e7-terminate/statistics
U   branches/j-dannberg/tests/box-first-time-step/screen-output
U   branches/j-dannberg/tests/box-first-time-step-alternate-bc/screen-output
U   branches/j-dannberg/tests/box-steady-state-terminate/screen-output
U   branches/j-dannberg/tests/composition-active/screen-output
U   branches/j-dannberg/tests/composition-active/statistics
U   branches/j-dannberg/tests/composition-passive-tracers/particles.pvd
U   branches/j-dannberg/tests/composition-passive-tracers/screen-output
U   branches/j-dannberg/tests/composition-passive-tracers/statistics
A   branches/j-dannberg/tests/compressibility/
D   branches/j-dannberg/tests/compressibility/screen-output
A   branches/j-dannberg/tests/compressibility/screen-output
A   branches/j-dannberg/tests/compressibility.cc
A   branches/j-dannberg/tests/compressibility.prm
A   branches/j-dannberg/tests/compressibility_iterated_stokes/
D   branches/j-dannberg/tests/compressibility_iterated_stokes/screen-output
A   branches/j-dannberg/tests/compressibility_iterated_stokes/screen-output
A   branches/j-dannberg/tests/compressibility_iterated_stokes.cc
A   branches/j-dannberg/tests/compressibility_iterated_stokes.prm
U   branches/j-dannberg/tests/conservative_with_mpi/screen-output
U   branches/j-dannberg/tests/diffusion/screen-output
U   branches/j-dannberg/tests/diffusion/statistics
U   branches/j-dannberg/tests/diffusion-velocity/screen-output
U   branches/j-dannberg/tests/diffusion-velocity/statistics
U   branches/j-dannberg/tests/graphical_output/screen-output
U   branches/j-dannberg/tests/graphical_output/solution-00000.0000.gnuplot
U   branches/j-dannberg/tests/inclusion_2/screen-output
U   branches/j-dannberg/tests/inclusion_4/screen-output
U   branches/j-dannberg/tests/inclusion_adaptive/screen-output
U   branches/j-dannberg/tests/no_flow/screen-output
U   branches/j-dannberg/tests/no_flow/statistics
U   branches/j-dannberg/tests/non_conservative_with_mpi/screen-output
U   branches/j-dannberg/tests/passive_comp/depthaverage.plt
U   branches/j-dannberg/tests/passive_comp/screen-output
U   branches/j-dannberg/tests/passive_comp/statistics
U   branches/j-dannberg/tests/plugin/screen-output
A   branches/j-dannberg/tests/run_testsuite.cmake
U   branches/j-dannberg/tests/sol_cx_2/screen-output
U   branches/j-dannberg/tests/sol_cx_2_conservative/screen-output
U   branches/j-dannberg/tests/sol_cx_2_normalized_pressure/screen-output
U   branches/j-dannberg/tests/sol_cx_2_q3/screen-output
U   branches/j-dannberg/tests/sol_cx_2_q3/statistics
U   branches/j-dannberg/tests/sol_cx_4/screen-output
U   branches/j-dannberg/tests/sol_cx_4_conservative/screen-output
U   branches/j-dannberg/tests/sol_cx_4_normalized_pressure/screen-output
U   branches/j-dannberg/tests/sol_cx_4_normalized_pressure_large_static_pressure/screen-output
U   branches/j-dannberg/tests/sol_cx_4_normalized_pressure_low_solver_tolerance/screen-output
U   branches/j-dannberg/tests/sol_cx_tracers/screen-output
U   branches/j-dannberg/tests/sol_cx_tracers/statistics
U   branches/j-dannberg/tests/sol_kz_2/screen-output
U   branches/j-dannberg/tests/sol_kz_2_cheaper_first_phase_solver/screen-output
U   branches/j-dannberg/tests/sol_kz_2_conservative/screen-output
U   branches/j-dannberg/tests/sol_kz_2_no_first_phase_solver/screen-output
U   branches/j-dannberg/tests/sol_kz_2_q3/screen-output
U   branches/j-dannberg/tests/sol_kz_2_q3/statistics
U   branches/j-dannberg/tests/sol_kz_4/screen-output
U   branches/j-dannberg/tests/sol_kz_4_conservative/screen-output
U   branches/j-dannberg/tests/stokes/screen-output
U   branches/j-dannberg/tests/temperature-dependent-stokes-matrix/screen-output
U   branches/j-dannberg/tests/temperature-dependent-stokes-matrix/statistics
U   branches/j-dannberg/tests/terminate_user_request/screen-output
U   branches/j-dannberg/tests/velocity_in_years/screen-output


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

Diff:
Index: branches/j-dannberg
===================================================================
--- branches/j-dannberg	2013-11-22 21:35:22 UTC (rev 2039)
+++ branches/j-dannberg	2013-11-22 22:06:38 UTC (rev 2040)

Property changes on: branches/j-dannberg
___________________________________________________________________
Modified: svn:mergeinfo
## -1,4 +1,4 ##
 /branches/active_compositions:1257-1296
 /branches/compositional:1141-1251
 /branches/fully-nonlinear:542-728
-/trunk/aspect:1677-2012
+/trunk/aspect:1677-2039
\ No newline at end of property
Copied: branches/j-dannberg/CTestConfig.cmake (from rev 2039, trunk/aspect/CTestConfig.cmake)
===================================================================
--- branches/j-dannberg/CTestConfig.cmake	                        (rev 0)
+++ branches/j-dannberg/CTestConfig.cmake	2013-11-22 22:06:38 UTC (rev 2040)
@@ -0,0 +1,54 @@
+## ---------------------------------------------------------------------
+## $Id$
+##
+## Copyright (C) 2013 by the authors of the ASPECT code.
+##
+##  This file is part of ASPECT.
+##
+##  ASPECT is free software; you can redistribute it and/or modify
+##  it under the terms of the GNU General Public License as published by
+##  the Free Software Foundation; either version 2, or (at your option)
+##  any later version.
+##
+##  ASPECT is distributed in the hope that it will be useful,
+##  but WITHOUT ANY WARRANTY; without even the implied warranty of
+##  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+##  GNU General Public License for more details.
+##
+##  You should have received a copy of the GNU General Public License
+##  along with ASPECT; see the file doc/COPYING.  If not see
+##  <http://www.gnu.org/licenses/>.
+##
+## ---------------------------------------------------------------------
+
+#
+# Dashboard configuration:
+#
+
+SET(CTEST_PROJECT_NAME "aspect")
+
+SET(CTEST_DROP_METHOD "http")
+SET(CTEST_DROP_SITE "cdash.kyomu.43-1.org")
+SET(CTEST_DROP_LOCATION "/submit.php?project=aspect")
+SET(CTEST_DROP_SITE_CDASH TRUE)
+
+SET(CTEST_CUSTOM_MAXIMUM_NUMBER_OF_ERRORS   100)
+SET(CTEST_CUSTOM_MAXIMUM_NUMBER_OF_WARNINGS 300)
+
+# number of lines to submit before an error:
+SET(CTEST_CUSTOM_ERROR_PRE_CONTEXT            5)
+# number of lines to submit after an error:
+SET(CTEST_CUSTOM_ERROR_POST_CONTEXT          20)
+
+#
+# Coverage options:
+#
+
+SET(CTEST_EXTRA_COVERAGE_GLOB
+  # These files should have executable lines and therefore coverage:
+  # source/**/*.cc
+  )
+
+SET(CTEST_CUSTOM_COVERAGE_EXCLUDE
+  "/tests/"
+  )

Modified: branches/j-dannberg/Makefile
===================================================================
--- branches/j-dannberg/Makefile	2013-11-22 21:35:22 UTC (rev 2039)
+++ branches/j-dannberg/Makefile	2013-11-22 22:06:38 UTC (rev 2040)
@@ -30,10 +30,10 @@
 SHELL = /bin/sh
 
 # The CMake executable.
-CMAKE_COMMAND = /w/bangerth/share/software/cmake-2.8.9/bin/cmake
+CMAKE_COMMAND = /w/bangerth/share/software/bin/cmake
 
 # The command to remove a file.
-RM = /w/bangerth/share/software/cmake-2.8.9/bin/cmake -E remove -f
+RM = /w/bangerth/share/software/bin/cmake -E remove -f
 
 # Escaping for special characters.
 EQUALS = =
@@ -63,7 +63,7 @@
 # Special rule for the target rebuild_cache
 rebuild_cache:
 	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running CMake to regenerate build system..."
-	/w/bangerth/share/software/cmake-2.8.9/bin/cmake -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
+	/w/bangerth/share/software/bin/cmake -H$(CMAKE_SOURCE_DIR) -B$(CMAKE_BINARY_DIR)
 .PHONY : rebuild_cache
 
 # Special rule for the target rebuild_cache
@@ -73,7 +73,7 @@
 # Special rule for the target test
 test:
 	@$(CMAKE_COMMAND) -E cmake_echo_color --switch=$(COLOR) --cyan "Running tests..."
-	/w/bangerth/share/software/cmake-2.8.9/bin/ctest --force-new-ctest-process $(ARGS)
+	/w/bangerth/share/software/bin/ctest --force-new-ctest-process $(ARGS)
 .PHONY : test
 
 # Special rule for the target test
@@ -671,6 +671,84 @@
 .PHONY : composition-passive-tracers.statistics.diff/fast
 
 #=============================================================================
+# Target rules for targets named compressibility
+
+# Build rule for target.
+compressibility: cmake_check_build_system
+	$(MAKE) -f CMakeFiles/Makefile2 compressibility
+.PHONY : compressibility
+
+# fast build rule for target.
+compressibility/fast:
+	$(MAKE) -f tests/CMakeFiles/compressibility.dir/build.make tests/CMakeFiles/compressibility.dir/build
+.PHONY : compressibility/fast
+
+#=============================================================================
+# Target rules for targets named compressibility.run
+
+# Build rule for target.
+compressibility.run: cmake_check_build_system
+	$(MAKE) -f CMakeFiles/Makefile2 compressibility.run
+.PHONY : compressibility.run
+
+# fast build rule for target.
+compressibility.run/fast:
+	$(MAKE) -f tests/CMakeFiles/compressibility.run.dir/build.make tests/CMakeFiles/compressibility.run.dir/build
+.PHONY : compressibility.run/fast
+
+#=============================================================================
+# Target rules for targets named compressibility.screen-output.diff
+
+# Build rule for target.
+compressibility.screen-output.diff: cmake_check_build_system
+	$(MAKE) -f CMakeFiles/Makefile2 compressibility.screen-output.diff
+.PHONY : compressibility.screen-output.diff
+
+# fast build rule for target.
+compressibility.screen-output.diff/fast:
+	$(MAKE) -f tests/CMakeFiles/compressibility.screen-output.diff.dir/build.make tests/CMakeFiles/compressibility.screen-output.diff.dir/build
+.PHONY : compressibility.screen-output.diff/fast
+
+#=============================================================================
+# Target rules for targets named compressibility_iterated_stokes
+
+# Build rule for target.
+compressibility_iterated_stokes: cmake_check_build_system
+	$(MAKE) -f CMakeFiles/Makefile2 compressibility_iterated_stokes
+.PHONY : compressibility_iterated_stokes
+
+# fast build rule for target.
+compressibility_iterated_stokes/fast:
+	$(MAKE) -f tests/CMakeFiles/compressibility_iterated_stokes.dir/build.make tests/CMakeFiles/compressibility_iterated_stokes.dir/build
+.PHONY : compressibility_iterated_stokes/fast
+
+#=============================================================================
+# Target rules for targets named compressibility_iterated_stokes.run
+
+# Build rule for target.
+compressibility_iterated_stokes.run: cmake_check_build_system
+	$(MAKE) -f CMakeFiles/Makefile2 compressibility_iterated_stokes.run
+.PHONY : compressibility_iterated_stokes.run
+
+# fast build rule for target.
+compressibility_iterated_stokes.run/fast:
+	$(MAKE) -f tests/CMakeFiles/compressibility_iterated_stokes.run.dir/build.make tests/CMakeFiles/compressibility_iterated_stokes.run.dir/build
+.PHONY : compressibility_iterated_stokes.run/fast
+
+#=============================================================================
+# Target rules for targets named compressibility_iterated_stokes.screen-output.diff
+
+# Build rule for target.
+compressibility_iterated_stokes.screen-output.diff: cmake_check_build_system
+	$(MAKE) -f CMakeFiles/Makefile2 compressibility_iterated_stokes.screen-output.diff
+.PHONY : compressibility_iterated_stokes.screen-output.diff
+
+# fast build rule for target.
+compressibility_iterated_stokes.screen-output.diff/fast:
+	$(MAKE) -f tests/CMakeFiles/compressibility_iterated_stokes.screen-output.diff.dir/build.make tests/CMakeFiles/compressibility_iterated_stokes.screen-output.diff.dir/build
+.PHONY : compressibility_iterated_stokes.screen-output.diff/fast
+
+#=============================================================================
 # Target rules for targets named conservative_with_mpi.run
 
 # Build rule for target.
@@ -2205,6 +2283,32 @@
 .PHONY : tests.composition-passive-tracers/fast
 
 #=============================================================================
+# Target rules for targets named tests.compressibility
+
+# Build rule for target.
+tests.compressibility: cmake_check_build_system
+	$(MAKE) -f CMakeFiles/Makefile2 tests.compressibility
+.PHONY : tests.compressibility
+
+# fast build rule for target.
+tests.compressibility/fast:
+	$(MAKE) -f tests/CMakeFiles/tests.compressibility.dir/build.make tests/CMakeFiles/tests.compressibility.dir/build
+.PHONY : tests.compressibility/fast
+
+#=============================================================================
+# Target rules for targets named tests.compressibility_iterated_stokes
+
+# Build rule for target.
+tests.compressibility_iterated_stokes: cmake_check_build_system
+	$(MAKE) -f CMakeFiles/Makefile2 tests.compressibility_iterated_stokes
+.PHONY : tests.compressibility_iterated_stokes
+
+# fast build rule for target.
+tests.compressibility_iterated_stokes/fast:
+	$(MAKE) -f tests/CMakeFiles/tests.compressibility_iterated_stokes.dir/build.make tests/CMakeFiles/tests.compressibility_iterated_stokes.dir/build
+.PHONY : tests.compressibility_iterated_stokes/fast
+
+#=============================================================================
 # Target rules for targets named tests.conservative_with_mpi
 
 # Build rule for target.
@@ -4208,6 +4312,30 @@
 	$(MAKE) -f CMakeFiles/aspect.dir/build.make CMakeFiles/aspect.dir/source/postprocess/visualization/friction_heating.cc.s
 .PHONY : source/postprocess/visualization/friction_heating.cc.s
 
+source/postprocess/visualization/melt_fraction.o: source/postprocess/visualization/melt_fraction.cc.o
+.PHONY : source/postprocess/visualization/melt_fraction.o
+
+# target to build an object file
+source/postprocess/visualization/melt_fraction.cc.o:
+	$(MAKE) -f CMakeFiles/aspect.dir/build.make CMakeFiles/aspect.dir/source/postprocess/visualization/melt_fraction.cc.o
+.PHONY : source/postprocess/visualization/melt_fraction.cc.o
+
+source/postprocess/visualization/melt_fraction.i: source/postprocess/visualization/melt_fraction.cc.i
+.PHONY : source/postprocess/visualization/melt_fraction.i
+
+# target to preprocess a source file
+source/postprocess/visualization/melt_fraction.cc.i:
+	$(MAKE) -f CMakeFiles/aspect.dir/build.make CMakeFiles/aspect.dir/source/postprocess/visualization/melt_fraction.cc.i
+.PHONY : source/postprocess/visualization/melt_fraction.cc.i
+
+source/postprocess/visualization/melt_fraction.s: source/postprocess/visualization/melt_fraction.cc.s
+.PHONY : source/postprocess/visualization/melt_fraction.s
+
+# target to generate assembly for a file
+source/postprocess/visualization/melt_fraction.cc.s:
+	$(MAKE) -f CMakeFiles/aspect.dir/build.make CMakeFiles/aspect.dir/source/postprocess/visualization/melt_fraction.cc.s
+.PHONY : source/postprocess/visualization/melt_fraction.cc.s
+
 source/postprocess/visualization/nonadiabatic_pressure.o: source/postprocess/visualization/nonadiabatic_pressure.cc.o
 .PHONY : source/postprocess/visualization/nonadiabatic_pressure.o
 

Modified: branches/j-dannberg/doc/modules/changes.h
===================================================================
--- branches/j-dannberg/doc/modules/changes.h	2013-11-22 21:35:22 UTC (rev 2039)
+++ branches/j-dannberg/doc/modules/changes.h	2013-11-22 22:06:38 UTC (rev 2040)
@@ -8,6 +8,27 @@
 </p>
 
 <ol>
+  <li>Fixed: When using compressible models with nonlinear iterations
+  such as "Stokes", "iterated IMPES" or "iterated Stokes" and prescribed
+  boundary values, there were numerous bugs that should now be fixed.
+  <br>
+  (Wolfgang Bangerth 2013/11/21)
+
+  <li>Changed: When the user selects to terminate by end time, the
+  final time step is adjusted to hit the final time exactly.
+  <br>
+  (Ryan Grove 2013/11/19)
+
+  <li>Fixed: When using compressible models, we fixed up the right hand side
+  vector in a way so that the mean divergence is zero (even though it is of
+  course locally nonzero due to the compressibility). This is necessary to ensure
+  the solvability of the linear system, but it is wrong if the domain has open
+  boundaries through which fluid can escape or enter. We now only perform
+  this correction if there are no open boundaries and no boundaries with a
+  prescribed velocity.
+  <br>
+  (Wolfgang Bangerth 2013/11/19)
+
   <li>New: It is now possible to prescribe the velocity only for certain
   components in the 'Prescribed velocity boundary indicators' parameter.
   <br>

Copied: branches/j-dannberg/include/aspect/postprocess/visualization/melt_fraction.h (from rev 2039, trunk/aspect/include/aspect/postprocess/visualization/melt_fraction.h)
===================================================================
--- branches/j-dannberg/include/aspect/postprocess/visualization/melt_fraction.h	                        (rev 0)
+++ branches/j-dannberg/include/aspect/postprocess/visualization/melt_fraction.h	2013-11-22 22:06:38 UTC (rev 2040)
@@ -0,0 +1,124 @@
+/*
+  Copyright (C) 2013 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file doc/COPYING.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+/*  $Id: thermodynamic_phase.h 1433 2012-12-08 08:24:55Z bangerth $  */
+
+
+#ifndef __aspect__postprocess_visualization_melt_fraction_h
+#define __aspect__postprocess_visualization_melt_fraction_h
+
+#include <aspect/postprocess/visualization.h>
+#include <aspect/simulator_access.h>
+
+#include <deal.II/numerics/data_out.h>
+
+
+namespace aspect
+{
+  namespace Postprocess
+  {
+    namespace VisualizationPostprocessors
+    {
+      /**
+       * A class derived from DataPostprocessor that takes an output vector and
+       * computes a variable that represents the melt fraction
+       * at every point.
+       *
+       * The member functions are all implementations of those declared in the base
+       * class. See there for their meaning.
+       */
+      template <int dim>
+      class MeltFraction
+        : public DataPostprocessorScalar<dim>,
+          public SimulatorAccess<dim>,
+          public Interface<dim>
+      {
+        public:
+          MeltFraction ();
+
+          virtual
+          void
+          compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                             const std::vector<std::vector<Tensor<1,dim> > > &duh,
+                                             const std::vector<std::vector<Tensor<2,dim> > > &dduh,
+                                             const std::vector<Point<dim> >                  &normals,
+                                             const std::vector<Point<dim> >                  &evaluation_points,
+                                             std::vector<Vector<double> >                    &computed_quantities) const;
+
+          /**
+           * Declare the parameters this class takes through input files.
+           */
+          static
+          void
+          declare_parameters (ParameterHandler &prm);
+
+          /**
+           * Read the parameters this class declares from the parameter
+           * file.
+           */
+          virtual
+          void
+          parse_parameters (ParameterHandler &prm);
+
+        private:
+          /**
+           * Parameters for anhydrous melting of peridotite after Katz, 2003
+           */
+
+          // for the solidus temperature
+          double A1;   // °C
+          double A2; // °C/Pa
+          double A3; // °C/(Pa^2)
+
+          // for the lherzolite liquidus temperature
+          double B1;   // °C
+          double B2;   // °C/Pa
+          double B3; // °C/(Pa^2)
+
+          // for the liquidus temperature
+          double C1;   // °C
+          double C2;  // °C/Pa
+          double C3; // °C/(Pa^2)
+
+          // for the reaction coefficient of pyroxene
+          double r1;     // cpx/melt
+          double r2;     // cpx/melt/GPa
+          double M_cpx;  // mass fraction of pyroxenite
+
+          // melt fraction exponent
+          double beta;
+
+          /**
+           * Parameters for melting of pyroxenite after Sobolev et al., 2011
+           */
+
+          // for the melting temperature
+          double D1;    // °C
+          double D2;  // °C/Pa
+          double D3; // °C/(Pa^2)
+
+          // for the melt-fraction dependence of productivity
+          double E1;
+          double E2;
+      };
+    }
+  }
+}
+
+#endif

Modified: branches/j-dannberg/include/aspect/simulator.h
===================================================================
--- branches/j-dannberg/include/aspect/simulator.h	2013-11-22 21:35:22 UTC (rev 2039)
+++ branches/j-dannberg/include/aspect/simulator.h	2013-11-22 22:06:38 UTC (rev 2040)
@@ -251,7 +251,7 @@
          */
 
         /**
-         * @name Parameters that have to do with compositional field
+         * @name Parameters that have to do with compositional fields
          * @{
          */
         unsigned int                   n_compositional_fields;
@@ -1105,7 +1105,7 @@
       double                                                    global_Omega_diameter;
       double                                                    global_volume;
 
-      MeshRefinement::Manager<dim>        mesh_refinement_manager;
+      MeshRefinement::Manager<dim>                              mesh_refinement_manager;
 
       const MappingQ<dim>                                       mapping;
 
@@ -1116,10 +1116,28 @@
       ConstraintMatrix                                          constraints;
       ConstraintMatrix                                          current_constraints;
 
+      /**
+       * The latest correction computed by normalize_pressure(). We store this so
+       * we can undo the correction in denormalize_pressure().
+       */
       double                                                    pressure_adjustment;
+
+      /**
+       * Scaling factor for the pressure as explained in the Kronbichler/Heister/Bangerth
+       * paper to ensure that the linear system that results from the Stokes equations
+       * is well conditioned.
+       */
       double                                                    pressure_scaling;
 
       /**
+       * A variable that determines whether we need to do the correction of
+       * the Stokes right hand side vector to ensure that the average divergence
+       * is zero. This is necessary for compressible models, but only if there
+       * are no in/outflow boundaries.
+       */
+      bool                           do_pressure_rhs_compatibility_modification;
+
+      /**
        * @}
        */
 

Modified: branches/j-dannberg/include/aspect/termination_criteria/end_time.h
===================================================================
--- branches/j-dannberg/include/aspect/termination_criteria/end_time.h	2013-11-22 21:35:22 UTC (rev 2039)
+++ branches/j-dannberg/include/aspect/termination_criteria/end_time.h	2013-11-22 22:06:38 UTC (rev 2040)
@@ -41,7 +41,13 @@
     class EndTime : public Interface<dim>, public SimulatorAccess<dim>
     {
       public:
+       
         /**
+         * Check this termination criterion and possibly reduce time step size
+         */
+        double check_for_last_time_step (double time_step) const;
+
+        /**
          * Evaluate this termination criterion.
          *
          * @return Whether to terminate the simulation (true) or continue (false).

Modified: branches/j-dannberg/include/aspect/termination_criteria/interface.h
===================================================================
--- branches/j-dannberg/include/aspect/termination_criteria/interface.h	2013-11-22 21:35:22 UTC (rev 2039)
+++ branches/j-dannberg/include/aspect/termination_criteria/interface.h	2013-11-22 22:06:38 UTC (rev 2040)
@@ -91,6 +91,19 @@
         bool
         execute () = 0;
 
+         /**
+         * Check for last time step and if so reduce the time step to user 
+         * specified end time
+         *
+         * @return Reduced or current time step size
+         *
+         * @note A reduced time step size will be returned for the
+         * last time step but the current time step size will be
+         * returned otherwise. Also, the default implementation will
+         * always return the full time step size.
+         */
+        virtual double check_for_last_time_step (double time_step) const;
+
         /**
          * Declare the parameters this class takes through input files.
          * Derived classes should overload this function if they actually
@@ -162,6 +175,20 @@
         std::pair<bool,bool>
         execute () const;
 
+         /**
+         * Check all of the termination criteria objects that have
+         * been requested in the input file for criteria regarding
+         * last time step and if so get the minimum of these values.
+         *
+         * @return Reduced or current time step size
+         *
+         * @note A reduced time step size will be returned for the
+         * last time step but the current time step size will be
+         * returned otherwise. The time step will be greater than zero
+         * as well as less than or equal to the inputted time step
+         */
+        double check_for_last_time_step (double time_step) const;
+
         /**
          * Declare the parameters of all known termination criteria plugins, as
          * well as of ones this class has itself.

Modified: branches/j-dannberg/source/geometry_model/box.cc
===================================================================
--- branches/j-dannberg/source/geometry_model/box.cc	2013-11-22 21:35:22 UTC (rev 2039)
+++ branches/j-dannberg/source/geometry_model/box.cc	2013-11-22 22:06:38 UTC (rev 2040)
@@ -122,8 +122,18 @@
     {
       const double d = maximal_depth()-position(dim-1);
 
-      Assert (d >= -1e-14*std::fabs(maximal_depth()), ExcInternalError());
-      Assert (d <= (1.+1e-14)*maximal_depth(), ExcInternalError());
+      // if we violate the bounds, check that we do so only very slightly and
+      // then just return maximal or minimal depth
+      if (d < 0)
+        {
+          Assert (d >= -1e-14*std::fabs(maximal_depth()), ExcInternalError());
+          return 0;
+        }
+      if (d > maximal_depth())
+        {
+          Assert (d <= (1.+1e-14)*maximal_depth(), ExcInternalError());
+          return maximal_depth();
+        }
 
       return d;
     }

Copied: branches/j-dannberg/source/postprocess/visualization/melt_fraction.cc (from rev 2039, trunk/aspect/source/postprocess/visualization/melt_fraction.cc)
===================================================================
--- branches/j-dannberg/source/postprocess/visualization/melt_fraction.cc	                        (rev 0)
+++ branches/j-dannberg/source/postprocess/visualization/melt_fraction.cc	2013-11-22 22:06:38 UTC (rev 2040)
@@ -0,0 +1,312 @@
+/*
+  Copyright (C) 2013 by the authors of the ASPECT code.
+
+  This file is part of ASPECT.
+
+  ASPECT is free software; you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation; either version 2, or (at your option)
+  any later version.
+
+  ASPECT is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with ASPECT; see the file doc/COPYING.  If not see
+  <http://www.gnu.org/licenses/>.
+*/
+/*  $Id: thermodynamic_phase.cc 1433 2012-12-08 08:24:55Z bangerth $  */
+
+
+#include <aspect/postprocess/visualization/melt_fraction.h>
+#include <aspect/simulator_access.h>
+
+#include <deal.II/numerics/data_out.h>
+#include <deal.II/base/parameter_handler.h>
+
+
+namespace aspect
+{
+  namespace Postprocess
+  {
+    namespace VisualizationPostprocessors
+    {
+      template <int dim>
+      MeltFraction<dim>::
+      MeltFraction ()
+        :
+        DataPostprocessorScalar<dim> ("melt_fraction",
+                                      update_values | update_q_points)
+      {}
+
+
+
+      template <int dim>
+      void
+      MeltFraction<dim>::
+      compute_derived_quantities_vector (const std::vector<Vector<double> >              &uh,
+                                         const std::vector<std::vector<Tensor<1,dim> > > &duh,
+                                         const std::vector<std::vector<Tensor<2,dim> > > &dduh,
+                                         const std::vector<Point<dim> >                  &normals,
+                                         const std::vector<Point<dim> >                  &evaluation_points,
+                                         std::vector<Vector<double> >                    &computed_quantities) const
+      {
+        const unsigned int n_quadrature_points = uh.size();
+        Assert (computed_quantities.size() == n_quadrature_points,    ExcInternalError());
+        Assert (computed_quantities[0].size() == 1,                   ExcInternalError());
+        Assert (uh[0].size() == dim+2+this->n_compositional_fields(), ExcInternalError());
+
+        for (unsigned int q=0; q<n_quadrature_points; ++q)
+          {
+            const double pressure    = uh[q][dim];
+            const double temperature = uh[q][dim+1];
+            std::vector<double> composition(this->n_compositional_fields());
+
+            for (unsigned int c=0; c<this->n_compositional_fields(); ++c)
+              composition[c] = uh[q][dim+2+c];
+
+            // anhydrous melting of peridotite after Katz, 2003
+            const double T_solidus  = A1 + 273.15
+            		                + A2 * pressure
+            		                + A3 * pressure * pressure;
+            const double T_lherz_liquidus = B1 + 273.15
+            		                      + B2 * pressure
+            		                      + B3 * pressure * pressure;
+            const double T_liquidus = C1 + 273.15
+            		                + C2 * pressure
+            		                + C3 * pressure * pressure;
+
+            // melt fraction for peridotite with clinopyroxene
+            double peridotite_melt_fraction;
+            if (temperature < T_solidus || pressure > 1.3e10)
+              peridotite_melt_fraction = 0.0;
+            else if (temperature > T_liquidus)
+              peridotite_melt_fraction = 1.0;
+            else
+              peridotite_melt_fraction = std::pow((temperature - T_solidus) / (T_lherz_liquidus - T_solidus),beta);
+
+            // melt fraction after melting of all clinopyroxene
+            const double R_cpx = r1 + r2 * pressure;
+            const double F_max = M_cpx / R_cpx;
+
+            if(peridotite_melt_fraction > F_max && peridotite_melt_fraction < 1.0)
+            {
+              const double T_max = std::pow(F_max,1/beta) * (T_lherz_liquidus - T_solidus) + T_solidus;
+              peridotite_melt_fraction = F_max + (1 - F_max) * (temperature - T_max) / (T_liquidus - T_max);
+            }
+
+            // melting of pyroxenite after Sobolev et al., 2011
+            const double T_melting = D1 + 273.15
+  		                           + D2 * pressure
+  		                           + D3 * pressure * pressure;
+
+            const double discriminant = E1*E1/(E2*E2*4) + (temperature-T_melting)/E2;
+
+            double pyroxenite_melt_fraction;
+            if (temperature < T_melting || pressure > 1.3e10)
+              pyroxenite_melt_fraction = 0.0;
+            else if (discriminant < 0)
+              pyroxenite_melt_fraction = 0.5429;
+            else
+              pyroxenite_melt_fraction = -E1/(2*E2) - std::sqrt(discriminant);
+
+            double melt_fraction;
+            if(this->n_compositional_fields()>0)
+              melt_fraction = composition[0] * pyroxenite_melt_fraction +
+                              (1-composition[0]) * peridotite_melt_fraction;
+            else
+              melt_fraction = peridotite_melt_fraction;
+
+            computed_quantities[q](0) = melt_fraction;
+          }
+      }
+
+
+
+      template <int dim>
+      void
+      MeltFraction<dim>::declare_parameters (ParameterHandler &prm)
+	  {
+		prm.enter_subsection("Postprocess");
+		{
+		  prm.enter_subsection("Visualization");
+		  {
+			prm.enter_subsection("Melt fraction");
+			{
+			  prm.declare_entry ("A1", "1085.7",
+							     Patterns::Double (),
+							     "Constant parameter in the quadratic "
+							     "function that approximates the solidus "
+							     "of peridotite. "
+							     "Units: $°C$.");
+			  prm.declare_entry ("A2", "1.329e-7",
+							     Patterns::Double (),
+							     "Prefactor of the linear pressure term "
+							     "in the quadratic function that approximates "
+							     "the solidus of peridotite. "
+							     "Units: $°C/Pa$.");
+			  prm.declare_entry ("A3", "-5.1e-18",
+							     Patterns::Double (),
+							     "Prefactor of the quadratic pressure term "
+							     "in the quadratic function that approximates "
+							     "the solidus of peridotite. "
+							     "Units: $°C/(Pa^2)$.");
+			  prm.declare_entry ("B1", "1475.0",
+							     Patterns::Double (),
+							     "Constant parameter in the quadratic "
+							     "function that approximates the lherzolite "
+							     "liquidus used for calculating the fraction "
+							     "of peridotite-derived melt. "
+							     "Units: $°C$.");
+			  prm.declare_entry ("B2", "8.0e-8",
+							     Patterns::Double (),
+							     "Prefactor of the linear pressure term "
+							     "in the quadratic function that approximates "
+							     "the  lherzolite liquidus used for "
+							     "calculating the fraction of peridotite-"
+							     "derived melt. "
+							     "Units: $°C/Pa$.");
+			  prm.declare_entry ("B3", "-3.2e-18",
+							     Patterns::Double (),
+							     "Prefactor of the quadratic pressure term "
+							     "in the quadratic function that approximates "
+							     "the  lherzolite liquidus used for "
+							     "calculating the fraction of peridotite-"
+							     "derived melt. "
+							     "Units: $°C/(Pa^2)$.");
+			  prm.declare_entry ("C1", "1780.0",
+							     Patterns::Double (),
+							     "Constant parameter in the quadratic "
+							     "function that approximates the liquidus "
+							     "of peridotite. "
+							     "Units: $°C$.");
+			  prm.declare_entry ("C2", "4.50e-8",
+							     Patterns::Double (),
+							     "Prefactor of the linear pressure term "
+							     "in the quadratic function that approximates "
+							     "the liquidus of peridotite. "
+							     "Units: $°C/Pa$.");
+			  prm.declare_entry ("C3", "-2.0e-18",
+							     Patterns::Double (),
+							     "Prefactor of the quadratic pressure term "
+							     "in the quadratic function that approximates "
+							     "the liquidus of peridotite. "
+							     "Units: $°C/(Pa^2)$.");
+			  prm.declare_entry ("r1", "0.4",
+							     Patterns::Double (),
+							     "Constant in the linear function that "
+							     "approximates the clinopyroxene reaction "
+							     "coefficient. "
+							     "Units: non-dimensional.");
+			  prm.declare_entry ("r2", "8e-11",
+							     Patterns::Double (),
+							     "Prefactor of the linear pressure term "
+							     "in the linear function that approximates "
+							     "the clinopyroxene reaction coefficient. "
+							     "Units: $1/Pa$.");
+			  prm.declare_entry ("beta", "1.5",
+							     Patterns::Double (),
+							     "Exponent of the melting temperature in "
+							     "the melt fraction calculation. "
+							     "Units: non-dimensional.");
+			  prm.declare_entry ("M_cpx", "0.3",
+							     Patterns::Double (),
+							     "Mass fraction of clinopyroxene in the "
+							     "peridotite to be molten. "
+							     "Units: non-dimensional.");
+			  prm.declare_entry ("D1", "976.0",
+							     Patterns::Double (),
+							     "Constant parameter in the quadratic "
+							     "function that approximates the solidus "
+							     "of pyroxenite. "
+							     "Units: $°C$.");
+			  prm.declare_entry ("D2", "1.23e-7",
+							     Patterns::Double (),
+							     "Prefactor of the linear pressure term "
+							     "in the quadratic function that approximates "
+							     "the solidus of pyroxenite. "
+							     "Units: $°C/Pa$.");
+			  prm.declare_entry ("D3", "-5.1e-18",
+							     Patterns::Double (),
+							     "Prefactor of the quadratic pressure term "
+							     "in the quadratic function that approximates "
+							     "the solidus of pyroxenite. "
+							     "Units: $°C/(Pa^2)$.");
+			  prm.declare_entry ("E1", "633.8",
+							     Patterns::Double (),
+							     "Prefactor of the linear depletion term "
+							     "in the quadratic function that approximates "
+							     "the melt fraction of pyroxenite. "
+							     "Units: $°C/Pa$.");
+			  prm.declare_entry ("E2", "-611.4",
+							     Patterns::Double (),
+							     "Prefactor of the quadratic depletion term "
+							     "in the quadratic function that approximates "
+							     "the melt fraction of pyroxenite. "
+							     "Units: $°C/(Pa^2)$.");
+	        }
+	        prm.leave_subsection();
+		  }
+		  prm.leave_subsection();
+		}
+		prm.leave_subsection();
+	  }
+
+      template <int dim>
+      void
+      MeltFraction<dim>::parse_parameters (ParameterHandler &prm)
+      {
+  		prm.enter_subsection("Postprocess");
+  		{
+  		  prm.enter_subsection("Visualization");
+  		  {
+  			prm.enter_subsection("Melt fraction");
+  			{
+              A1              = prm.get_double ("A1");
+              A2              = prm.get_double ("A2");
+              A3              = prm.get_double ("A3");
+              B1              = prm.get_double ("B1");
+              B2              = prm.get_double ("B2");
+              B3              = prm.get_double ("B3");
+              C1              = prm.get_double ("C1");
+              C2              = prm.get_double ("C2");
+              C3              = prm.get_double ("C3");
+              r1              = prm.get_double ("r1");
+              r2              = prm.get_double ("r2");
+              beta            = prm.get_double ("beta");
+              M_cpx           = prm.get_double ("M_cpx");
+              D1              = prm.get_double ("D1");
+              D2              = prm.get_double ("D2");
+              D3              = prm.get_double ("D3");
+              E1              = prm.get_double ("E1");
+              E2              = prm.get_double ("E2");
+            }
+            prm.leave_subsection();
+          }
+          prm.leave_subsection();
+        }
+        prm.leave_subsection();
+      }
+    }
+  }
+}
+
+
+// explicit instantiations
+namespace aspect
+{
+  namespace Postprocess
+  {
+    namespace VisualizationPostprocessors
+    {
+      ASPECT_REGISTER_VISUALIZATION_POSTPROCESSOR(MeltFraction,
+                                                  "melt fraction",
+                                                  "A visualization output object that generates output "
+                                                  "for the melt fraction at the temperature and "
+                                                  "pressure of the current point (batch melting). "
+                                                  "Does not take into account latent heat.")
+    }
+  }
+}

Modified: branches/j-dannberg/source/postprocess/visualization.cc
===================================================================
--- branches/j-dannberg/source/postprocess/visualization.cc	2013-11-22 21:35:22 UTC (rev 2039)
+++ branches/j-dannberg/source/postprocess/visualization.cc	2013-11-22 22:06:38 UTC (rev 2040)
@@ -590,6 +590,10 @@
         prm.leave_subsection();
       }
       prm.leave_subsection();
+
+      // now declare the parameters of each of the registered
+      // visualization postprocessors in turn
+      std_cxx1x::get<dim>(registered_plugins).declare_parameters (prm);
     }
 
 
@@ -599,6 +603,7 @@
     {
       Assert (std_cxx1x::get<dim>(registered_plugins).plugins != 0,
               ExcMessage ("No postprocessors registered!?"));
+      std::vector<std::string> viz_names;
 
       prm.enter_subsection("Postprocess");
       {
@@ -609,8 +614,7 @@
           group_files     = prm.get_integer("Number of grouped files");
 
           // now also see which derived quantities we are to compute
-          std::vector<std::string> viz_names
-            = Utilities::split_string_list(prm.get("List of output variables"));
+          viz_names = Utilities::split_string_list(prm.get("List of output variables"));
 
           // see if 'all' was selected (or is part of the list). if so
           // simply replace the list with one that contains all names
@@ -624,37 +628,37 @@
                    p != std_cxx1x::get<dim>(registered_plugins).plugins->end(); ++p)
                 viz_names.push_back (std_cxx1x::get<0>(*p));
             }
-
-          // then go through the list, create objects and let them parse
-          // their own parameters
-          for (unsigned int name=0; name<viz_names.size(); ++name)
-            {
-              VisualizationPostprocessors::Interface<dim> *
-              viz_postprocessor = std_cxx1x::get<dim>(registered_plugins)
-                                  .create_plugin (viz_names[name],
-                                                  "Visualization plugins",
-                                                  prm);
-
-              // make sure that the postprocessor is indeed of type
-              // dealii::DataPostprocessor or of type
-              // VisualizationPostprocessors::CellDataVectorCreator
-              Assert ((dynamic_cast<DataPostprocessor<dim>*>(viz_postprocessor)
-                       != 0)
-                      ||
-                      (dynamic_cast<VisualizationPostprocessors::CellDataVectorCreator<dim>*>(viz_postprocessor)
-                       != 0)
-                      ,
-                      ExcMessage ("Can't convert visualization postprocessor to type "
-                                  "dealii::DataPostprocessor or "
-                                  "VisualizationPostprocessors::CellDataVectorCreator!?"));
-
-              postprocessors.push_back (std_cxx1x::shared_ptr<VisualizationPostprocessors::Interface<dim> >
-                                        (viz_postprocessor));
-            }
         }
         prm.leave_subsection();
       }
       prm.leave_subsection();
+
+	  // then go through the list, create objects and let them parse
+	  // their own parameters
+	  for (unsigned int name=0; name<viz_names.size(); ++name)
+		{
+		  VisualizationPostprocessors::Interface<dim> *
+		  viz_postprocessor = std_cxx1x::get<dim>(registered_plugins)
+							  .create_plugin (viz_names[name],
+											  "Visualization plugins",
+											  prm);
+
+		  // make sure that the postprocessor is indeed of type
+		  // dealii::DataPostprocessor or of type
+		  // VisualizationPostprocessors::CellDataVectorCreator
+		  Assert ((dynamic_cast<DataPostprocessor<dim>*>(viz_postprocessor)
+				   != 0)
+				  ||
+				  (dynamic_cast<VisualizationPostprocessors::CellDataVectorCreator<dim>*>(viz_postprocessor)
+				   != 0)
+				  ,
+				  ExcMessage ("Can't convert visualization postprocessor to type "
+							  "dealii::DataPostprocessor or "
+							  "VisualizationPostprocessors::CellDataVectorCreator!?"));
+
+		  postprocessors.push_back (std_cxx1x::shared_ptr<VisualizationPostprocessors::Interface<dim> >
+									(viz_postprocessor));
+		}
     }
 
 

Modified: branches/j-dannberg/source/simulator/assembly.cc
===================================================================
--- branches/j-dannberg/source/simulator/assembly.cc	2013-11-22 21:35:22 UTC (rev 2039)
+++ branches/j-dannberg/source/simulator/assembly.cc	2013-11-22 22:06:38 UTC (rev 2040)
@@ -384,7 +384,8 @@
         template <int dim>
         struct StokesSystem : public StokesPreconditioner<dim>
         {
-          StokesSystem (const FiniteElement<dim> &finite_element);
+          StokesSystem (const FiniteElement<dim> &finite_element,
+                        const bool                do_pressure_rhs_compatibility_modification);
           StokesSystem (const StokesSystem<dim> &data);
 
           Vector<double> local_rhs;
@@ -395,11 +396,15 @@
 
         template <int dim>


More information about the CIG-COMMITS mailing list