[cig-commits] r1350 - in trunk/aspect: doc/modules include/aspect source/simulator
heien at dealii.org
heien at dealii.org
Tue Nov 6 18:44:43 PST 2012
Author: heien
Date: 2012-11-06 19:44:43 -0700 (Tue, 06 Nov 2012)
New Revision: 1350
Modified:
trunk/aspect/doc/modules/changes.h
trunk/aspect/include/aspect/simulator.h
trunk/aspect/source/simulator/core.cc
trunk/aspect/source/simulator/helper_functions.cc
trunk/aspect/source/simulator/parameters.cc
Log:
Added heat conduction based timestep and parameter to enable it - the current calculation method may not be optimal for the implicit methods used in Aspect
Modified: trunk/aspect/doc/modules/changes.h
===================================================================
--- trunk/aspect/doc/modules/changes.h 2012-11-07 02:37:48 UTC (rev 1349)
+++ trunk/aspect/doc/modules/changes.h 2012-11-07 02:44:43 UTC (rev 1350)
@@ -11,6 +11,14 @@
<ol>
<li>
+New: If the "Use conduction timestep" parameter is true, the timestep
+is calculated as the minimum of the convection *and* heat conduction
+timesteps. This is to allow models where flow is very slow relative
+to conduction.
+<br>
+(Eric Heien, 2012/11/06)
+
+<li>
New: Aspect now catches if the output directory specified in
the parameter file does not exist, rather than providing a
cryptic error message upon first attempt to write to this
Modified: trunk/aspect/include/aspect/simulator.h
===================================================================
--- trunk/aspect/include/aspect/simulator.h 2012-11-07 02:37:48 UTC (rev 1349)
+++ trunk/aspect/include/aspect/simulator.h 2012-11-07 02:44:43 UTC (rev 1350)
@@ -366,7 +366,7 @@
* Return a set of boudary indicators that describes which of the boundaries
* have a fixed temperature.
*/
- const std::set<types::boundary_id_t>&
+ const std::set<types::boundary_id_t> &
get_fixed_temperature_boundary_indicators () const;
@@ -474,6 +474,7 @@
double start_time;
double end_time;
double CFL_number;
+ bool use_conduction_timestep;
bool convert_to_years;
std::string output_directory;
double surface_pressure;
@@ -1170,12 +1171,16 @@
* Compute the size of the next time step from the mesh size and
* the velocity on each cell. The computed time step has to satisfy
* the CFL number chosen in the input parameter file on each cell
- * of the mesh.
+ * of the mesh. If specified in the parameter file, the time step
+ * will be the minimum of the convection *and* conduction time
+ * steps. Also returns whether the timestep is dominated by
+ * convection or conduction.
*
* This function is implemented in
* <code>source/simulator/helper_functions.cc</code>.
*/
- double compute_time_step () const;
+ void compute_time_step (double &new_time_step,
+ bool &convection_dominant) const;
/**
* Compute the artificial diffusion coefficient value on a cell
Modified: trunk/aspect/source/simulator/core.cc
===================================================================
--- trunk/aspect/source/simulator/core.cc 2012-11-07 02:37:48 UTC (rev 1349)
+++ trunk/aspect/source/simulator/core.cc 2012-11-07 02:44:43 UTC (rev 1350)
@@ -196,15 +196,15 @@
// or another in the member initializer list above
// if any plugin wants access to the Simulator by deriving from SimulatorAccess, initialize it:
- if (SimulatorAccess<dim>* sim = dynamic_cast<SimulatorAccess<dim>*>(&*geometry_model))
+ if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(&*geometry_model))
sim->initialize (*this);
- if (SimulatorAccess<dim>* sim = dynamic_cast<SimulatorAccess<dim>*>(&*material_model))
+ if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(&*material_model))
sim->initialize (*this);
- if (SimulatorAccess<dim>* sim = dynamic_cast<SimulatorAccess<dim>*>(&*gravity_model))
+ if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(&*gravity_model))
sim->initialize (*this);
- if (SimulatorAccess<dim>* sim = dynamic_cast<SimulatorAccess<dim>*>(&*boundary_temperature))
+ if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(&*boundary_temperature))
sim->initialize (*this);
- if (SimulatorAccess<dim>* sim = dynamic_cast<SimulatorAccess<dim>*>(&*initial_conditions))
+ if (SimulatorAccess<dim> *sim = dynamic_cast<SimulatorAccess<dim>*>(&*initial_conditions))
sim->initialize (*this);
postprocess_manager.parse_parameters (prm);
@@ -1336,8 +1336,9 @@
pcout << std::endl;
// update the time step size
+ bool convection_dominant; // for now this is unused, will be added to statistics later
old_time_step = time_step;
- time_step = compute_time_step();
+ compute_time_step(time_step, convection_dominant);
if (parameters.convert_to_years == true)
statistics.add_value("Time step size (years)", time_step / year_in_seconds);
Modified: trunk/aspect/source/simulator/helper_functions.cc
===================================================================
--- trunk/aspect/source/simulator/helper_functions.cc 2012-11-07 02:37:48 UTC (rev 1349)
+++ trunk/aspect/source/simulator/helper_functions.cc 2012-11-07 02:44:43 UTC (rev 1350)
@@ -207,18 +207,22 @@
template <int dim>
- double Simulator<dim>::compute_time_step () const
+ void Simulator<dim>::compute_time_step (double &new_time_step, bool &convection_dominant) const
{
const QIterated<dim> quadrature_formula (QTrapez<1>(),
parameters.stokes_velocity_degree);
const unsigned int n_q_points = quadrature_formula.size();
- FEValues<dim> fe_values (mapping, finite_element, quadrature_formula, update_values);
+ FEValues<dim> fe_values (mapping, finite_element, quadrature_formula, update_values | (parameters.use_conduction_timestep ? update_quadrature_points : update_default));
std::vector<Tensor<1,dim> > velocity_values(n_q_points);
+ std::vector<double> pressure_values(n_q_points), temperature_values(n_q_points);
const FEValuesExtractors::Vector velocities (0);
+ const FEValuesExtractors::Scalar pressure (dim);
+ const FEValuesExtractors::Scalar temperature (dim+1);
double max_local_speed_over_meshsize = 0;
+ double min_local_conduction_timestep = std::numeric_limits<double>::max(), min_conduction_timestep;
typename DoFHandler<dim>::active_cell_iterator
cell = dof_handler.begin_active(),
@@ -238,19 +242,53 @@
max_local_velocity
/
cell->minimum_vertex_distance());
+ if (parameters.use_conduction_timestep)
+ {
+ fe_values[pressure].get_function_values (solution,
+ pressure_values);
+ fe_values[temperature].get_function_values (solution,
+ temperature_values);
+ // In the future we may want to evaluate thermal diffusivity at
+ // each point in the mesh, but for now we just use the reference value
+ for (unsigned int q=0; q<n_q_points; ++q)
+ {
+ std::vector<double> composition_values_at_q_point (parameters.n_compositional_fields);
+ double thermal_diffusivity;
+
+ // TODO: calculate composition field as well
+ thermal_diffusivity = material_model->thermal_diffusivity(temperature_values[q],
+ pressure_values[q],
+ composition_values_at_q_point,
+ fe_values.quadrature_point(q));
+
+ min_local_conduction_timestep = std::min(min_local_conduction_timestep,
+ pow(cell->minimum_vertex_distance(),2)
+ / thermal_diffusivity);
+ }
+ }
}
const double max_global_speed_over_meshsize
= Utilities::MPI::max (max_local_speed_over_meshsize, mpi_communicator);
+ if (parameters.use_conduction_timestep)
+ MPI_Allreduce (&min_local_conduction_timestep, &min_conduction_timestep, 1, MPI_DOUBLE, MPI_MIN, mpi_communicator);
+ else
+ min_conduction_timestep = std::numeric_limits<double>::max();
// if the velocity is zero, then it is somewhat arbitrary what time step
// we should choose. in that case, do as if the velocity was one
- if (max_global_speed_over_meshsize == 0)
- return (parameters.CFL_number / (parameters.temperature_degree *
- 1));
+ if (max_global_speed_over_meshsize == 0 && !parameters.use_conduction_timestep)
+ {
+ new_time_step = (parameters.CFL_number /
+ (parameters.temperature_degree * 1));
+ convection_dominant = false;
+ }
else
- return (parameters.CFL_number / (parameters.temperature_degree *
- max_global_speed_over_meshsize));
+ {
+ new_time_step = std::min(min_conduction_timestep,
+ (parameters.CFL_number / (parameters.temperature_degree * max_global_speed_over_meshsize)));
+ convection_dominant = (new_time_step < min_conduction_timestep);
+ }
}
@@ -1121,7 +1159,7 @@
template void Simulator<dim>::extract_composition_values_at_q_point (const std::vector<std::vector<double> > &composition_values, \
const unsigned int q, \
std::vector<double> &composition_values_at_q_point) const; \
- template double Simulator<dim>::compute_time_step () const; \
+ template void Simulator<dim>::compute_time_step (double &new_time_step, bool &convection_dominant) const; \
template void Simulator<dim>::make_pressure_rhs_compatible(LinearAlgebra::BlockVector &vector); \
template void Simulator<dim>::compute_depth_average_field(const unsigned int index, std::vector<double> &values) const; \
template void Simulator<dim>::compute_depth_average_viscosity(std::vector<double> &values) const; \
Modified: trunk/aspect/source/simulator/parameters.cc
===================================================================
--- trunk/aspect/source/simulator/parameters.cc 2012-11-07 02:37:48 UTC (rev 1349)
+++ trunk/aspect/source/simulator/parameters.cc 2012-11-07 02:44:43 UTC (rev 1350)
@@ -99,6 +99,14 @@
"one can choose $c>1$) though a CFL number significantly larger than "
"one will yield rather diffusive solutions. Units: None.");
+ prm.declare_entry ("Use conduction timestep", "false",
+ Patterns::Bool (),
+ "Mantle convection simulations are often focused on convection "
+ "dominated systems. However, these codes can also be used to "
+ "investigate systems where heat conduction plays a dominant role. "
+ "This parameter indicates whether the simulator should also use "
+ "heat conduction in determining the length of each time step.");
+
prm.declare_entry ("Nonlinear solver scheme", "IMPES",
Patterns::Selection ("IMPES|iterated IMPES|iterated Stokes"),
"The kind of scheme used to resolve the nonlinearity in the system. "
@@ -442,6 +450,7 @@
resume_computation = prm.get_bool ("Resume computation");
CFL_number = prm.get_double ("CFL number");
+ use_conduction_timestep = prm.get_bool ("Use conduction timestep");
convert_to_years = prm.get_bool ("Use years in output instead of seconds");
timing_output_frequency = prm.get_integer ("Timing output frequency");
More information about the CIG-COMMITS
mailing list