Reading an ascii file as initial composition

ASPECT
1

Hello, when I give ASPECT an ascii file as an initial composition for density, the density is interpolated at the quadrature points; ASPECT finds the closest density in the provided dataset at the cells boundaries only and then makes the interpolation onto the quadrature points:

ASPECTascii-interp

This means that the resolution of my domain has to overfit the resolution of the dataset that I provided to be able to catch (to some extent) details of the dataset.
Therefore, would it be possible to have interpolation alternatives so that the setting of initial conditions from reading an ascii file is similar to the setting of the initial conditions using a material model?

ASPECTascii-alter

At the moment, there is a subsection interpolation scheme in the initial composition subsection but this only concerns what the value at the cells boundaries is; there is still an interpolation onto the quadrature points.

And just to be sure, if I use a material model (for instance to have PREM density in my domain), the density would then be written directly at the quadrature points, right?

Many thanks,
Ludovic

2

That’s not quite true: When you interpolate, it evaluates the density at
the node points of the finite element (which happen to be the vertices
if you use a Q1 element) and then these values define the finite
element field because the finite element field is piecewise polynomial.
In your alternative scenario, you want a field that happens to have a
kink – but that’s not polynomial inside the one cell, and you can’t
represent it as a finite element field.

So far, there are no quadrature points. Only when you use this
interpolated compositional (or temperature) field in other equations is
the piecewise polynomial evaluated at the quadrature points. I think
that what you want is to not evaluate the (interpolated) finite element
field, but instead to use the original initial composition as
represented by the ASCIIData file. But that’s an awkward choice, because
you could only possibly do that in the very first time step: After that,
you only have the current compositional field, whereas the initial
conditions no longer apply.

Best
W.

3

Thanks Wolfgang for your reply and clarification of my messy words…

That actually makes sense and I forgot about that. The thing is, I do instantaneous modeling and I do not need to solve Stokes.

That’s clearer, but what happen if I define a new material model where I specify a certain density with depth (for instance PREM), does the density also evaluated only at the nodes points of the finite element?

Here my problem and experiment: (still work in progress)
I have an ASCIIData file of 167 layers of density values PREM and I try to figure out why I get an excess of mass of my domain (hollow sphere filled with PREM density) compared with the analytical solution. At the moment I can reach a minimum mass error of 0.06% (it seems nothing but in terms of gravity prediction or the Earth, it is quite a lot.

I did the following experiment: I created a customized hollow sphere using a list of radial values to have independent layers for each region described in PREM. When I use the ASCIIData file in Q1 and I output the density for my topmost layer (3km), I have the following densities: 1198, 1910, 2422 from top to bottom. The 4 first layers of the ASCIIData file are every km (1020, 1020, 1020 and 2600). So I understand my output: during the setup of the initial condition, ASPECT looks at the layer 0km (density 1020) and 3km (density 2600) and these values only are used to define the finite element field.
However, when I am in Q2 (compositional field), I have densities 882, 1020, 2105 for the same layer. And that I do not understand because the densities go below the minimum density of the dataset.

From this single region cells customized spherical shell, I also refined the mesh 1 to 3 times - I would expect to come closer to the theoretical mass but no, I still have a this 0.06% error.

Finally, In the gravity_point_value postprocessor, I did the experiment of calling the material model and assign the PREM value directly at the quadrature points. I performed tests by increasing the resolution of the domain and I have a convergence towards the analytical mass (the error converges to 1.7e-5%)

There are 2 reasons why that could happen: the ASCIIData file I created is somewhat wrong (I continue checking it). Or this is about the interpolation during the evaluation which cannot represent the details of the dataset?