[cig-commits] commit 1937 by buerg to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Fri Oct 4 10:21:49 PDT 2013
Revision 1937
Add Strang splitting.
U trunk/aspire/include/aspect/material_model/cantera.h
U trunk/aspire/include/aspect/material_model/interface.h
U trunk/aspire/include/aspect/material_model/mixture_fraction.h
U trunk/aspire/include/aspect/material_model/simple.h
U trunk/aspire/include/aspect/material_model/three_species.h
U trunk/aspire/include/aspect/material_model/three_species_simple.h
U trunk/aspire/source/material_model/cantera.cc
U trunk/aspire/source/material_model/mixture_fraction.cc
U trunk/aspire/source/material_model/simple.cc
U trunk/aspire/source/material_model/three_species.cc
U trunk/aspire/source/material_model/three_species_simple.cc
U trunk/aspire/source/simulator/assembly.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1937&peg=1937
Diff:
Modified: trunk/aspire/include/aspect/material_model/cantera.h
===================================================================
--- trunk/aspire/include/aspect/material_model/cantera.h 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/include/aspect/material_model/cantera.h 2013-10-04 17:21:10 UTC (rev 1937)
@@ -109,6 +109,7 @@
virtual void evaluate (const MaterialModelInputs &in, MaterialModelOutputs &out) const;
+ virtual void react (const MaterialModelInputs &in, MaterialModelOutputs &out) const;
/**
Modified: trunk/aspire/include/aspect/material_model/interface.h
===================================================================
--- trunk/aspire/include/aspect/material_model/interface.h 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/include/aspect/material_model/interface.h 2013-10-04 17:21:10 UTC (rev 1937)
@@ -484,6 +484,11 @@
struct MaterialModelInputs
{
MaterialModelInputs(unsigned int n_points, unsigned int n_comp);
+
+ /**
+ * Current time_step of the simulation.
+ */
+ double time_step;
/**
* Vector with global positions where the material has to be evaluated in evaluate().
@@ -574,9 +579,18 @@
* @param in
* @param out
*/
- virtual void evaluate(const MaterialModelInputs &in, MaterialModelOutputs &out) const = 0;
+ virtual void evaluate (const MaterialModelInputs &in, MaterialModelOutputs &out) const = 0;
/**
+ * The react() function is implemented to call the individual
+ * functions in this class, so there is no need to implement this
+ * in your material model derived from InterfaceCompatibility.
+ * @param in
+ * @param out
+ */
+ virtual void react (const MaterialModelInputs &in, MaterialModelOutputs &out) const = 0;
+
+ /**
* @name Functions used in dealing with run-time parameters
* @{
*/
Modified: trunk/aspire/include/aspect/material_model/mixture_fraction.h
===================================================================
--- trunk/aspire/include/aspect/material_model/mixture_fraction.h 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/include/aspect/material_model/mixture_fraction.h 2013-10-04 17:21:10 UTC (rev 1937)
@@ -106,7 +106,8 @@
typedef typename MaterialModel::Interface<dim>::MaterialModelOutputs MaterialModelOutputs;
- virtual void evaluate(const MaterialModelInputs &in, MaterialModelOutputs &out) const;
+ virtual void evaluate (const MaterialModelInputs &in, MaterialModelOutputs &out) const;
+ virtual void react (const MaterialModelInputs &in, MaterialModelOutputs &out) const;
/**
Modified: trunk/aspire/include/aspect/material_model/simple.h
===================================================================
--- trunk/aspire/include/aspect/material_model/simple.h 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/include/aspect/material_model/simple.h 2013-10-04 17:21:10 UTC (rev 1937)
@@ -107,7 +107,8 @@
typedef typename MaterialModel::Interface<dim>::MaterialModelOutputs MaterialModelOutputs;
- virtual void evaluate(const MaterialModelInputs &in, MaterialModelOutputs &out) const;
+ virtual void evaluate (const MaterialModelInputs &in, MaterialModelOutputs &out) const;
+ virtual void react (const MaterialModelInputs &in, MaterialModelOutputs &out) const;
/**
Modified: trunk/aspire/include/aspect/material_model/three_species.h
===================================================================
--- trunk/aspire/include/aspect/material_model/three_species.h 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/include/aspect/material_model/three_species.h 2013-10-04 17:21:10 UTC (rev 1937)
@@ -102,6 +102,7 @@
virtual void evaluate (const MaterialModelInputs &in, MaterialModelOutputs &out) const;
+ virtual void react (const MaterialModelInputs &in, MaterialModelOutputs &out) const;
/**
Modified: trunk/aspire/include/aspect/material_model/three_species_simple.h
===================================================================
--- trunk/aspire/include/aspect/material_model/three_species_simple.h 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/include/aspect/material_model/three_species_simple.h 2013-10-04 17:21:10 UTC (rev 1937)
@@ -102,6 +102,7 @@
virtual void evaluate (const MaterialModelInputs &in, MaterialModelOutputs &out) const;
+ virtual void react (const MaterialModelInputs &in, MaterialModelOutputs &out) const;
/**
Modified: trunk/aspire/source/material_model/cantera.cc
===================================================================
--- trunk/aspire/source/material_model/cantera.cc 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/source/material_model/cantera.cc 2013-10-04 17:21:10 UTC (rev 1937)
@@ -39,7 +39,7 @@
template <int dim>
void
- Cantera<dim>::evaluate(const MaterialModelInputs &in, MaterialModelOutputs &out) const
+ Cantera<dim>::evaluate (const MaterialModelInputs &in, MaterialModelOutputs &out) const
{
Assert (in.composition[0].size () >= 1, ExcMessage ("Need compositional field!"));
@@ -53,23 +53,37 @@
out.specific_heat[q] = ideal_gas_mix->cp_mass ();
out.thermal_conductivities[q] = transport->thermalConductivity ();
out.diffusivities[q] = out.thermal_conductivities[q] / (out.densities[q] * out.specific_heat[q]);
+ out.thermal_expansion_coefficients[q] = ideal_gas_mix->thermalExpansionCoeff ();
+ out.viscosities[q] = transport->viscosity ();
+ }
+ }
+
+ template <int dim>
+ void
+ Cantera<dim>::react (const MaterialModelInputs &in, MaterialModelOutputs &out) const
+ {
+ Assert (in.composition[0].size () >= 1, ExcMessage ("Need compositional field!"));
+
+ const unsigned int n_species = ideal_gas_mix->nSpecies ();
+
+ for (unsigned int q = 0; q < out.viscosities.size (); ++q)
+ {
+ ideal_gas_mix->setState_TPY (in.temperature[q], reference_pressure, &in.composition[q][0]);
+ out.densities[q] = ideal_gas_mix->density ();
ideal_gas_mix->getEnthalpy_RT (enthalpy_RT);
ideal_gas_mix->getNetProductionRates (net);
out.thermal_sources[q] = 0.0;
for (unsigned int i = 0; i < n_species; ++i)
{
- out.compositional_sources[q][i] = molar_mass[i] * net[i] / out.densities[q];
- out.thermal_sources[q] += enthalpy_RT[i] * out.compositional_sources[q][i];
+ out.compositional_sources[q][i] = in.time_step * molar_mass[i] * net[i] / out.densities[q];
+ out.thermal_sources[q] += out.compositional_sources[q][i];
}
- out.thermal_sources[q] *= -1.0 * ideal_gas_mix->enthalpy_mass ();
- out.thermal_expansion_coefficients[q] = ideal_gas_mix->thermalExpansionCoeff ();
- out.viscosities[q] = transport->viscosity ();
+ out.thermal_sources[q] *= -1.0 * ideal_gas_mix->enthalpy_mass () / ideal_gas_mix->entropy_mass ();
}
}
-
template <int dim>
bool
Cantera<dim>::is_compressible () const
Modified: trunk/aspire/source/material_model/mixture_fraction.cc
===================================================================
--- trunk/aspire/source/material_model/mixture_fraction.cc 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/source/material_model/mixture_fraction.cc 2013-10-04 17:21:10 UTC (rev 1937)
@@ -61,6 +61,14 @@
template <int dim>
+ void
+ MixtureFraction<dim>::react (const MaterialModelInputs &in, MaterialModelOutputs &out) const
+ {
+ Assert(in.composition[0].size()>=1, ExcMessage("need compositional field!"));
+ }
+
+
+ template <int dim>
bool
MixtureFraction<dim>::is_compressible () const
{
Modified: trunk/aspire/source/material_model/simple.cc
===================================================================
--- trunk/aspire/source/material_model/simple.cc 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/source/material_model/simple.cc 2013-10-04 17:21:10 UTC (rev 1937)
@@ -48,6 +48,12 @@
}
}
+ template <int dim>
+ void
+ Simple<dim>::react (const MaterialModelInputs &in, MaterialModelOutputs &out) const
+ {
+ Assert(in.composition[0].size()>=1, ExcMessage("need compositional field!"));
+ }
template <int dim>
bool
Modified: trunk/aspire/source/material_model/three_species.cc
===================================================================
--- trunk/aspire/source/material_model/three_species.cc 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/source/material_model/three_species.cc 2013-10-04 17:21:10 UTC (rev 1937)
@@ -31,7 +31,7 @@
{
template <int dim>
void
- ThreeSpecies<dim>::evaluate(const MaterialModelInputs &in, MaterialModelOutputs &out) const
+ ThreeSpecies<dim>::evaluate (const MaterialModelInputs &in, MaterialModelOutputs &out) const
{
Assert (in.composition[0].size () >= 1, ExcMessage ("Need compositional field!"));
@@ -41,14 +41,24 @@
out.specific_heat[q] = 1.0;
out.thermal_conductivities[q] = 0.0;
out.diffusivities[q] = out.thermal_conductivities[q] / (out.densities[q] * out.specific_heat[q]);
- out.thermal_sources[q] = 1e2 * in.composition[q][0] * in.composition[q][1];
- out.compositional_sources[q][0] = -1.0 * in.composition[q][0] * in.composition[q][1];
- out.compositional_sources[q][1] = out.compositional_sources[q][0];
- out.compositional_sources[q][2] = in.composition[q][0] * in.composition[q][1];
out.viscosities[q] = 1e-3;
}
}
+ template <int dim>
+ void
+ ThreeSpecies<dim>::react (const MaterialModelInputs &in, MaterialModelOutputs &out) const
+ {
+ Assert (in.composition[0].size () >= 1, ExcMessage ("Need compositional field!"));
+
+ for (unsigned int q = 0; q < out.viscosities.size (); ++q)
+ {
+ out.thermal_sources[q] = 1e2 * in.composition[q][0] * in.composition[q][1] * in.time_step;
+ out.compositional_sources[q][0] = in.composition[q][0] * (std::exp (-1.0 * in.composition[q][1] * in.time_step) - 1.0);
+ out.compositional_sources[q][1] = in.composition[q][1] * (std::exp (-1.0 * in.composition[q][0] * in.time_step) - 1.0);
+ out.compositional_sources[q][2] = in.composition[q][0] * in.composition[q][1] * in.time_step;
+ }
+ }
template <int dim>
bool
Modified: trunk/aspire/source/material_model/three_species_simple.cc
===================================================================
--- trunk/aspire/source/material_model/three_species_simple.cc 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/source/material_model/three_species_simple.cc 2013-10-04 17:21:10 UTC (rev 1937)
@@ -31,7 +31,7 @@
{
template <int dim>
void
- ThreeSpeciesSimple<dim>::evaluate(const MaterialModelInputs &in, MaterialModelOutputs &out) const
+ ThreeSpeciesSimple<dim>::evaluate (const MaterialModelInputs &in, MaterialModelOutputs &out) const
{
Assert (in.composition[0].size () >= 1, ExcMessage ("Need compositional field!"));
@@ -41,14 +41,24 @@
out.specific_heat[q] = 1.0;
out.thermal_conductivities[q] = 0.0;
out.diffusivities[q] = out.thermal_conductivities[q] / (out.densities[q] * out.specific_heat[q]);
- out.thermal_sources[q] = 0.0;
- out.compositional_sources[q][0] = -1.0 * in.composition[q][0] * in.composition[q][1];
- out.compositional_sources[q][1] = out.compositional_sources[q][0];
- out.compositional_sources[q][2] = in.composition[q][0] * in.composition[q][1];
out.viscosities[q] = 1e-3;
}
}
+ template <int dim>
+ void
+ ThreeSpeciesSimple<dim>::react (const MaterialModelInputs &in, MaterialModelOutputs &out) const
+ {
+ Assert (in.composition[0].size () >= 1, ExcMessage ("Need compositional field!"));
+
+ for (unsigned int q = 0; q < out.viscosities.size (); ++q)
+ {
+ out.thermal_sources[q] = 0.0;
+ out.compositional_sources[q][0] = in.composition[q][0] * (std::exp (-1.0 * in.composition[q][1] * in.time_step) - 1.0);
+ out.compositional_sources[q][1] = in.composition[q][1] * (std::exp (-1.0 * in.composition[q][0] * in.time_step) - 1.0);
+ out.compositional_sources[q][2] = in.composition[q][0] * in.composition[q][1] * in.time_step;
+ }
+ }
template <int dim>
bool
Modified: trunk/aspire/source/simulator/assembly.cc
===================================================================
--- trunk/aspire/source/simulator/assembly.cc 2013-10-03 15:12:46 UTC (rev 1936)
+++ trunk/aspire/source/simulator/assembly.cc 2013-10-04 17:21:10 UTC (rev 1937)
@@ -698,6 +698,8 @@
for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
material_model_inputs.composition[q][c] = input_solution_composition.block (c) (dof);
}
+
+ material_model_inputs.time_step = 0.5 * time_step;
}
@@ -1001,7 +1003,6 @@
const double lambda = scratch.material_model_outputs.thermal_conductivities[q];
const double old_value = scratch.old_temperature_values[q];
const double rho = scratch.material_model_outputs.densities[q];
- const double source_term = scratch.material_model_outputs.thermal_sources[q] / c_P;
const Tensor<1, dim> current_u = scratch.current_velocity_values[q];
for (unsigned int i = 0; i < dofs_per_cell; ++i)
@@ -1017,7 +1018,7 @@
scratch.fe_values_advection.quadrature_point (q) (0) : 1.0);
}
- data.local_rhs (i) += scratch.phi_field[i] * (time_step * source_term + old_value) * scratch.fe_values_advection.JxW (q)
+ data.local_rhs (i) += scratch.phi_field[i] * old_value * scratch.fe_values_advection.JxW (q)
* (axisymmetric_terms ? scratch.fe_values_advection.quadrature_point (q) (0) : 1.0);
}
}
@@ -1076,7 +1077,6 @@
{
const double diffusivity = rho * scratch.material_model_outputs.diffusivities[q];
const double old_value = scratch.old_composition_values[c][q];
- const double source_term = scratch.material_model_outputs.compositional_sources[q][c];
for (unsigned int i = 0; i < dofs_per_cell; ++i)
{
@@ -1091,7 +1091,7 @@
scratch.fe_values_advection.quadrature_point (q) (0) : 1.0);
}
- data.local_rhs.block (c) (i) += scratch.phi_field[i] * (time_step * source_term + old_value) * scratch.fe_values_advection.JxW (q)
+ data.local_rhs.block (c) (i) += scratch.phi_field[i] * old_value * scratch.fe_values_advection.JxW (q)
* (axisymmetric_terms ? scratch.fe_values_advection.quadrature_point (q) (0) : 1.0);
}
}
@@ -1263,6 +1263,8 @@
Simulator<dim>::advance_chemistry ()
{
computing_timer.enter_section (" Perform chemical reactions");
+ solution_temperature = old_solution_temperature;
+ solution_composition = old_solution_composition;
const IndexSet& locally_owned_dofs = dof_handler_advection.locally_owned_dofs ();
typename MaterialModel::Interface<dim>::MaterialModelInputs material_model_inputs (locally_owned_dofs.n_elements (),
@@ -1273,6 +1275,30 @@
locally_owned_dofs,
material_model_inputs);
+ typename MaterialModel::Interface<dim>::MaterialModelOutputs material_model_outputs (locally_owned_dofs.n_elements (),
+ parameters.n_compositional_fields);
+
+ material_model->react (material_model_inputs, material_model_outputs);
+
+ for (unsigned int i = 0; i < locally_owned_dofs.n_elements (); ++i)
+ {
+ const types::global_dof_index dof = locally_owned_dofs.nth_index_in_set (i);
+
+ solution_temperature (dof) += material_model_outputs.thermal_sources[i];
+
+ for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
+ solution_composition.block (c) (dof) += material_model_outputs.compositional_sources[i][c];
+ }
+
+ solution_temperature.compress (VectorOperation::add);
+ constraints_temperature.distribute (solution_temperature);
+ current_linearization_point_temperature = solution_temperature;
+ solution_composition.compress (VectorOperation::add);
+
+ for (unsigned int c = 0; c < parameters.n_compositional_fields; ++c)
+ current_constraints_composition[c].distribute (solution_composition.block (c));
+
+ current_linearization_point_composition = solution_composition;
computing_timer.exit_section ();
}
}
More information about the CIG-COMMITS
mailing list