Skip to content

Conversation

buchanankerswell
Copy link
Contributor

This PR adds a coordinate system option to source/particle/property/function.cc. It follows from the starter pack issue #5626.

I attempted to mimic @gassmoeller's examples in #5626, but

ctest -R particle_property_multiple_functions_spherical_coordinates -V

is failing. Apparently the parameter set Coordinate system = spherical is not declared correctly in, but I cannot figure out why. Here is the relevant part of source/particle/property/function.cc:

      template <int dim>
      void
      Function<dim>::declare_parameters (ParameterHandler &prm)
      {
        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 "
                             "by a function expression delimited by a ';'.");
          Functions::ParsedFunction<dim>::declare_parameters (prm, 1);
        }
        prm.leave_subsection();
      }


      template <int dim>
      void
      Function<dim>::parse_parameters (ParameterHandler &prm)
      {
        prm.enter_subsection("Function");
        {
          n_components = prm.get_integer ("Number of components");
          try
          {
            function = std::make_unique<Functions::ParsedFunction<dim>>(n_components);
            function->parse_parameters (prm);
          }
          catch (...)
          {
            std::cerr << "ERROR: FunctionParser failed to parse\n"
              << "\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();
      }

Here is the relevant part of the prm file tests/particle_property_multiple_functions_spherical_coordinates.prm:

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

Here is part of the ctest output on my local machine:

790: ----------------------------------------------------
790: Exception 'ExcCannotParseLine(current_line_n, input_filename, ("No entry with name <" + entry_name + "> was declared in the current subsection."))' on rank 0 on processing:
790:
790: --------------------------------------------------------
790: An error occurred in line <2035> of file </Users/bck/Working/kerswell_et_al_dynp/dealii/tmp/unpack/deal.II-v9.6.1/source/base/parameter_handler.cc> in function
790:     void dealii::ParameterHandler::scan_line(std::string, const std::string &, const unsigned int, const bool)
790: The violated condition was:
790:     skip_undefined
790: Additional information:
790:     Line <105> of file <input string>: No entry with name <Coordinate
790:     system> was declared in the current subsection.
790:
790: Stacktrace:
790: -----------
790: #0  2   libdeal_II.g.9.6.1.dylib            0x0000000114177914 _ZN6dealii16ParameterHandler9scan_lineENSt3__112basic_stringIcNS1_11char_traitsIcEENS1_9allocatorIcEEEERKS7_jb + 6244: 2   libdeal_II.g.9.6.1.dylib            0x0000000114177914 _ZN6dealii16ParameterHandler9scan_lineENSt3__112basic_stringIcNS1_11char_traitsIcEENS1_9allocatorIcEEEERKS7_jb
790: #1  3   libdeal_II.g.9.6.1.dylib            0x0000000114168794 _ZZN6dealii16ParameterHandler11parse_inputERNSt3__113basic_istreamIcNS1_11char_traitsIcEEEERKNS1_12basic_stringIcS4_NS1_9allocatorIcEEEESC_bENK3$_0clESC_SC_j + 120: 3   libdeal_II.g.9.6.1.dylib            0x0000000114168794 _ZZN6dealii16ParameterHandler11parse_inputERNSt3__113basic_istreamIcNS1_11char_traitsIcEEEERKNS1_12basic_stringIcS4_NS1_9allocatorIcEEEESC_bENK3$_0clESC_SC_j
790: #2  4   libdeal_II.g.9.6.1.dylib            0x00000001141681fc _ZN6dealii16ParameterHandler11parse_inputERNSt3__113basic_istreamIcNS1_11char_traitsIcEEEERKNS1_12basic_stringIcS4_NS1_9allocatorIcEEEESC_b + 676: 4   libdeal_II.g.9.6.1.dylib            0x00000001141681fc _ZN6dealii16ParameterHandler11parse_inputERNSt3__113basic_istreamIcNS1_11char_traitsIcEEEERKNS1_12basic_stringIcS4_NS1_9allocatorIcEEEESC_b
790: #3  5   libdeal_II.g.9.6.1.dylib            0x0000000114169730 _ZN6dealii16ParameterHandler23parse_input_from_stringERKNSt3__112basic_stringIcNS1_11char_traitsIcEENS1_9allocatorIcEEEES9_b + 116: 5   libdeal_II.g.9.6.1.dylib            0x0000000114169730 _ZN6dealii16ParameterHandler23parse_input_from_stringERKNSt3__112basic_stringIcNS1_11char_traitsIcEENS1_9allocatorIcEEEES9_b
790: #4  6   aspect-debug                        0x000000010124f2ec _Z16parse_parametersRKNSt3__112basic_stringIcNS_11char_traitsIcEENS_9allocatorIcEEEERN6dealii16ParameterHandlerE + 256: 6   aspect-debug                        0x000000010124f2ec _Z16parse_parametersRKNSt3__112basic_stringIcNS_11char_traitsIcEENS_9allocatorIcEEEERN6dealii16ParameterHandlerE
790: #5  7   aspect-debug                        0x0000000101250500 _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbbb + 268: 7   aspect-debug                        0x0000000101250500 _Z13run_simulatorILi2EEvRKNSt3__112basic_stringIcNS0_11char_traitsIcEENS0_9allocatorIcEEEES8_bbbb
790: #6  8   aspect-debug                        0x000000010124fe80 main + 2180: 8   aspect-debug                        0x000000010124fe80 main
790: #7  9   dyld                                0x000000019435ab98 start + 6076: 9   dyld                                0x000000019435ab98 start
790: --------------------------------------------------------
790:
790: Aborting!
790: ----------------------------------------------------

Lastly, I copy/pasted the test files from tests/particle_property_multiple_functions, so you can ignore those until the test actually runs.

Copy link
Member

@gassmoeller gassmoeller left a comment

Choose a reason for hiding this comment

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

Thanks for tackling this. There are a few issues with the PR as you already noticed. I left comments for how to fix them. Afterwards the test will run without issues (I tested that locally and it works).

After fixing my comments please:

  • add a changelog entry
  • update the test results
  • make sure you run make indent

Then this should make a nice addition.


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

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);
}

Copy link
Member

Choose a reason for hiding this comment

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

delete this file, it is one of the Paraview files we do not usually include in test output. the same for all the other paraview output files, replace them with the gnuplot output files.


subsection Particles
set Time between data output = 70
set Data output format = vtu
Copy link
Member

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).

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

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.

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

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

@buchanankerswell buchanankerswell force-pushed the add-coordinate-system-particle-property-function branch from 306c097 to d9307ab Compare July 22, 2025 11:52
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

2 participants