[cig-commits] [commit] master: Evaluate time dependent boundary conditions in every time step, in just the same way as we do for the velocity boundary conditions. (10c017c)

cig_noreply at geodynamics.org cig_noreply at geodynamics.org
Wed May 21 14:30:01 PDT 2014


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

On branch  : master
Link       : https://github.com/geodynamics/aspect/compare/58a3ca45d49c7d331fa02d1fce0b2a74705241eb...7c9b74ed7f6e149b86519d3abf0c555383cba169

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

commit 10c017cd2260ed4cec4e62a6797d1585b4b920e9
Author: Wolfgang Bangerth <bangerth at math.tamu.edu>
Date:   Wed May 21 15:11:41 2014 -0500

    Evaluate time dependent boundary conditions in every time step, in just the same way as we do for the velocity boundary conditions.


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

10c017cd2260ed4cec4e62a6797d1585b4b920e9
 source/simulator/core.cc                           | 118 ++++++++++-----------
 tests/time-dependent-temperature-bc-2.cc           |   1 +
 ...-bc.prm => time-dependent-temperature-bc-2.prm} |  11 +-
 .../screen-output                                  |  36 ++++---
 4 files changed, 83 insertions(+), 83 deletions(-)

diff --git a/source/simulator/core.cc b/source/simulator/core.cc
index 2de57d4..9288c2c 100644
--- a/source/simulator/core.cc
+++ b/source/simulator/core.cc
@@ -505,12 +505,11 @@ namespace aspect
           * introspection.system_dofs_per_block[introspection.block_indices.compositional_fields[0]]);
 
     // then interpolate the current boundary velocities. this adds to
-    // the current_constraints object we already have
+    // the 'constraints' object we already have
+    current_constraints.clear ();
+    current_constraints.reinit (introspection.index_sets.system_relevant_set);
+    current_constraints.merge (constraints);
     {
-      current_constraints.clear ();
-      current_constraints.reinit (introspection.index_sets.system_relevant_set);
-      current_constraints.merge (constraints);
-
       // set the current time and do the interpolation
       // for the prescribed velocity fields
       for (typename std::map<types::boundary_id,std_cxx1x::shared_ptr<VelocityBoundaryConditions::Interface<dim> > >::iterator
@@ -554,9 +553,61 @@ namespace aspect
                                                     current_constraints,
                                                     mask);
         }
-      current_constraints.close();
     }
 
+    // now do the same for the temperature variable
+    {
+      // obtain the boundary indicators that belong to Dirichlet-type
+      // temperature boundary conditions and interpolate the temperature
+      // there
+      for (std::set<types::boundary_id>::const_iterator
+           p = parameters.fixed_temperature_boundary_indicators.begin();
+           p != parameters.fixed_temperature_boundary_indicators.end(); ++p)
+        {
+          Assert (is_element (*p, geometry_model->get_used_boundary_indicators()),
+                  ExcInternalError());
+          VectorTools::interpolate_boundary_values (dof_handler,
+                                                    *p,
+                                                    VectorFunctionFromScalarFunctionObject<dim>(std_cxx1x::bind (&BoundaryTemperature::Interface<dim>::temperature,
+                                                        std_cxx1x::cref(*boundary_temperature),
+                                                        std_cxx1x::cref(*geometry_model),
+                                                        *p,
+                                                        std_cxx1x::_1),
+                                                        introspection.component_masks.temperature.first_selected_component(),
+                                                        introspection.n_components),
+                                                    current_constraints,
+                                                    introspection.component_masks.temperature);
+
+        }
+
+      // now do the same for the composition variable:
+      // obtain the boundary indicators that belong to Dirichlet-type
+      // composition boundary conditions and interpolate the composition
+      // there
+      for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
+        for (std::set<types::boundary_id>::const_iterator
+             p = parameters.fixed_composition_boundary_indicators.begin();
+             p != parameters.fixed_composition_boundary_indicators.end(); ++p)
+          {
+            Assert (is_element (*p, geometry_model->get_used_boundary_indicators()),
+                    ExcInternalError());
+            VectorTools::interpolate_boundary_values (dof_handler,
+                                                      *p,
+                                                      VectorFunctionFromScalarFunctionObject<dim>(std_cxx1x::bind (&BoundaryComposition::Interface<dim>::composition,
+                                                          std_cxx1x::cref(*boundary_composition),
+                                                          std_cxx1x::cref(*geometry_model),
+                                                          *p,
+                                                          std_cxx1x::_1,
+                                                          c),
+                                                          introspection.component_masks.compositional_fields[c].first_selected_component(),
+                                                          introspection.n_components),
+                                                      current_constraints,
+                                                      introspection.component_masks.compositional_fields[c]);
+          }
+    }
+    current_constraints.close();
+
+
     //TODO: do this in a more efficient way (TH)? we really only need
     // to make sure that the time dependent velocity boundary conditions
     // end up in the right hand side in the right way; we currently do
