[cig-commits] commit 1884 by heister to /var/svn/dealii/aspect
dealii.demon at gmail.com
dealii.demon at gmail.com
Wed Sep 11 08:05:04 PDT 2013
Revision 1884
add support for periodic domains.
A trunk/aspect/cookbooks/periodic_box.prm
U trunk/aspect/include/aspect/geometry_model/box.h
U trunk/aspect/include/aspect/geometry_model/interface.h
U trunk/aspect/source/geometry_model/box.cc
U trunk/aspect/source/geometry_model/interface.cc
U trunk/aspect/source/simulator/core.cc
http://www.dealii.org/websvn/revision.php?repname=Aspect+Repository&path=%2F&rev=1884&peg=1884
Diff:
Added: trunk/aspect/cookbooks/periodic_box.prm
===================================================================
--- trunk/aspect/cookbooks/periodic_box.prm (rev 0)
+++ trunk/aspect/cookbooks/periodic_box.prm 2013-09-11 15:03:13 UTC (rev 1884)
@@ -0,0 +1,105 @@
+set Dimension = 2
+set CFL number = 1.0
+set End time = 1e9
+set Output directory = output
+set Resume computation = false
+set Start time = 0
+set Adiabatic surface temperature = 0
+set Surface pressure = 0
+set Pressure normalization = no
+set Linear solver tolerance = 1.e-7
+set Use years in output instead of seconds = true
+set Nonlinear solver scheme = IMPES
+
+subsection Boundary temperature model
+ set Model name = box
+ subsection Box
+ set Top temperature = 0.0
+ set Bottom temperature = 1000.0
+ end
+end
+
+
+subsection Discretization
+ set Stokes velocity polynomial degree = 2
+ set Temperature polynomial degree = 2
+ set Use locally conservative discretization = false
+ subsection Stabilization parameters
+ set alpha = 2
+ set beta = 0.078
+ set cR = 0.5 # default: 0.11
+ end
+end
+
+
+subsection Geometry model
+ set Model name = box
+ subsection Box
+ set X periodic = true
+ set X extent = 1.e6
+ set Y extent = 5.e5
+ set Z extent = 5.e5
+ end
+end
+
+
+subsection Gravity model
+ set Model name = vertical
+ subsection Vertical
+ set Magnitude = 10.0
+ end
+end
+
+
+subsection Initial conditions
+ set Model name = function
+ subsection Function
+ set Variable names = x,y
+ set Function expression = if((sqrt((x-1.e5)^2+(y-4.0e5)^2)<5.0e4) | (sqrt((x-3.e5)^2+(y-2.e5)^2)<1.0e5) , 800.0, 0)
+ end
+end
+
+
+subsection Material model
+ set Model name = simple
+ subsection Simple model
+ set Reference density = 3300
+ set Reference specific heat = 1250
+ set Reference temperature = 0.0
+ set Thermal conductivity = 4.7
+ set Thermal expansion coefficient = 4e-5
+ set Viscosity = 1.e20
+ end
+end
+
+
+subsection Mesh refinement
+ set Additional refinement times =
+ set Initial adaptive refinement = 2 # default: 2
+ set Initial global refinement = 6 # default: 2
+ set Refinement fraction = 0.3
+ set Coarsening fraction = 0.03
+ set Strategy = thermal energy density
+ set Time steps between mesh refinement = 10 # default: 10
+end
+
+
+subsection Model settings
+ set Include adiabatic heating = false
+ set Include shear heating = false
+ set Radiogenic heating rate = 0
+ set Fixed temperature boundary indicators = 2,3
+ set Prescribed velocity boundary indicators =
+ set Tangential velocity boundary indicators =
+ set Zero velocity boundary indicators = 2,3
+end
+
+subsection Postprocess
+ set List of postprocessors = visualization
+ subsection Visualization
+ set List of output variables =
+ set Number of grouped files = 1
+ set Output format = vtu
+ set Time between graphical output = 1.e5
+ end
+end
Modified: trunk/aspect/include/aspect/geometry_model/box.h
===================================================================
--- trunk/aspect/include/aspect/geometry_model/box.h 2013-09-11 15:01:51 UTC (rev 1883)
+++ trunk/aspect/include/aspect/geometry_model/box.h 2013-09-11 15:03:13 UTC (rev 1884)
@@ -84,6 +84,13 @@
get_used_boundary_indicators () const;
/**
+ * Return the set of periodic boundaries
+ */
+ virtual
+ std::set< std::pair< std::pair<types::boundary_id, types::boundary_id>, unsigned int> >
+ get_periodic_boundary_pairs () const;
+
+ /**
* Declare the parameters this class takes through input files.
*/
static
@@ -103,6 +110,11 @@
* Extent of the box in x-, y-, and z-direction (in 3d).
*/
Point<dim> extents;
+
+ /**
+ * Flag whether the box is periodic in the x-, y-, and z-direction.
+ */
+ bool periodic[dim];
};
}
}
Modified: trunk/aspect/include/aspect/geometry_model/interface.h
===================================================================
--- trunk/aspect/include/aspect/geometry_model/interface.h 2013-09-11 15:01:51 UTC (rev 1883)
+++ trunk/aspect/include/aspect/geometry_model/interface.h 2013-09-11 15:03:13 UTC (rev 1884)
@@ -114,6 +114,16 @@
get_used_boundary_indicators () const = 0;
/**
+ * Returns a set of periodic boundary pairs. The elements of the set are
+ a pair of boundary ids and a cartesian unit direction each. The base class
+ returns an empty set, so this does nothing unless you specifically use a
+ geometry model with periodic boundary conditions
+ **/
+ virtual
+ std::set< std::pair< std::pair<types::boundary_id, types::boundary_id>, unsigned int> >
+ get_periodic_boundary_pairs () const;
+
+ /**
* Declare the parameters this class takes through input files.
* The default implementation of this function does not describe
* any parameters. Consequently, derived classes do not have to
Modified: trunk/aspect/source/geometry_model/box.cc
===================================================================
--- trunk/aspect/source/geometry_model/box.cc 2013-09-11 15:01:51 UTC (rev 1883)
+++ trunk/aspect/source/geometry_model/box.cc 2013-09-11 15:03:13 UTC (rev 1884)
@@ -25,6 +25,7 @@
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_accessor.h>
+#include <deal.II/grid/grid_tools.h>
namespace aspect
@@ -41,6 +42,21 @@
extents);
for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f)
coarse_grid.begin_active()->face(f)->set_boundary_indicator(f);
+
+ //Tell p4est about the periodicity of the mesh.
+#if (DEAL_II_MAJOR*100 + DEAL_II_MINOR) >= 810
+ std::vector<std_cxx1x::tuple< typename parallel::distributed::Triangulation<dim>::cell_iterator, unsigned int,
+ typename parallel::distributed::Triangulation<dim>::cell_iterator, unsigned int> >
+ periodicity_vector;
+ for( unsigned int i=0; i<dim; ++i)
+ if (periodic[i])
+ GridTools::identify_periodic_face_pairs(coarse_grid, 2*i, 2*i+1, i, periodicity_vector);
+
+ coarse_grid.add_periodicity(periodicity_vector);
+#else
+ for( unsigned int i=0; i<dim; ++i)
+ AssertThrow(!periodic[i],"please update deal.II to the latest version to get support for periodic domains.");
+#endif
}
@@ -56,6 +72,17 @@
return s;
}
+ template <int dim>
+ std::set< std::pair< std::pair<types::boundary_id, types::boundary_id>, unsigned int> >
+ Box<dim>::
+ get_periodic_boundary_pairs () const
+ {
+ std::set< std::pair< std::pair<types::boundary_id, types::boundary_id>, unsigned int> > periodic_boundaries;
+ for( unsigned int i=0; i<dim; ++i)
+ if (periodic[i])
+ periodic_boundaries.insert( std::make_pair( std::pair<types::boundary_id, types::boundary_id>(2*i, 2*i+1), i) );
+ return periodic_boundaries;
+ }
template <int dim>
Point<dim>
@@ -130,6 +157,15 @@
Patterns::Double (0),
"Extent of the box in z-direction. This value is ignored "
"if the simulation is in 2d Units: m.");
+ prm.declare_entry ("X periodic", "false",
+ Patterns::Bool (),
+ "set to true to have the box be periodic in X direction");
+ prm.declare_entry ("Y periodic", "false",
+ Patterns::Bool (),
+ "set to true to have the box be periodic in Y direction");
+ prm.declare_entry ("Z periodic", "false",
+ Patterns::Bool (),
+ "set to true to have the box be periodic in Z direction");
}
prm.leave_subsection();
}
@@ -147,12 +183,19 @@
prm.enter_subsection("Box");
{
extents[0] = prm.get_double ("X extent");
+ periodic[0] = prm.get_bool ("X periodic");
if (dim >= 2)
- extents[1] = prm.get_double ("Y extent");
+ {
+ extents[1] = prm.get_double ("Y extent");
+ periodic[1] = prm.get_bool ("Y periodic");
+ }
if (dim >= 3)
- extents[2] = prm.get_double ("Z extent");
+ {
+ extents[2] = prm.get_double ("Z extent");
+ periodic[2] = prm.get_bool ("Z periodic");
+ }
}
prm.leave_subsection();
}
Modified: trunk/aspect/source/geometry_model/interface.cc
===================================================================
--- trunk/aspect/source/geometry_model/interface.cc 2013-09-11 15:01:51 UTC (rev 1883)
+++ trunk/aspect/source/geometry_model/interface.cc 2013-09-11 15:03:13 UTC (rev 1884)
@@ -34,8 +34,15 @@
{}
+ template<int dim>
+ std::set< std::pair< std::pair<types::boundary_id, types::boundary_id>, unsigned int > >
+ Interface<dim>::get_periodic_boundary_pairs() const
+ {
+ //return an empty set in the base class
+ return std::set< std::pair< std::pair< types::boundary_id, types::boundary_id>, unsigned int > >();
+ }
+
template <int dim>
-
void
Interface<dim>::
declare_parameters (dealii::ParameterHandler &prm)
Modified: trunk/aspect/source/simulator/core.cc
===================================================================
--- trunk/aspect/source/simulator/core.cc 2013-09-11 15:01:51 UTC (rev 1883)
+++ trunk/aspect/source/simulator/core.cc 2013-09-11 15:03:13 UTC (rev 1884)
@@ -613,17 +613,51 @@
pcout.get_stream().imbue(s);
}
+ //reinit the constraints matrix and make hanging node constraints
+ constraints.clear();
+ constraints.reinit(introspection.index_sets.system_relevant_set);
+ DoFTools::make_hanging_node_constraints (dof_handler,
+ constraints);
+
+ //Now set up the constraints for periodic boundary conditions
+ {
+ typedef std::set< std::pair< std::pair< types::boundary_id, types::boundary_id>, unsigned int> >
+ periodic_boundary_set;
+ periodic_boundary_set pbs = geometry_model->get_periodic_boundary_pairs();
+
+ for(periodic_boundary_set::iterator p = pbs.begin(); p != pbs.end(); ++p)
+ {
+ //Throw error if we are trying to use the same boundary for more than one boundary condition
+ Assert( is_element( (*p).first.first, parameters.fixed_temperature_boundary_indicators ) == false &&
+ is_element( (*p).first.second, parameters.fixed_temperature_boundary_indicators ) == false &&
+ is_element( (*p).first.first, parameters.zero_velocity_boundary_indicators ) == false &&
+ is_element( (*p).first.second, parameters.zero_velocity_boundary_indicators ) == false &&
+ is_element( (*p).first.first, parameters.tangential_velocity_boundary_indicators ) == false &&
+ is_element( (*p).first.second, parameters.tangential_velocity_boundary_indicators ) == false &&
+ parameters.prescribed_velocity_boundary_indicators.find( (*p).first.first)
+ == parameters.prescribed_velocity_boundary_indicators.end() &&
+ parameters.prescribed_velocity_boundary_indicators.find( (*p).first.second)
+ == parameters.prescribed_velocity_boundary_indicators.end(),
+ ExcInternalError());
+
+#if (DEAL_II_MAJOR*100 + DEAL_II_MINOR) >= 810
+ DoFTools::make_periodicity_constraints(dof_handler,
+ (*p).first.first, //first boundary id
+ (*p).first.second, //second boundary id
+ (*p).second, //cartesian direction for translational symmetry
+ constraints);
+#endif
+ }
+
+
+ }
// then compute constraints for the velocity. the constraints we compute
// here are the ones that are the same for all following time steps. in
// addition, we may be computing constraints from boundary values for the
// velocity that are different between time steps. these are then put
// into current_constraints in start_timestep().
{
- constraints.clear();
- constraints.reinit(introspection.index_sets.system_relevant_set);
- DoFTools::make_hanging_node_constraints (dof_handler,
- constraints);
// do the interpolation for zero velocity
for (std::set<types::boundary_id>::const_iterator
@@ -672,8 +706,8 @@
// we do nothing with the compositional fields: homogeneous Neumann boundary conditions
- constraints.close();
}
+ constraints.close();
// finally initialize vectors, matrices, etc.
More information about the CIG-COMMITS
mailing list