Grain size model with density compositional fields

Hi,

I have been running some Rayleigh instability analysis models with two different density fluids (ice). I am trying to incorporate the grain size material model but can’t prescribe the two different densities for the compositional fields with that model. I see the material files, but the composition shouldn’t depend on pressure/temperature (and temperature dependencies are ignored) since I am not basing the fields on depth-dependent phase transitions. Is there a way around this?

Thanks!
Max

Hi Max,

Could you be a bit more specific about what the issue is and post the relevant sections of your input file? I think you should be able to read in different material files for different compositional fields (Have you made sure the files are actually being used, for example, that “Use table properties” is set to true?), and these material files can specify different densities.

Or is your question about how you have to change the material files so that the density no longer depends on temperature/pressure? If you want your density difference to be constant, can’t you simply make a first table with the properties you want for ice, and then create a second table that simply is the same as the first but with a constant value added?

Best,
Juliane

Hi Juliane,

Thanks for your reply. Yes, your second point is what I am trying to get at. I’m just not sure how to apply look up tables for the different compositions that are independent of temperature and pressure. Would the tables be the material files?

I would have something like

subsection Compositional fields
  set Number of fields = 3
  set Names of fields   = pure, saline, grain_size
end

subsection Initial composition model
  set Model name =   function

  subsection Function
    set Function constants  = A=50, W=500, C=40e3
    set Function expression = if(y>29e3+A*cos(pi*(x-C)/W),0,1) ; if(y>29e3+A*cos(pi*(x-C)/W),1,0); 1e-3
  end
end

where the first composition layer is the saline ice and the underlying region is the pure ice.

Then I would have something like:

subsection Material model
  set Model name = grain size


  subsection Grain size model
    set Data directory                   = ...
    set Material file names              = saline_ice.txt, pure_ice.txt
    ...
    end
end

Although, the only material file I can find is “pyr-ringwood88.txt” which includes a range of temperatures and pressures and I’m not sure how the material files are applied to the compositions, if it’s tied to a naming convention or just the order of the material files given.

Hi Max,

Yes, the look-up tables are the material files and we provide the “pyr-ringwood88.txt” file as an example. Looking at the “pyr-ringwood88.txt” shows you the format of the material table, and the format is also documented in this cookbook. Although if you want all your properties to be independent of temperature and pressure, a better starting point might be this file since it only has four lines (for all combinations of minimum temperature, max temperature, min pressure, max pressure). In your case, simply set all density, alpha, cp values to be the same in each line—making them independent of pressure and temperature. (But note that this file uses Hefesto and not Perplex format, which is also something you need to specify in the input file.)

As explained in the documentation of the “Material file names” parameter, it is a list with as many components as active compositional fields and material data is assumed to be in order with the ordering of the fields. So in your example, you would still need to swap the file names since your compositional fields list the pure ice first.

Best,
Juliane

That all makes sense, thanks for the advice!
Max

I keep recieving this error:

Exception 'ExcMessage("parse_map_to_double_array can only check the structure " "of the parsed map for " + options.property_name + " if an expected number of values for each key is given.")' on rank 0 on processing: 

--------------------------------------------------------
An error occurred in line <526> of file </Users/maxcollins/Documents/Research/Ice Dynamics/aspect/source/utilities.cc> in function
    std::vector<double> aspect::Utilities::MapParsing::parse_map_to_double_array(const std::string &, Options &)
The violated condition was: 
    options.check_values_per_key == false || options.n_values_per_key.size() == options.list_of_required_keys.size()
Additional information: 
    parse_map_to_double_array can only check the structure of the parsed
    map for Angles of internal friction if an expected number of values
    for each key is given.

for both hefesto or perplex formatted material files. I haven’t touched the angles of internal friction. I have also tried just one compositional field (other than grain size) and still recieve the same issues.

Hi Max,

This is a problem that I think we fixed last month. So if you update your ASPECT to a more recent version, it should work. (Otherwise, if you want to keep the version you have, I think you can specify as many “Angles of internal friction” and “Cohesions” as you have chemical fields (i.e. specifying composition rather than grain size), which is what the error message complains about).

