Unable to update the densityField through array manipulation

[bknight1]

Hi team,

I'm trying to load in data for density derived from Vp velocities.

I am usually able to update the swarm or mesh variable array through array manipulation, e.g. Model.materialField.data[:] = crust1.index, works fine for updating the material field.

However, when trying to update the density field using the same method (Model.densityField.data[:,0] = 0.) I am unable to change the values within the densityField array.

I have added in an additional swarm variable and I am able to import the density data that way so the ordering of the data is fine.

Cheers

[rbeucher]

Hi Ben,

I think that is the correct way to do things. You need to revert back to some core underworld functionalities for that. An additional swarm variable is good. You can define the material density to be that variable.

R

[bknight1]

Hi @rbeucher

Thanks for that, whats the best way of assigning the new density variable to the density/buoyancy function?

Cheers

[rbeucher]

Hi @bknight1 ,

I have added the option to define the densities as UW functions. (see development branch. This will be folded into version 2.11 once we have a release) You can also do what you were describing but accessing the density field via Model._densityField.data...

I think @NengLu should be able to help if needed.

I am closing this for now. R

[bknight1]

Thanks @rbeucher

I did try the Model._densityField.data approach but it wouldn't update the density field on v2.10. I did manage to get it to work by converting density to temperature. I'll give it a go on v2.11.

Cheers

[bknight1]

Hi @rbeucher

I'm interested in this again as I want to include some thermodynamic modelling calculations into my models to replicate a more realistic density with the geodynamic models. Once I have the density working I'll be looking at other parameters that could be implemented.

I have a PT space with a corresponding density (fig below) created using perpleX. I have then created an array for the density based on the pressure and temperature at each particle. I then set the Model._densityField.data in pre and post processing loops to the calculated density array, however this does not appear to alter the density used in the stokes solver, which (I believe) is linked to Model.densityField, which cannot be updated (?) through Model.densityField.data[:,0] = ...

If you have any ideas please let me know.

[rbeucher]

Hi @bknight1,

In UWGeo you will have to update the density per Material. Have a look at the PhaseChange Class that is used on a Material.

Your PT space could be used as a condition. But you will have to Add it to the materials. Technically the Model class inherits from Material so a PhaseChange could be applied to the entire model.

Let me know how you go. I'll be interested to see what you come up with.

Romain

[rbeucher]

I don't think that changing the values via the .data handler is a good. idea. Better to define a function that is applied to the Model.density. You can use a condition based on P,T.

[bknight1]

Thanks for the help @rbeucher

I'll look into it over the next coming days and see if I can come up with something!

Cheers

[bknight1]

Hi @rbeucher

Quick update, I made a load of new materials which corresponds to each unique density for each material I want to perform changes in density (e.g. crust and mantle). I then modify the materialField based on the P-T to give the material index that corresponds to that density.

I'm not sure how many different materials have been used before? I have ~140 (~90 for the crust and ~50 for the mantle) which appears to slow down the Model.init_model() and Model.run_for(...) significantly.

Cheers

[PatriceFRey]

Hello @bknight, the Earth isn't in thermodynamic equilibrium, that's why we find granulite, and fresh peridotite at the Earth's surface. Adapting densities based on a thermodynamic model makes sense at high temperature, but perhaps not at temperatures prevalent in the mid to upper crust. Cheers

[bknight1]

I agree @patrice-rey, my aim is to develop a workflow so whoever is running the model can input their own thermodynamic model and which materials to apply it to. I am not a petrologist so I can't comment on the PT pseudosections, that's why I aim to create a workflow for others to be able to import their own. I recently came across a paper that used a lookup table to modify the density as the orogenic wedge evolved (https://doi.org/10.5194/se-12-1749-2021) so I am trying to see if this method can be used in UW. Cheers

[PatriceFRey]

Sounds great @bknight1. Thinking ahead ... temperature, fluid and strain rate allow rocks to change their mineralogical composition. Putting aside fluid, having temperature and strain rate conditions would be very useful. For instance densities of the lookup table could become progressively important as temperature increases from say T <= 250ºC (i.e. densities remain unaffected by that in the lookup table) to T >= 750ºC where densities of the lookup table are fully realised. In addition we could also add a condition on strain rates as we know that strain rates facilitate greatly thermodynamic equilibration. As @rbeucher did mentioned earlier, we can use phase changes to introduce key density variations (eclogitisation, partial melting, etc), but I see how a lookup table is a much better approach. Best of luck!

[rbeucher]

I expect it to be very slow. Maybe we should use the mesh instead of the particles...

[bknight1]

I am always looking at implementing new stuff @patrice-rey so would be interested to look into it in the future!

I have a fairly fast implementation @rbeucher, the code needs a bit of tidying up but I'll share what I have implemented so far in the coming days. I project the temperature and pressure onto the particles, but you might right that using the mesh to do the lookup then project the materials from the mesh to the particles might be a better approach.

[bknight1]

I have invited you @rbeucher to a private repo where the files are hosted. At the moment its just a PoC to show my implementation. Keen to hear if you have any feedback on the workflow.

Cheers

[rbeucher]

Great. I'll see if I can have a look asap. I'm teaching this week but will do my best

[bknight1]

There's no rush @rbeucher I've got some other stuff to work on for now so all good!