@@ -564,6 +615,8 @@ namespace aspect
     if (!velocity_boundary_conditions.empty())
       rebuild_stokes_matrix = rebuild_stokes_preconditioner = true;
 
+
+
     // notify different system components that we started the next time step
     material_model->update();
     gravity_model->update();
@@ -787,8 +840,6 @@ namespace aspect
     // velocity that are different between time steps. these are then put
     // into current_constraints in start_timestep().
     {
-
-
       // do the interpolation for zero velocity
       for (std::set<types::boundary_id>::const_iterator
            p = parameters.zero_velocity_boundary_indicators.begin();
@@ -808,57 +859,6 @@ namespace aspect
                                                        constraints,
                                                        mapping);
     }
-
-    // now do the same for the temperature variable
-    {
-      // obtain the boundary indicators that belong to Dirichlet-type
-      // temperature boundary conditions and interpolate the temperature
-      // there
-      for (std::set<types::boundary_id>::const_iterator
-           p = parameters.fixed_temperature_boundary_indicators.begin();
-           p != parameters.fixed_temperature_boundary_indicators.end(); ++p)
-        {
-          Assert (is_element (*p, geometry_model->get_used_boundary_indicators()),
-                  ExcInternalError());
-          VectorTools::interpolate_boundary_values (dof_handler,
-                                                    *p,
-                                                    VectorFunctionFromScalarFunctionObject<dim>(std_cxx1x::bind (&BoundaryTemperature::Interface<dim>::temperature,
-                                                        std_cxx1x::cref(*boundary_temperature),
-                                                        std_cxx1x::cref(*geometry_model),
-                                                        *p,
-                                                        std_cxx1x::_1),
-                                                        introspection.component_masks.temperature.first_selected_component(),
-                                                        introspection.n_components),
-                                                    constraints,
-                                                    introspection.component_masks.temperature);
-
-        }
-
-      // now do the same for the composition variable:
-      // obtain the boundary indicators that belong to Dirichlet-type
-      // composition boundary conditions and interpolate the composition
-      // there
-      for (unsigned int c=0; c<parameters.n_compositional_fields; ++c)
-        for (std::set<types::boundary_id>::const_iterator
-             p = parameters.fixed_composition_boundary_indicators.begin();
-             p != parameters.fixed_composition_boundary_indicators.end(); ++p)
-          {
-            Assert (is_element (*p, geometry_model->get_used_boundary_indicators()),
-                    ExcInternalError());
-            VectorTools::interpolate_boundary_values (dof_handler,
-                                                      *p,
-                                                      VectorFunctionFromScalarFunctionObject<dim>(std_cxx1x::bind (&BoundaryComposition::Interface<dim>::composition,
-                                                          std_cxx1x::cref(*boundary_composition),
-                                                          std_cxx1x::cref(*geometry_model),
-                                                          *p,
-                                                          std_cxx1x::_1,
-                                                          c),
-                                                          introspection.component_masks.compositional_fields[c].first_selected_component(),
-                                                          introspection.n_components),
-                                                      constraints,
-                                                      introspection.component_masks.compositional_fields[c]);
-          }
-    }
     constraints.close();
 
     // finally initialize vectors, matrices, etc.
diff --git a/tests/time-dependent-temperature-bc-2.cc b/tests/time-dependent-temperature-bc-2.cc
new file mode 100644
index 0000000..44d5456
--- /dev/null
+++ b/tests/time-dependent-temperature-bc-2.cc
@@ -0,0 +1 @@
+#include "time-dependent-temperature-bc.cc"
diff --git a/tests/time-dependent-temperature-bc.prm b/tests/time-dependent-temperature-bc-2.prm
similarity index 87%
copy from tests/time-dependent-temperature-bc.prm
copy to tests/time-dependent-temperature-bc-2.prm
index 208261f..2d6b7c8 100644
--- a/tests/time-dependent-temperature-bc.prm
+++ b/tests/time-dependent-temperature-bc-2.prm
@@ -1,19 +1,16 @@
 # originally taken from diffusion.prm, but intended to test that we can deal
 # with time dependent temperature boundary conditions
 #