Best,
Juliane

Ahh okay I see, thank you. If I were to put that into the input file what would that look like?

As a list or keys? I work on an HPC so it will take awhile to get updated.

    set Angles of internal friction          = pure:0.0, saline:0.0
    set Cohesions                            = pure:1e20, saline:1e20

or

    set Angles of internal friction          = 0.0, 0.0
    set Cohesions                            = 1e20, 1e20

or

    set Angles of internal friction          = pure:0.0 | saline:0.0
    set Cohesions                            = pure:1e20 | saline:1e20

These all still give me the error.

Quick update, installed Aspect locally with the most recently pulled version from github. I’ve gotten through the other errors but just continuously get a segmentation fault when I use my provided data files. Is there something I’ve entered wrong or something is getting divided and returning a NaN that I’m missing?

pure_ice.txt (1005 Bytes)

saline_ice_940.txt (1005 Bytes)

input_940_grainsize_v2.prm (5.2 KB)

@maxqcollins I assume you’re running in debug mode when you get the segmentation fault. Can you run the program in a debugger and get a backtrace that shows where the segmentation fault happens?

Best

W.

@bangerth Thanks for the reply! The output I get is:

❯ ./aspect-debug ../runs/input_940_grainsize_v2.prm
-----------------------------------------------------------------------------
--                             This is ASPECT                              --
-- The Advanced Solver for Planetary Evolution, Convection, and Tectonics. --
-----------------------------------------------------------------------------
--     . version 3.1.0-pre (main, 17e66766c)
--     . using deal.II 9.6.2
--     .       with 32 bit indices
--     .       with vectorization level 1 (SSE2, 2 doubles, 128 bits)
--     . using Trilinos 13.4.1
--     . using p4est 2.8.0
--     . using Geodynamic World Builder 1.0.0
--     . running in DEBUG mode
--     . running with 1 MPI process
-----------------------------------------------------------------------------

[1]    53211 segmentation fault  ./aspect-debug ../runs/input_940_grainsize_v2.prm

It seems like running it on the computing cluster gives a different output with the same input file:

-----------------------------------------------------------------------------
--                             This is ASPECT                              --
-- The Advanced Solver for Planetary Evolution, Convection, and Tectonics. --
-----------------------------------------------------------------------------
--     . version 3.1.0-pre (main, 0a1784a42)
--     . using deal.II 9.5.1
--     .       with 32 bit indices
--     .       with vectorization level 0 (disabled, 1 doubles, 64 bits)
--     . using Trilinos 13.2.0
--     . using p4est 2.3.6
--     . using Geodynamic World Builder 1.0.0
--     . running in DEBUG mode
--     . running with 1 MPI process
-----------------------------------------------------------------------------

