There is a problem when use "Initial temperature model.Adiabatic.half-space cooling" with number of composition greater than 1(ascii model)

Hello everyone,

As the title, when I use

subsection Compositional fields
  set Number of fields = 2
  set Names of fields  = continent, ocean
end

subsection Initial composition model
  set List of model names =  ascii data
  
  subsection Ascii data model
    set Data directory = ./Gplates_seafloorAge/
    set Data file name = Double-Composition140Ma_64800.txt
  end
end

And the ascii data is

# POINTS: 4 360 180
# r phi theta proportion_material_1 proportion_material_2
3471000	0.008727 	0.008727 	0.00 	0.00 
3471000	0.026180 	0.008727 	0.00 	0.00 
3471000	0.043633 	0.008727 	0.00 	0.00 
...

with the initial temperature subsection as:

subsection Initial temperature model
  set Model name =  adiabatic

  subsection Adiabatic
    # set Age top boundary layer = 5e7
    set Age bottom boundary layer = 10e7
    set Cooling model = half-space cooling
    set Lithosphere thickness = 120000
    set Top boundary layer age model = constant
    set Data directory = ./Gplates_seafloorAge/
    set Data file name = SeafloorAge140Ma_64800.txt
  end

this will be wrong with the following information:

ERROR: FunctionParser failed to parse
	'Initial temperature model.Adiabatic.Function'
with expression
	'0'More information about the cause of the parse error 
is shown below.
---
An error occurred in line <80> of file </public5/home/sch8017/software/dealii-9.4.1/source/base/function_parser.cc> in function
    void dealii::FunctionParser<dim>::initialize(const string&, const std::vector<std::__cxx11::basic_string<char> >&, const std::map<std::__cxx11::basic_string<char>, double>&, bool) [with int dim = 1; std::__cxx11::string = std::__cxx11::basic_string<char>]
The violated condition was: 
    this->n_components == expressions.size()
Additional information: 
    The number of components (2) is not equal to the number of expressions
    (1).

It seems wrong with dealii source code not ASPECT. But I don’t know how to change the code to make it running well with “Adiabatic” and “two or more ascii composition fields”, and I can’t find why this will work with one composition field but trouble with two composition field.
Can anyone give me some suggests? thank you very much !

Dear Shaohui,

The adiabatic initial temperature model has a Function subsection, with the function expression set to zero by default. This is needed because the material properties needed to compute the profile (density, thermal expansivity etc.) depend on the composition, so we need to let the model know what values to use for the composition (i.e. what is the composition along the adiabatic reference profile?).

By default, this function expression is set to zero, i.e.

subsection Function
    set Function constants = 
    set Function expression = 0
    set Variable names = x,t
end

Since you have two compositional fields, you also have to provide two components in your function expression rather than just one, i.e.

subsection Function
    set Function constants = 
    set Function expression = 0; 0
    set Variable names = x,t
end

(if 0 is a reasonable assumption for your average 1d profile.)
You have to add this as a subsection into your subsection Adiabatic.

This is what the error message is trying to tell you when it says The number of components (2) (you should provide) is not equal to the number of expressions (1) (which is what is currently in the input file).

Best, 
Juliane

Dear Prof. Juliane Dannberg,

After following your suggestion, it was solved very well.
Learning ASPECT independently will always encounter some problems, thank you very much for your help.
I recently read your article submitted to GJI “Changes in core-mantle boundary heat flux patterns
throughout the supercontinent cycle”.
I saw that you used the particle method not field to track the composition, why do you choose this?
I used to use the “particle” but I found that " Maximum/Minimum particles per cell" is useless, I’m confused about that.
If it’s possible, could this pre-print paper’s .prm file be public for reference? I want to find some 3D-spherical .prm files to learn. It also doesn’t matter if it’s inconvenient to make this file public, I can also try to find other .prm files to reference.

Sincerely,
Shaohui Liu

Hi Shaohui,

Response to your questions below, but Juliane or others may also weigh in.

I saw that you used the particle method not field to track the composition, why do you choose this?

Tracking composition (lithology) and other properties on particles minimizes the amount of diffusion relative to tracking material with compositional fields (although using DG fields + limiters can help).

One drawback is that a very large number of particles per cell is needed in 3D to accurately track interfaces, but recently added particle interpolation schemes (bilinear and quadratic least squares + limiters) allow using a tractable number of particles per cell in 3D for most problems. This is work in progress and I believe a paper is forthcoming, but those features are currently available in the main branch of ASPECT.

I used to use the “particle” but I found that " Maximum/Minimum particles per cell" is useless, I’m confused about that.

In what way is this feature not working? Is the min/max particles per cell not being enforced, or are your simulation results not satisfactory with this feature enabled? If the former, is the following parameter specified in the particle subsection of your prm file?

set Load balancing strategy = remove and add particles

If it’s possible, could this pre-print paper’s .prm file be public for reference? I want to find some 3D-spherical .prm files to learn. It also doesn’t matter if it’s inconvenient to make this file public, I can also try to find other .prm files to reference.

There is a link at the bottom of the Arxiv version of the paper to the following zenodo repository, which contains the prm files and other data required to reproduce the simulations.

Cheers,
John

Dear John,

Thank you for your reply! You are really patient and kind.
I have two more questions, could you give me some advice.
First, whether should I add “weak zone” in the subduction area in my 3D spherical model, it’s difficult for me now. But I saw that the model without weak zone, the subduction shape will be not normal. Or yield stress should be add instead.
Second, I want to track the “sediment carbonates” and “CO2 degassing” in my future work. But I’m confused now. Should I use the “Steinberger model” with Perple_X to track the thermo-mechanical CO2 degassing in my model? Could you please give some advice about achieve this in ASPECT? It will be useful for me.
Looking forward to your reply.

Cheers,
Shaohui

Hi Shaouhi,

Happy to help and hopefully the particles in your simulations are not working as intended. With that said, both of your questions above are quite difficulty to answers.

Indeed, most models of subduction (instantaneous or long-term) do include either a weak zone, yielding, or both features. To decide how best to proceed, I recommend examining setups from previously published models that investigate similar subduction processes to those you wish to study.

Tracking Sediment carbonates and CO2 degassing - Depending on the level of complexity and intended goals of the simulation - that will likely involve quite a bit of underlying development work, even if it is just tracking CO2 evolution (i.e., no reactive flow). We would need to know quite a bit more about your specific goals to advise further.

Cheers,
John