How to calculate the thickness of a specific compositional field?

[Xie]

Hi all,

In ASPECT, compositional fields are used to represent geological layers in the real world, and their initial thickness is limited by manually input parameters in the ‘Subsection Initial composition model’. However, in most cases, both the layer’s shape and thickness will change over time, for instance, the lithosphere is deformed by compression. Therefore, i would like to know how to calcualte the thickness at different timestep of a specific compositional field. Ideally, use ‘Postprocess’ to output the results for subsequent calculation.

Assuming a compositional field named ‘lower_crust’, I think 'get_geometry_model().depth(position) ’ should be used, it seems that the key point to this problem is how to accurately locate the ‘position’ in the lower_crust, so, any suggestions?

Best,
Xie

[jperryhouts]

I think this is probably going to be complicated to do in a postprocessor. I’d suggest either: add flat layer of passive particles along the interface you’d like to track and analyze their locations over time or do this postprocessing via Paraview with an isosurface of “lower_crust” or something similar.

[bangerth]

Xie,
I don’t think there is anything that’s already done. You’d have to loop over all cells, evaluate the compositional fields at points on each cell, determine whether the value of the compositional field indicates that you are in layer X, and if so determine the depth of that point.

Then you have to decide what you want to do with this information. The question is how exactly you define “thickness”. Is it the maximal depth over the entire domain? Or is it a quantity that is different for every lat/long point?

If you answer this question, you also still have to make this work for the parallel case by doing a max/min/average over all processors.

Best
Wolfgang

[jbnaliboff]

Hi Xie,

I agree with Jonathan and Wolfgang that this will not easily be done in a postprocessor. Exactly how you want to define the ‘thickness’ of each layer may change through time, and you do not want to rerun models simply to reprocess the results.

Rather, I suggest importing the data into python or Matlab and then writing a script to calculate the layer thickness. I typically do this by evaluating which compositional field representing rock types has the largest values at each point. From there, you can do things like find the minimum and maximum depth at each horizontal and/or longitudinal position.

One can import solution data into python or Matlab by exporting data from Paraview in a csv file format, or using the vtk python package. I’ll make a separate post about the latter method.

Cheers,
John

[jbnaliboff]

Hi Xie,

See the following post for importing ASPECT VTU files into numpy arrays:

Loading VTK data into Python

Hi all, This topic has come up a number of times in recent conversations, so I thought I would share the general workflow for importing VTK data directly into python (numpy arrays). Below is a snippet of python code that should largely work out of the box for 3D cartesian models, where all one has to do is provide the full path and name of the .pvtu in a string. Please note that lines in bold are comments, and you would need to add a # symbol before them in python. I excluded them as they prod…

Cheers,
John

[Xie]

@all,

Thanks for your reply. As for the definition of ‘thickness’, i’d like to define it as the distance between the upper and lower boundary of a compositional filed corresponding to the same X-coordinate.

I think it would be meaningful to do so, because in some cases, for example, some criterions only applicable to some certain layers, such as applying Byerlee’s law to the lithosphere.

I’ll try what John mentioned.

Best,
Xie

[bangerth]

Xie,
this definition (distance between top and bottom) is something that is easy to compute on a uniform lat/long/depth mesh stored on one processor, because you can walk down in depth direction from every point and just check for where you enter and leave a particular layer. But it illustrates the difficulty with more complicated mesh data structures:

• In parallel, each processor only stores part of the mesh. It is not guaranteed that the processor that has the top part of a layer also has the bottom part. Some processors may have no part of the layer at all.

• Let’s assume that you have found a grid point at a particular x-coordinate that indicates the top part of the layer. It is not guaranteed that at that x-value there is another grid point at all in unstructured meshes, and even if there is, it may not be where the interface is.

Like John already suggested, the way to deal with these sorts of issues is to work in a viz tool like Paraview, extract isosurfaces of the compositional field that correspond to layer boundaries (e.g., “compositional field for layer X has value 0.5”), and extract those in the form of CSV files. Then you want to postprocess them in python or matlab somehow.

Best
Wolfgang