-# this is just a preliminary test. we used to crash because we forgot to set
-# the time before we evaluate the temperature boundary conditions for the
-# first time, and plugins providing such boundary conditions consequently
-# had no way to tell the time when evaluated when computing the initial
-# solution
+# we use a time dependent boundary condition at the top boundary and verify
+# that the temperature at the indeed changes with time as expected
 
 
 set Dimension = 2
 
 
-set CFL number                             = 0.01
+set CFL number                             = 1
 
-set End time                               = 0
+set End time                               = 1
 
 
 set Resume computation                     = false
diff --git a/tests/zero-temperature/screen-output b/tests/time-dependent-temperature-bc-2/screen-output
similarity index 50%
copy from tests/zero-temperature/screen-output
copy to tests/time-dependent-temperature-bc-2/screen-output
index 3c8947a..b617721 100644
--- a/tests/zero-temperature/screen-output
+++ b/tests/time-dependent-temperature-bc-2/screen-output
@@ -1,12 +1,14 @@
 -----------------------------------------------------------------------------
 -- This is ASPECT, the Advanced Solver for Problems in Earth's ConvecTion.
---     . running in OPTIMIZED mode
+--     . running in DEBUG mode
 --     . running with 1 MPI process
 --     . using Trilinos
 -----------------------------------------------------------------------------
 
-Number of active cells: 1,024 (on 6 levels)
-Number of degrees of freedom: 13,764 (8,450+1,089+4,225)
+Loading shared library <./libtime-dependent-temperature-bc-2.so>
+
+Number of active cells: 16 (on 3 levels)
+Number of degrees of freedom: 268 (162+25+81)
 
 *** Timestep 0:  t=0 seconds
    Solving temperature system... 0 iterations.
@@ -17,35 +19,35 @@ Number of degrees of freedom: 13,764 (8,450+1,089+4,225)
      Temperature min/avg/max: 0 K, 0 K, 0 K
 
 *** Timestep 1:  t=0.5 seconds
-   Solving temperature system... 0 iterations.
+   Solving temperature system... 11 iterations.
    Solving Stokes system... 0 iterations.
 
    Postprocessing:
-     Temperature min/avg/max: 0 K, 0 K, 0 K
+     Temperature min/avg/max: -0.07318 K, 0.01471 K, 0.5 K
 
 *** Timestep 2:  t=1 seconds
-   Solving temperature system... 0 iterations.
+   Solving temperature system... 8 iterations.
    Solving Stokes system... 0 iterations.
 
    Postprocessing:
-     Temperature min/avg/max: 0 K, 0 K, 0 K
+     Temperature min/avg/max: -0.1463 K, 0.02943 K, 1 K
 
 Termination requested by criterion: end time
 
 
 +---------------------------------------------+------------+------------+
-| Total wallclock time elapsed since start    |     0.458s |            |
+| Total wallclock time elapsed since start    |     0.247s |            |
 |                                             |            |            |
 | Section                         | no. calls |  wall time | % of total |
 +---------------------------------+-----------+------------+------------+
-| Assemble Stokes system          |         3 |    0.0874s |        19% |
-| Assemble temperature system     |         3 |     0.129s |        28% |
-| Build Stokes preconditioner     |         1 |    0.0779s |        17% |
-| Build temperature preconditioner|         3 |    0.0444s |       9.7% |
-| Solve Stokes system             |         3 |    0.0101s |       2.2% |
-| Solve temperature system        |         3 |   0.00318s |       0.7% |
-| Initialization                  |         2 |    0.0216s |       4.7% |
-| Postprocessing                  |         3 |   0.00506s |       1.1% |
-| Setup dof systems               |         1 |    0.0644s |        14% |
+| Assemble Stokes system          |         3 |    0.0309s |        13% |
+| Assemble temperature system     |         3 |    0.0823s |        33% |
+| Build Stokes preconditioner     |         1 |     0.024s |       9.7% |
+| Build temperature preconditioner|         3 |   0.00207s |      0.84% |
+| Solve Stokes system             |         3 |   0.00508s |       2.1% |
+| Solve temperature system        |         3 |   0.00208s |      0.84% |
+| Initialization                  |         2 |    0.0404s |        16% |
+| Postprocessing                  |         3 |   0.00616s |       2.5% |
+| Setup dof systems               |         1 |    0.0234s |       9.5% |
 +---------------------------------+-----------+------------+------------+
 



More information about the CIG-COMMITS mailing list