/usr/include/c++/11/bits/stl_vector.h:1045: std::vector<_Tp, _Alloc>::reference std::vector<_Tp, _Alloc>::operator[](std::vector<_Tp, _Alloc>::size_type) [with _Tp = std::__cxx11::basic_string<char>; _Alloc = std::allocator<std::__cxx11::basic_string<char> >; std::vector<_Tp, _Alloc>::reference = std::__cxx11::basic_string<char>&; std::vector<_Tp, _Alloc>::size_type = long unsigned int]: Assertion '__n < this->size()' failed.
[quser42:2814811] *** Process received signal ***
[quser42:2814811] Signal: Aborted (6)
[quser42:2814811] Signal code:  (-6)
[quser42:2814811] [ 0] /lib/x86_64-linux-gnu/libc.so.6(+0x42520)[0x7f0a4ed15520]
[quser42:2814811] [ 1] /lib/x86_64-linux-gnu/libc.so.6(pthread_kill+0x12c)[0x7f0a4ed699fc]
[quser42:2814811] [ 2] /lib/x86_64-linux-gnu/libc.so.6(raise+0x16)[0x7f0a4ed15476]
[quser42:2814811] [ 3] /lib/x86_64-linux-gnu/libc.so.6(abort+0xd3)[0x7f0a4ecfb7f3]
[quser42:2814811] [ 4] /opt/aspect/bin/aspect(_ZSt17__size_to_integerm+0x0)[0x55b5ce455053]
[quser42:2814811] [ 5] /opt/aspect/bin/aspect(_ZNSt6vectorINSt7__cxx1112basic_stringIcSt11char_traitsIcESaIcEEESaIS5_EEixEm+0x63)[0x55b5ce553aed]
[quser42:2814811] [ 6] /opt/aspect/bin/aspect(_ZN6aspect13MaterialModel9GrainSizeILi2EE10initializeEv+0x22f)[0x55b5cfec1295]
[quser42:2814811] [ 7] /opt/aspect/bin/aspect(_ZN6aspect9SimulatorILi2EEC2EP19ompi_communicator_tRN6dealii16ParameterHandlerE+0xe55)[0x55b5cf321e29]
[quser42:2814811] [ 8] /opt/aspect/bin/aspect(+0x48c93da)[0x55b5d03e33da]
[quser42:2814811] [ 9] /opt/aspect/bin/aspect(main+0x6b5)[0x55b5d02eef8e]
[quser42:2814811] [10] /lib/x86_64-linux-gnu/libc.so.6(+0x29d90)[0x7f0a4ecfcd90]
[quser42:2814811] [11] /lib/x86_64-linux-gnu/libc.so.6(__libc_start_main+0x80)[0x7f0a4ecfce40]
[quser42:2814811] [12] /opt/aspect/bin/aspect(_start+0x25)[0x55b5ce449955]
[quser42:2814811] *** End of error message ***
Aborted

@maxqcollins In both cases, as I said above, run the program and see if you can get a backtrace that helps you pinpoint in which line the issue happens!

Best

W.

Just to clarify, are you asking in the input file which line is responsible for causing the segmentation fault? It occurs when I try to read from the data files “pure_ice.txt” and “saline_ice_940.txt” for the compositional fields.

@maxqcollins No, which line in the ASPECT source code. That may give us an idea in which way the input file is wrong or read incorrectly.

Best

W>

I am not exactly sure how else to run aspect in debug mode. The previous runs I posted were all run with aspect-debug.

@maxqcollins As I mentioned, running ASPECT in a debugger is what will help you narrow down what the problem is. If you’ve never done that, you may want to take a look at lectures 25 (and perhaps 41.25, if you prefer working on the command line) here: Page Not Found There are also many other resources out there that show how to run a program in a debugger to query the state of a program and find out why it crashes. You will find that there is a learning curve associated with using a debugger, but from experience I can guarantee that you will in the long run consider it a major boost of productivity if you know how to use debuggers.

Best

W.

@bangerth Thanks for the direction! It seems like the link isn’t working, but I’ll do some googling myself as well.

@bangerth This is the backtrace:

