Skip to content

[WIP] Add coordinate system to particle/property/function #6544

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
wants to merge 1 commit into
base: main
Choose a base branch
from
Open

[WIP] Add coordinate system to particle/property/function #6544

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
6 changes: 6 additions & 0 deletions include/aspect/particle/property/function.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -88,6 +88,12 @@ namespace aspect
*/
std::unique_ptr<Functions::ParsedFunction<dim>> function;

/**
* The coordinate representation to evaluate the function. Possible
* choices are depth, cartesian and spherical.
*/
Utilities::Coordinates::CoordinateSystem coordinate_system;
Copy link
Member

@gassmoeller gassmoeller Jun 24, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In order for this to work you need to add #include <aspect/utilities.h> at the top of the file. Otherwise CoordinateSystem is unknown in this file.


/**
* A private variable that stores the number of particle property
* function components.
Expand Down
32 changes: 25 additions & 7 deletions source/particle/property/function.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
*/

#include <aspect/particle/property/function.h>
#include <aspect/geometry_model/interface.h>

namespace aspect
{
Expand All @@ -37,8 +38,11 @@ namespace aspect
Function<dim>::initialize_one_particle_property(const Point<dim> &position,
std::vector<double> &data) const
{
// convert the position into the selected coordinate system
const Utilities::NaturalCoordinate<dim> point = this->get_geometry_model().cartesian_to_other_coordinates(position, coordinate_system);
Copy link
Member

@gassmoeller gassmoeller Jun 24, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

In order for this to work you also need to derive the Function class from SimulatorAccess(and include#include <aspect/simulator_access.h>in the header file). SimulatorAccess is the class that provides theget_geometry_model` function. You can see how that works for the initial temperature function object here: https://github.com/geodynamics/aspect/blob/main/include/aspect/initial_temperature/function.h#L42


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<dim>(point.get_coordinates()), i));
Copy link
Member

@gassmoeller gassmoeller Jun 24, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

this is correct, but a bit hard to read. what about:

Suggested change
data.push_back(function->value(Utilities::convert_array_to_point<dim>(point.get_coordinates()), i));
{
const double function_value = function->value(Utilities::convert_array_to_point<dim>(point.get_coordinates()), i)
data.push_back(function_value);
}

}

template <int dim>
Expand All @@ -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 "
Expand All @@ -71,20 +85,24 @@ namespace aspect
Function<dim>::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<Functions::ParsedFunction<dim>>(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();
}
}
Expand Down
118 changes: 118 additions & 0 deletions tests/particle_property_multiple_functions_spherical_coordinates.prm
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
Copy link
Member

@gassmoeller gassmoeller Jun 24, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

can you set this to gnuplot and include the created gnuplot data files instead of the Paraview files? The benefit is that gnuplot files are text files and it is therefore easier to see how the output changes (compared to the binary files of Paraview).

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
10 changes: 10 additions & 0 deletions tests/particle_property_multiple_functions_spherical_coordinates/particles.pvd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,10 @@
<?xml version="1.0"?>
<!--
#This file was generated by the deal.II library on 2016/9/20 at 14:39:51
-->
<VTKFile type="Collection" version="0.1" ByteOrder="LittleEndian">
<Collection>
<DataSet timestep="0" group="" part="0" file="particles/particles-00000.pvtu"/>
<DataSet timestep="70" group="" part="0" file="particles/particles-00001.pvtu"/>
</Collection>
</VTKFile>
47 changes: 47 additions & 0 deletions ...icle_property_multiple_functions_spherical_coordinates/particles/particles-00000.0000.vtu
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
<?xml version="1.0" ?>
<!--
# vtk DataFile Version 3.0
#This file was generated
<VTKFile type="UnstructuredGrid" version="0.1" compressor="vtkZLibDataCompressor" byte_order="LittleEndian">
<UnstructuredGrid>
<FieldData>
<DataArray type="Float32" Name="TIME" NumberOfTuples="1" format="ascii">0</DataArray>
</FieldData>
<Piece NumberOfPoints="3" NumberOfCells="3" >
<Points>
<DataArray type="Float32" NumberOfComponents="3" format="binary">
AQAAACQAAAAkAAAAJQAAAA==eAH7d0PWXrX8gh0DENzxSref7LUUzJ630cB+SfwzMBsA6zYLTg==
</DataArray>
</Points>

