diff --git a/include/aspect/particle/property/function.h b/include/aspect/particle/property/function.h index f45d52de16c..84b7fc4ff09 100644 --- a/include/aspect/particle/property/function.h +++ b/include/aspect/particle/property/function.h @@ -88,6 +88,12 @@ namespace aspect */ std::unique_ptr> function; + /** + * The coordinate representation to evaluate the function. Possible + * choices are depth, cartesian and spherical. + */ + Utilities::Coordinates::CoordinateSystem coordinate_system; + /** * A private variable that stores the number of particle property * function components. diff --git a/source/particle/property/function.cc b/source/particle/property/function.cc index eb5747d9e0d..78bd437c014 100644 --- a/source/particle/property/function.cc +++ b/source/particle/property/function.cc @@ -19,6 +19,7 @@ */ #include +#include namespace aspect { @@ -37,8 +38,11 @@ namespace aspect Function::initialize_one_particle_property(const Point &position, std::vector &data) const { + // convert the position into the selected coordinate system + const Utilities::NaturalCoordinate point = this->get_geometry_model().cartesian_to_other_coordinates(position, coordinate_system); + for (unsigned int i = 0; i < n_components; ++i) - data.push_back(function->value(position, i)); + data.push_back(function->value(Utilities::convert_array_to_point(point.get_coordinates()), i)); } template @@ -56,6 +60,16 @@ namespace aspect { prm.enter_subsection("Function"); { + prm.declare_entry ("Coordinate system", "cartesian", + Patterns::Selection ("cartesian|spherical|depth"), + "A selection that determines the assumed coordinate " + "system for the function variables. Allowed values " + "are `cartesian', `spherical', and `depth'. `spherical' coordinates " + "are interpreted as r,phi or r,phi,theta in 2d/3d " + "respectively with theta being the polar angle. `depth' " + "will create a function, in which only the first " + "parameter is non-zero, which is interpreted to " + "be the depth of the point."); prm.declare_entry ("Number of components", "1", Patterns::Integer (0), "The number of function components where each component is described " @@ -71,20 +85,24 @@ namespace aspect Function::parse_parameters (ParameterHandler &prm) { prm.enter_subsection("Function"); - n_components = prm.get_integer ("Number of components"); - try + { + n_components = prm.get_integer ("Number of components"); + try { function = std::make_unique>(n_components); function->parse_parameters (prm); } - catch (...) + catch (...) { std::cerr << "ERROR: FunctionParser failed to parse\n" - << "\t'Particles.Function'\n" - << "with expression\n" - << "\t'" << prm.get("Function expression") << "'"; + << "\t'Particles.Function'\n" + << "with expression\n" + << "\t'" << prm.get("Function expression") << "'"; throw; } + + coordinate_system = Utilities::Coordinates::string_to_coordinate_system(prm.get("Coordinate system")); + } prm.leave_subsection(); } } diff --git a/tests/particle_property_multiple_functions_spherical_coordinates.prm b/tests/particle_property_multiple_functions_spherical_coordinates.prm new file mode 100644 index 00000000000..c9523e1aae9 --- /dev/null +++ b/tests/particle_property_multiple_functions_spherical_coordinates.prm @@ -0,0 +1,118 @@ +# This test is identical to particle_property_multiple_functions.prm +# but with spherical coordinates. +# +# A description of the van Keken et al. benchmark using smooth initial +# conditions instead of the discontinuous ones in van-keken-discontinuous.prm. +# See the manual for more information. This test adds particles to the cookbook +# parameter file with two particle properties defined by two different functions. +# One is a discontinuous function and the other is a smooth function. + +# MPI: 2 + +set Dimension = 2 +set End time = 70 +set Use years in output instead of seconds = false + +subsection Geometry model + set Model name = box + + subsection Box + set X extent = 0.9142 + set Y extent = 1.0000 + end +end + +subsection Boundary velocity model + set Tangential velocity boundary indicators = left, right + set Zero velocity boundary indicators = bottom, top +end + +subsection Material model + set Model name = simple + + subsection Simple model + set Reference density = 1010 + set Viscosity = 1e2 + set Thermal expansion coefficient = 0 + set Density differential for compositional field 1 = -10 + end +end + +subsection Gravity model + set Model name = vertical + + subsection Vertical + set Magnitude = 10 + end +end + +############### Parameters describing the temperature field +# Note: The temperature plays no role in this model + +subsection Initial temperature model + set Model name = function + + subsection Function + set Function expression = 0 + end +end + +############### Parameters describing the compositional field +# Note: The compositional field is what drives the flow +# in this example + +subsection Compositional fields + set Number of fields = 1 +end + +subsection Initial composition model + set Model name = function + + subsection Function + set Variable names = x,z + set Function constants = pi=3.1415926 + set Function expression = 0.5*(1+tanh((0.2+0.02*cos(pi*x/0.9142)-z)/0.02)) + end +end + +############### Parameters describing the discretization + +subsection Mesh refinement + set Initial adaptive refinement = 0 + set Strategy = composition + set Initial global refinement = 4 + set Time steps between mesh refinement = 0 + set Coarsening fraction = 0.05 + set Refinement fraction = 0.3 +end + +############### Parameters describing what to do with the solution + +subsection Postprocess + set List of postprocessors = velocity statistics, composition statistics, particles + + subsection Particles + set Time between data output = 70 + set Data output format = vtu + end +end + +subsection Particles + set List of particle properties = velocity, function, initial composition, initial position + set Particle generator name = probability density function + + subsection Function + set Coordinate system = spherical + set Number of components = 2 + set Variable names = x,z + set Function expression = if( (z>0.2+0.02*cos(pi*x/0.9142)) , 0 , 1 ); 0.5*(1+tanh((0.2+0.02*cos(pi*x/0.9142)-z)/0.02)) + end + + subsection Generator + subsection Probability density function + set Variable names = x,z + set Function expression = x*x*z + set Number of particles = 10 + end + end +end diff --git a/tests/particle_property_multiple_functions_spherical_coordinates/particles.pvd b/tests/particle_property_multiple_functions_spherical_coordinates/particles.pvd new file mode 100644 index 00000000000..4df1e418fb4 --- /dev/null +++ b/tests/particle_property_multiple_functions_spherical_coordinates/particles.pvd @@ -0,0 +1,10 @@ + + + + + + + + diff --git a/tests/particle_property_multiple_functions_spherical_coordinates/particles/particles-00000.0000.vtu b/tests/particle_property_multiple_functions_spherical_coordinates/particles/particles-00000.0000.vtu new file mode 100644 index 00000000000..265632a84e1 --- /dev/null +++ b/tests/particle_property_multiple_functions_spherical_coordinates/particles/particles-00000.0000.vtu @@ -0,0 +1,47 @@ + +