(lldb) bt
* thread #1, queue = 'com.apple.main-thread', stop reason = EXC_BAD_ACCESS (code=1, address=0x17)
  * frame #0: 0x0000000100047ed8 aspect-debug`std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char>>::size[abi:ne190102]() const [inlined] std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char>>::__is_long[abi:ne190102](this=Summary Unavailable) const at string:1881:29
    frame #1: 0x0000000100047ed8 aspect-debug`std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char>>::size[abi:ne190102](this=Summary Unavailable) const at string:1285:12
    frame #2: 0x0000000100167ca0 aspect-debug`std::__1::basic_string<char, std::__1::char_traits<char>, std::__1::allocator<char>> std::__1::operator+[abi:ne190102]<char, std::__1::char_traits<char>, std::__1::allocator<char>>(__lhs="/Users/maxcollins/Documents/Research/Ice Dynamics/aspect_updated/data/material-model/steinberger/", __rhs=Summary Unavailable) at string:3994:25
    frame #3: 0x00000001009c1964 aspect-debug`aspect::MaterialModel::GrainSize<2>::initialize(this=0x000000012730a250) at grain_size.cc:97:37
    frame #4: 0x000000010050e620 aspect-debug`aspect::Simulator<2>::Simulator(this=0x000000016fdfab90, mpi_communicator_=<unavailable>, prm=0x000000016fdfa420) at core.cc:311:21
    frame #5: 0x0000000100b8a47c aspect-debug`main [inlined] void (anonymous namespace)::run_simulator<2>(raw_input_as_string="#########################################################\n# This is a variation of the composition-passive.prm file.\n# Here, we choose a higher Rayleigh number and make the\n# density depend on the compositional variable as well.\n#\n# See the manual for more information about this setup.\n\n\nset Dimension                              = 2\nset Start time                             = 0\nset End time                               = 1e6\n# set Maximum time step                      = 100e2\n# set CFL number                             = 0.1\nset Use years in output instead of seconds = true\nset Output directory                       = output-940-grainsize\nset Resume computation\t \t\t   = auto\n\n# Issues with converging if too large \nsubsection Time stepping\n  set Minimum time step size                 = 0.5\nend\n\nsubsection Geometry model\n  set Model name = box\n\n  subsection Box\n    set X extent = 60e3\n    set Y extent = 30e3\n  end\nend\n\n############### Parameters describing the temperature field\n# We ignore any temperature e"..., input_as_string="#########################################################\n# This is a variation of the composition-passive.prm file.\n# Here, we choose a higher Rayleigh number and make the\n# density depend on the compositional variable as well.\n#\n# See the manual for more information about this setup.\n\n\nset Dimension                              = 2\nset Start time                             = 0\nset End time                               = 1e6\n# set Maximum time step                      = 100e2\n# set CFL number                             = 0.1\nset Use years in output instead of seconds = true\nset Output directory                       = output-940-grainsize\nset Resume computation\t \t\t   = auto\n\n# Issues with converging if too large \nsubsection Time stepping\n  set Minimum time step size                 = 0.5\nend\n\nsubsection Geometry model\n  set Model name = box\n\n  subsection Box\n    set X extent = 60e3\n    set Y extent = 30e3\n  end\nend\n\n############### Parameters describing the temperature field\n# We ignore any temperature e"..., output_json=false, output_xml=false, output_plugin_graph=false, validate_only=false) at main.cc:632:32
    frame #6: 0x0000000100b8a46c aspect-debug`main(argc=<unavailable>, argv=<unavailable>) at main.cc:861:13
    frame #7: 0x0000000197696b98 dyld`start + 6076

It seems like the problem is coming from grain_size.cc at line 97 concatenating the directory and filenames, but the filenames are listed in the prm file and are located in the default directory -

    template <int dim>
    void
    GrainSize<dim>::initialize()
    {
      CitationInfo::add("grainsize");

      n_material_data = material_file_names.size();
      for (unsigned i = 0; i < n_material_data; ++i)
        {
          if (material_file_format == perplex)
            material_lookup
            .push_back(std::make_unique<MaterialModel::MaterialUtilities::Lookup::PerplexReader>(datadirectory+material_file_names[i],
                       use_bilinear_interpolation,
                       this->get_mpi_communicator()));
          else if (material_file_format == hefesto)
            material_lookup
            .push_back(std::make_unique<MaterialModel::MaterialUtilities::Lookup::HeFESToReader>(datadirectory+material_file_names[i], # ****here***
                       datadirectory+derivatives_file_names[i],
                       use_bilinear_interpolation,
                       this->get_mpi_communicator()));
          else
            AssertThrow (false, ExcNotImplemented());
        }
    }

Thanks,

Max

@maxqcollins This is helpful. In the loop, you are accessing material_file_names[i], which is fine because you are running a loop of the form

      n_material_data = material_file_names.size();
      for (unsigned i = 0; i < n_material_data; ++i)
 

but then you are also accessing derivatives_file_names[i], for which it is not clear that there are enough elements in that array. In the debugger, can you go to frame #3 after the program crashes and print both i and derivatives_file_names.size()?