<Cells>
<DataArray type="Int32" Name="connectivity" format="binary">
AQAAAAwAAAAMAAAAEQAAAA==eAFjYGBgYARiJiAGAAAcAAQ=
</DataArray>
<DataArray type="Int32" Name="offsets" format="binary">
AQAAAAwAAAAMAAAAEQAAAA==eAFjZGBgYAJiZiAGAAA0AAc=
</DataArray>
<DataArray type="UInt8" Name="types" format="binary">
AQAAAAMAAAADAAAACwAAAA==eAFjZGQEAAAJAAQ=
</DataArray>
</Cells>
<PointData Scalars="scalars">
<DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="binary">
AQAAACQAAAAkAAAACwAAAA==eAFjYCAMAAAkAAE=
</DataArray>
<DataArray type="Float32" Name="function" NumberOfComponents="3" format="binary">
AQAAACQAAAAkAAAAHwAAAA==eAFjYGBgmBN2WB9IgUFZhIYpjM2wMFMHxAYAWAkERg==
</DataArray>
<DataArray type="Float32" Name="initial position" NumberOfComponents="3" format="binary">
AQAAACQAAAAkAAAAJQAAAA==eAH7d0PWXrX8gh0DENzxSref7LUUzJ630cB+SfwzMBsA6zYLTg==
</DataArray>
<DataArray type="Float32" Name="id" format="binary">
AQAAAAwAAAAMAAAAEAAAAA==eAFjYACBBnsg4QAABIcBAA==
</DataArray>
<DataArray type="Float32" Name="initial C_1" format="binary">
AQAAAAwAAAAMAAAAFAAAAA==eAGbE3ZYvyxCw5RhYaYOAB9FBEY=
</DataArray>
</PointData>
</Piece>
</UnstructuredGrid>
</VTKFile>
47 changes: 47 additions & 0 deletions ...icle_property_multiple_functions_spherical_coordinates/particles/particles-00000.0001.vtu
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,47 @@
<?xml version="1.0" ?>
<!--
# vtk DataFile Version 3.0
#This file was generated
<VTKFile type="UnstructuredGrid" version="0.1" compressor="vtkZLibDataCompressor" byte_order="LittleEndian">
<UnstructuredGrid>
<FieldData>
<DataArray type="Float32" Name="TIME" NumberOfTuples="1" format="ascii">0</DataArray>
</FieldData>
<Piece NumberOfPoints="7" NumberOfCells="7" >
<Points>
<DataArray type="Float32" NumberOfComponents="3" format="binary">
AQAAAFQAAABUAAAASAAAAA==eAFjn6Rm77LF2p4BCFiffbNbrl0EZnflcNjbs5WA2cXxCfYr2z3B7MioSHvRc6Fg9rKgAHvuUylg9qK+KPt3E+rAbACDXRU2
</DataArray>
</Points>

<Cells>
<DataArray type="Int32" Name="connectivity" format="binary">
AQAAABwAAAAcAAAAGwAAAA==eAFjYGBgYARiJiBmBmIWIGYFYjYgBgAA/AAW
</DataArray>
<DataArray type="Int32" Name="offsets" format="binary">
AQAAABwAAAAcAAAAGwAAAA==eAFjZGBgYAJiZiBmAWJWIGYDYnYgBgABbAAd
</DataArray>
<DataArray type="UInt8" Name="types" format="binary">
AQAAAAcAAAAHAAAACwAAAA==eAFjZAQDAAAjAAg=
</DataArray>
</Cells>
<PointData Scalars="scalars">
<DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="binary">
AQAAAFQAAABUAAAADAAAAA==eAFjYKA+AAAAVAAB
</DataArray>
<DataArray type="Float32" Name="function" NumberOfComponents="3" format="binary">
AQAAAFQAAABUAAAADAAAAA==eAFjYKA+AAAAVAAB
</DataArray>
<DataArray type="Float32" Name="initial position" NumberOfComponents="3" format="binary">
AQAAAFQAAABUAAAASAAAAA==eAFjn6Rm77LF2p4BCFiffbNbrl0EZnflcNjbs5WA2cXxCfYr2z3B7MioSHvRc6Fg9rKgAHvuUylg9qK+KPt3E+rAbACDXRU2
</DataArray>
<DataArray type="Float32" Name="id" format="binary">
AQAAABwAAAAcAAAAHgAAAA==eAFjYHBwYGBoAOIFQHwAiB8AMYMjA4OAIwBHAgTT
</DataArray>
<DataArray type="Float32" Name="initial C_1" format="binary">
AQAAABwAAAAcAAAACwAAAA==eAFjYMANAAAcAAE=
</DataArray>
</PointData>
</Piece>
</UnstructuredGrid>
</VTKFile>
19 changes: 19 additions & 0 deletions ...particle_property_multiple_functions_spherical_coordinates/particles/particles-00000.pvtu
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,19 @@
<?xml version="1.0"?>
<!--
#This file was generated
<VTKFile type="PUnstructuredGrid" version="0.1" byte_order="LittleEndian">
<PUnstructuredGrid GhostLevel="0">
<PPointData Scalars="scalars">
<PDataArray type="Float32" Name="velocity" NumberOfComponents="3" format="ascii"/>
<PDataArray type="Float32" Name="function" NumberOfComponents="3" format="ascii"/>
<PDataArray type="Float32" Name="initial position" NumberOfComponents="3" format="ascii"/>
<PDataArray type="Float32" Name="id" format="ascii"/>
<PDataArray type="Float32" Name="initial C_1" format="ascii"/>
</PPointData>
<PPoints>
<PDataArray type="Float32" NumberOfComponents="3"/>
</PPoints>
<Piece Source="particles-00000.0000.vtu"/>
<Piece Source="particles-00000.0001.vtu"/>
</PUnstructuredGrid>
</VTKFile>
48 changes: 48 additions & 0 deletions ...icle_property_multiple_functions_spherical_coordinates/particles/particles-00001.0000.vtu
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,48 @@
<?xml version="1.0" ?>
<!--
# vtk DataFile Version 3.0
#This file was generated
<VTKFile type="UnstructuredGrid" version="0.1" compressor="vtkZLibDataCompressor" byte_order="LittleEndian">
<UnstructuredGrid>
<FieldData>
<DataArray type="Float32" Name="CYCLE" NumberOfTuples="1" format="ascii">1</DataArray>
<DataArray type="Float32" Name="TIME" NumberOfTuples="1" format="ascii">70</DataArray>
</FieldData>
<Piece NumberOfPoints="3" NumberOfCells="3" >
<Points>
<DataArray type="Float32" NumberOfComponents="3" format="binary">
AQAAACQAAAAkAAAAJQAAAA==eAHTuS5nP0v1tB0DEDxyTrcva5wDZk9/bGh/5Ot9MBsA2bMLhQ==
</DataArray>
</Points>

<Cells>
<DataArray type="Int32" Name="connectivity" format="binary">
AQAAAAwAAAAMAAAAEQAAAA==eAFjYGBgYARiJiAGAAAcAAQ=
</DataArray>
<DataArray type="Int32" Name="offsets" format="binary">
AQAAAAwAAAAMAAAAEQAAAA==eAFjZGBgYAJiZiAGAAA0AAc=
</DataArray>
<DataArray type="UInt8" Name="types" format="binary">
AQAAAAMAAAADAAAACwAAAA==eAFjZGQEAAAJAAQ=
</DataArray>
</Cells>
<PointData Scalars="scalars">
<DataArray type="Float32" Name="velocity" NumberOfComponents="3" format="binary">
AQAAACQAAAAkAAAAJQAAAA==eAGbbB1qMWWB7E4GIPDKZ91WMrsezF4m1GKxKN0OzAYA1EsKlQ==
</DataArray>
<DataArray type="Float32" Name="function" NumberOfComponents="3" format="binary">
AQAAACQAAAAkAAAAHwAAAA==eAFjYGBgmBN2WB9IgUFZhIYpjM2wMFMHxAYAWAkERg==
</DataArray>
<DataArray type="Float32" Name="initial position" NumberOfComponents="3" format="binary">
AQAAACQAAAAkAAAAJQAAAA==eAH7d0PWXrX8gh0DENzxSref7LUUzJ630cB+SfwzMBsA6zYLTg==
</DataArray>
<DataArray type="Float32" Name="id" format="binary">
AQAAAAwAAAAMAAAAEAAAAA==eAFjYACBBnsg4QAABIcBAA==
</DataArray>
<DataArray type="Float32" Name="initial C_1" format="binary">
AQAAAAwAAAAMAAAAFAAAAA==eAGbE3ZYvyxCw5RhYaYOAB9FBEY=
</DataArray>
</PointData>
</Piece>
</UnstructuredGrid>
</VTKFile>
Loading
Loading