Tutorial: Mandyoc, a finite element simulator for the mantle

live streaming is on all right wonderful well welcome back to uh this next installment of transform 2022 um this is a sort of festival of digital science or something like that um subsurface science and engineering we've got a whole host of tutorials lined up for you today or in this session we're going to be hearing about a finite element modeling tool for geodynamics so i'm pretty excited about that because i know very little indeed about that topic uh either of those topics actually so i'm uh very excited to have um averstina here who's gonna take us through that in a minute first i just want to mention that we've got lots more content coming up this week so check out software underground.org transform to see sort of what's going on um the chat is all happening in slack there's a slack group uh so if you want to take part in that not just during the tutorial but also afterwards and for the next weeks and months um we'll be chatting about this topic and others in software underground.org slack so that's completely free to sign up um you could just come and look you know you don't uh you don't have to be a vocal person um but just come and find out all the cool stuff that's going on in this space um what else did i want to say oh yeah i mean if you're watching live don't forget you can like dvr youtube so you can pause you can rewind uh you can even slow it down i think um so yeah take advantage of that the youtube video will be up for the foreseeable future so uh enjoy and share please and also subscribe to the channel if you want to get notifications of the next stuff that comes up um that's it that's all i have to say i'm going to hang over to augustine who's going gonna take us through this uh interesting tool so over to you know welcome hello let me share the screen to start with this tutorial great okay we'll start so uh thanks matt for the presentation and thanks for underground to let me to do this tutorial today in this opportunity i would like to present you a showdynamic code called mandioc and a showcase to simulate the evolution of the divergent machines using this code so i will explain you how to create an input file to make this simulation how to run the model using mandioc and plot some of the results for that i create um this i create this repository you can find the all the material here in this repository of github let me show you here you can find the same notebook that i will to do today you can this is more complete than the the notebook that i will show you today and here you have all the information about how to install the software and everything so let's go back to the notebook so but first to but first let me to introduce myself and the others instructors for this tutorial my name is authena pese from argentina i do i doing my phd i work in application physics processing gravity and magnetic data i am recently finished in my post-doctoral research where i studied subduction songs through numerical modeling using mandioc and contribute to development development of this software i am the main speaker today but you can find the other instructor in the slack channel answering your doubts they are viktor sasek shane sona sun sao and rafael montero da silva uh i would like to point out that the development of mandioc began some years ago by the doctor victor zasek and the institute of astronomia geophysica la universidad is paulo in brazil but in the last years we dedicate to creating the documentation the examples and improve everything to publish this code now mandioc is open source code in github and [Music] we up and finally we have published a paper in shows about this this code so let me show you here is the repository of mandior here you can find everything to run a simulation like i will show you in some minutes and you can find the documentation websites of of course with everything that you need so if you want to know more about this i you can visit this website and ask of course in slack so before to start what is mandioc mandioc is the aggrams of mandiot dynamic simulations code and it is designed to simulate the stoic flow type thermo chemical convection taking different compositional layers into a ground and it is a also a appropriate to simulate earth lithospheric dynamics on hyeologic theological times time scale sorry the co the code is fully parallelized using pepsi library different linear unknown linear rheologists can be adopted to simulate the strain and the stress pattern in the earth's crust and mandel both in extensional or conditional tectonics scenarios one advantage of mandiog is the possibility to create scenarios with velocity boundaries conditions variables in time and space a allowing the user the user to simulate different tectonics pools in the same uh in the same model additional the current version incorporate surface process imposing rate of erosion and sedimentation on the top of the free surface to simulate mandrel thermodynamics thermochemical sorry a convection mandioc abducts the formulation of a newtonian fluid together with the boussinesq approximation to solve the equation of conservation of mass momentum and energy these three equations that we can see here but let me show you to this picture to summarize the step that mandiox takes to solve the conservation's equations and perform a simulation this picture shows that once the goals start running and the input file are read but mandiog used the effected viscosity field to calculate the velocity field and check if the comprehensive conditions satisfies a given value while the minimum tolerance is not right mandio used the osawa method to integrate indirectly calculate a new velocity and pressure fields the updated fields modify the viscosity field which perturbs again the velocity field so this field will be updated until the tolerance is reached the compositional strain history is struck by particles presented in the interior of the finite element the exchange of particles almost in the sub along the subdomains of the model is efficiently paralyzed using betsy functions the free surface of the earth can be simulate an aesthetically numerically numerically stabilized using the free surface stabilization algorithms so here i will show you three papers that will that that they use mandiox to solve the the problems the economics problems so in the repository you can find the links to these papers so i hope i hope you have an idea of how mandiog works i know that is a little difficult to understand and the first in the first moment if you want to learn more about this code you can visit the mandiog website or read the paper that we publish in shows or ask in the slack channel of course so in this moment i am going to introduce you introduce you a particular model that can be made using mandioc but in the maniac repository there are another examples er to start to play with a mandioc and making simulations eh so let me show you this if you go to the repository of mandior and go to the examples webs example folder there you can find another examples that you can perform as examples and [Applause] play on how to make a similar make a model and make a simulation using mandioc so in this tutorial we want to simulate the evolution of diversion machines taking into account a plastic rheology and the syn roof shear dynamics i am going to explain you how to generate the initials the initial setup to obtain the resulting simulation that you can see in this gift so um there are a few ways to run this simulation the ones that i am going to present here kind came from the association that the geometry temperature and velocity at the beginning of the simulations are known because of that we must prove the velocity the temp the temperature and the geometry of the model as input files these input files are the interface with the boundaries between these different layers of the model the initial temperature file with the initial temperature fill and the initial velocity file obviously with the initial velocity field so the idea is to generate the input file that create this setup you can find different layers this is the profile of temperature that we want to generate and the velocity we only have velocity horizontal velocity in the borders of the model and you can see of course the the the domain of the model so first of all to do that we need to import some python package they are numpy pandas matplotlib and warnings only to turn off any warnings that we can we can generate in this in this notebook and also we will import some functions in this part this function that these functions was created by me only to make this posit this this tutorial more gradable so we need to run the cell and now we need to define the domain of the model in this case this model uh comprise a 1600 by 300 kilometers squares so we define ls sorry l x equal to one sig zero zero e three because we need we want everything in meters and l c equal to three zero zero e three and and enter and run the cell now we need to discretize this domain in element of one by one kilometers for that we we define the number of nodes in both axis as nx equal to 1601 and xc equal to 301 and run the cell now we need to generate the a point of the of the axis because we want to generate a mess for that we want to use the number line space a function so x equal to mp dot line space and put here we want um an arrive that go from go from goes from 0 to a lx with nx numbers value sorry and we'll generate c equal mp dot line space from minus um lc to 0 with nc numbers so now we run this and now generate the mercury they met the the mesh grid you see in the numpy function so now that we have the domain define it we can generate the input files so for this simulation the domain is a half different lithological units these units are the aston offer the mandelita the lithospheric mantle we have here a mantle seed a lower cross upper cross and the sdkr this is how they are our six different layers the weak seed or the mantle seed that we put there in the in the in the in the picture is a small layer into the lithospheric mantle to ensure the nucleation of the wrist team at the center of the numerical model so now before to start to generate a interface file we need to define the thickness of the these different layers for that we are going to define a h underscore s higher equal to 40 kilometers eight underscore um lower grad underscore trust equal to 20 kilometers like sorry here i have a typo 20 20 kilometers so we need to define the thickness of the upper grass upper grass as h underscore upper underscore crust equal to 20 kilometers and all the little tree will have and will has an thickness of 130 kilometers so age underscore little equal to 130 kilometers and here we are we define the location of the mantle seed this will be a 30 kilometers below the lower cross like we can see in the in the picture so run the cell and now we need to generate the interface of these layers um believe me when i tell you that this cell create a dictionary called interfaces which contain the depth between the different layers youth and the sickness that we defined before so here only generate that a dictionary with the different layers and use the the thickness to generate the depth of the different um interfaces so i will run this one and here only in this ah sorry i have a typo sorry p i have a typo in the name of the thickness of the little spray so i run again the cells so great now in this cell that we have here i'm only plot the different interfaces using f using a for loop for that and of course must plot leave so if we run and everything are in kilometers so you can find here you can see these interfaces the different interfaces on the depth of the interfaces here in blue we have the interface of the astronomer and the lithopheric mantle in red we have the interface between the upper lower crust and the lithospheric mantle this one is the interface of the upper and lower crust in brown the interface between the air and the upper crust here in green and the yellow interface are the interface of the weak layer seed you can see that they have the same the same depth because we define this in the interfaces dictionary in this moment we want to we are going to create this a this mantle seed that we don't have defined yet so for that we want that the these um this is seed will be in the center and the center of the model so the horizontal position will be um will be 1750 kilometers and this seed will have a thickness of six kilometers so i will run this and here in the cell and only use a conditional statement to add the signals to the different uh interface that generate the the weeks the weak layer so but i if i run this um plot again use the same cell that i used before to plot the interfaces again now you can see that we split the the interface the seed top and seed base for the for the for the seed mantle so we generate the interface for the for the the mantle seed so now we have all the interfaces created so now we need to define the properties of this layer of these different layers to generate the interface file because i think that i never said that the hindi sorry the interface files contain the depth between the different layers and the compositional properties of these layers sorry so to create the the to define the properties of this layer we are going to adopt the viscoplastic rheology when where the effective viscosity combines non-linear power law biscuit rheology and a plastic yale criterions so the biscuit component is given by by this equation so the different properties that we need to specify in the interface files are the compositional factor the density the rheological heat the exponential scale factor the power load component the activation energy and activation volume so in the table that you can see here we have um all the values that we are going to use for for each for all the interfaces so here are the parameters and in the vertical we have the parameters and in the original we have the names of the different layers in addition to simulate the free surface we adapt the stick air will be you use the stick air model so the 30 kilometer for the iron layer will have a low viscosity and density similar to the asthenospheric air here uh so taking this on because of ah sorry for so now here in this cell we need to define that then we are defined the density of the upper grass and the lower grass and generate the radio heat for these two layers please remember these parameters because we are going to use them to to generate the temperature file so now here you can see that we only transform the table into an string called layers properties because in this other cell that we run here we generate an interface file called interfaces dot txt that contain this the properties the layer properties here and here we save the depth between the different interfaces so if you go to the repository folder you can find find here the interface file so in this file you can see here the difference properties of the layers and here that that f uh the interface depth between the layers here is for the astronomer and the up the lithospheric mantle this one for and this one these two this one for the stick for sorry for the mantle seed this one is between the lithospheric mounted lithophoric mantle and the lower um lower grass these between the upper and lower grass and here between the air and the upper crust so so now we generate the interface files the next step is to generate the input temperature file for that the initial the initial temperature structure in the lithosphere is death depends and it's zero degrees celsius for in the surface and um once a thousand and three hundred kill a degree celsius at the base of the little tree at a 130 kilometers and the sublithospheric temperature will be there the adiabatic temperature that is given by this equation so before to the define the temporal structure here in our model uh we need to define some parameters they are the thermal divusity the kappa the c capacity is the specific heat capacity temp underscore p is the potential temperature of the mantle she is a gravity value that we hear put them and alpha is the volumetric expansion expression confusion coefficient and so now we are defined the structure and the temperature structure so here in this line we are defining the a linear temperature for the asteroid fair sorry for the lithospheric for all the lithospheric here we are defined the adiabatic temperature we calculate the adiabatic temperature for the asthenospheric here we put zero the value of the all the nodes over the into the higher layer and here we only combinate we only add the temperate adiabatic tempura in the in the temp uh in the temperature uh grid that we generate at first so we add the adiabatic temperature in the node that are located in the asteroid fair so now we so now we need to run we need to apply the 1d diffusion equation in the model for but before to do that we need to create the initial heat production model or so because of that we need to initialize h as zeros values for all the nodes in the grid but we change these values of zero for h uc and hlc because we have we want that the nodes into the upper cross and lower graphs have another values different to zero and these values are the values that we use in them to generate the interface file so i am in this two lines i do that i put these values into the h m h m a array so grind these cells and now we can apply the fusion equation for ten thousands of years using an equation that i create for this tutorial you can if you want to know more about this function you can go to their tutorial repository and there you can a read about this equation and and if you have any doubt you can ask in a slack channel of course so now we apply this equation and in this cell i only plot the thermal structure under and a thermal profile so because i that i use a contour f and a maple function to plot the temperature and use um and then i plot the temperature profile and the interface because it's more easy to see the change of temperature so here you can see this one is the temperature for the lithospheric structure and here we have the adiabatic lithospheric values for the aston offer and this one is the temperature for all the domain the structure the temperature structure for all the domain and we have for the all the depth over zero have a temperature of zero so now we only save the structure temperature into any files called input underscore temperature underscore zero point txt in all so if we go to open the initial temperature file here you can this file includes the temperature or each node starting at the bottom left of the model and going up in that direction so here we have all the values of timber for the nodes so let me go back to the notebook so now the last the last file that we need to create is the initial number the velocity file and for that we need to it's more easy to understand this file using this picture you can see that the we have only velocity in the borders of the model in the rest of the note of the model we don't have velocity so in this case the velocity field simulates the lithospheric stretching assuming a reference frame fixed on the lithospheric plate on the left side of the model and the plate on the right side moves a rightward with a velocity of one cent centimeter per year the velocity file in the left and right boundaries of the model is chosen to ensure ensure the conservation of mass and is symmetrical if it adopted reference thread moves to the right with a velocity of 0.5 centimeter per year relative to the left plate so taking this into account i need to say that in this cell i generate the velocity the velocity model here i only define some parameter this um i will don't show you these cells with a lot of detail because they are as very complicated to follow but here i define the velocities some parameters that we need to generate the velocities file and here initialize the velocities like zeros for all the nodes in the domains here use a conditional statement to generate the velocity for the left side of the model so i change the zero on the on the left side of the model for other value only for the horizontal component here i i do the same but in the right side um here i try and in this cell i check if the if the conservation is the inflow flow in the inflow on the vertical boundaries below the no sorry i check here that the inflow it's the same same value that the outflow only flow so i'm only going to prove this um and then if we have a some difference between the input and output flow i use i assign this value in the vertical velocity for the node in the for the node in the left side only that here so if you can see here i plot these velocities for the light for the left and right side and you can see that you can see that we don't have velocities in the vertical component component of the left and right side but we have different velocities for the left and right sides in the horizontal components and zero e for the i air one for now in the cell and only save this file these velocities this velocity structure as an input file called input input underscore velocity underscore zero point txt so now we have all the input files that we need to run this model with mandiog so in this moment i will make an stop and take 10 minutes as break and we can continue with this in third in 10 minutes if you want so let me put the clock here to take a break 14 minutes i will see you in 10 minutes okay thanks thanks aristina [Music] you can maybe turn off your video and your mic and then just leave it on that that's great love it okay i will [Music] the video i oh i closed the video sorry let me find with you perfect see you soon hey how's deana i'm just gonna i'll let you take it away when you're when you're ready okay looks cool by the way yeah thanks for using that i've seen that before thank you so we can generate the input files to obtain the the initial setup that i show you before show show you before in the picture so now we can run a simulation using manuk but prefused to do that we need to define we need to create another file called parameters file these files there are many parameters that we need to define and taking take into account when we're running a hydrodynamic simulation um depending on the type of simulation and you want to and you need to define some of these parameters or not sometimes this parameter file can be a a nightmare so in this moment i will show you this parameter file and explain some of these parameters that we need to define i only explain these parameters because they are the more important parameters that we need to define for this simulation in this case in this parameter file we need to define the number of lagrangian particles in each element so for that we need to put particles per element a underscore element equal to current a equal to 40 sorry we need to specify the maximum time time and a step for the simulation for that we define these two parameters the first one is step a underscores max equal to 4 40 for 40 south sorry 1400 and time underscore max equal to um see a million years one sorry 100 million years and we need to specify the temperatures from boundary conditions for the temperate the temperature values will be fixed on all the sides of the model during all the simulations and the nodes for the eye layer will be always zero temperature um we need to do the same with the velocity we need to specify the the velocity boundary condition in this case we'll adopt free free sleep on the top and the bottom of the simulator of the border of the simulation to avoid a asymmetry asymmetric lateral homogeneous model we adopt a random perturbation of the initial strain in each infinite element of the model for that we need to specify this a variable random underscore initial underscore strain equal to 0.3 and then we need to er predefine the real sheet set we need to set the predef finite rheology model to use during the simulation in this case we select rheology underscore model equals 9 because a underworld sorry a mandioc have a this type of rheology assigned in number nine this is for the viscoplastic reality that we define we when we create the the the interface file if you want to see more about these parameters or the others parameters that you can predefine in mandioc you can go to the magic documentations there you can find more information about it and and see the different examples that we create in the repository because depending of the simulation you have a different parameter file so this this big cell is the generate the parameter file you can see that we have a lotus variable here we here we define the geometry some simulations options that we need to define here we use the the number of particles here we can the time constant where we have the defined sorry define the step max enzyme max that we used before that we see before and here we define the sum the viscosity values that we are going to use in the simulation here is about the input and output files you don't worry easy to understand this you can find this information in the documentation these are we here define some physical physical parameters and and here we define the boundary velocity the velocity boundary conditions and the temperature boundary conditions again here that we see we use the same that we explained before and here we set up the relationship model equal to nine in the mandiogram documentation you can see the others types of rheological model that you can use with mandioc so if you run this um this cell you generate a parameter a file called param dot txt and if we open this file let me find here we have the information that you see before only that so now we have everything to run the model so to run a model you need to um to run the model you need to have mandioc installed on your pc in the tutorial repository or in the mandiog repository you you can find the all the instructions to install mandiox in your computers and explain they are how to run a simulation in this case to run this model that i create i need to open um a terminal sorry here i open the terminal and type in this case i will use a virtual a virtual machine so i need to log into this virtual machine sorry let me give me a moment i forgot to put my password oh and yes so i go to the tranform repository well because there i have all the files and s here we have all the input files and the parameter files so for run the model i need to type m mpe run minus sorry dash n equal to 2 because here i specify the number of cpus that i want to use for run this model and now i need to call mandiog mandiok and we here to set up two flags for mandiog here we put minus strain underscore seed and put 0.0 comma 3 comma 1.0 so really what are doing what are going to do this flag the first flag seed specify which layer will have an specific initial strain at this layer are the numbers zero and two for the aston offer and the mantle weak seat in the middle of the lithospheric mantle the second flag a strain seed specify the magnitude of the initial strain with respect to the each layer indicate in the previous flag in this case the asthenosphere will have no initial strain and the mantle seed width have an initial a constant strain of one the other layers will present a random perturbation of the value of the initial strain and the maximum magnitude of this initial strain is defined by the parameter random initial strain that we define in the parameter file so if i press enter now ah i have a typo here sorry it's mpe run no np mpe run so if i run star mandiak to process the initial files first and start to solve the stoke system give me a minute to start this take time is not so fast so they are start to calculate the velocities fields so i will stop here this emulation sorry i'm gonna stop the execution of mandior in these first time steps because the complete simulation can take many hours in in my case for this model this took me around 2 days using 8 cpu scores the time is relative to the complexity of the a particular scenario running 40 million years in time steps in a few tenths of thousands of years additionally this scenario has a nearly half million of nodes with 20 millions of particles tracking tracking the composition and the other physical properties so this is very complex i could have chosen a simple scenario for this tutorial more that we can to run in this tutorial but this type of model are not so geologically interesting because of that i prefer to do this model and then give you to you the homework to run this model in your computer when you have free time to see the the output so let me to come back to the notebook so you can you can't run the model in this moment but in the repository of the tutorial i put a folder called output done when there you can find some different files for different properties and for different timing steps that might generate for the simulations in computer so we use these files to plot the results and see the simulation and and see the temporal evolution of the model so you can see that the different file that mandus generate state the density the strain the temperature the time the velocity and the viscosity data and so to plot these file we only need to specify in this case the initial step that will be zero and the final step that will be fifty hundred and the step one 100 so here i use a for loop only to go to read the files and plot them to read the files i i used i use functions that i created for this tutorial you can find them in the repository tutorial so i only load the the time the density here the accumulative strain the temperature velocities and density file and then use mass plot leaf to plot the results for that i use contour f to plot the the density with the accumulative strain and he other plot that i will generate is the temperature and the velocity using weather and i also plot the viscosity using pc color mesh so if we run this cell now with generate five pictures five different um pictures with the indifference times this is four times zero million years here 8.5 million years fifth uh 16.9 million years 25 4 million years and 33 million years and the last one is for 20 22 million years so and then if you see the yeah we'll see the last one in the last one you can see that the dark and light orange color represent the upper and lower crust and while the dark and light green colors represent the lithosphere and the asthenosphere and the shadows of greed of green sorry the shadows of gray indicate the magnitude of the accumulated accumulative strain and here we have the density with the accumulative strain here the temperature with the velocities arrows and here the viscosity and so and finally here if i here i generate a gif with all these pictures with more pictures to see this the evolution but this is the gif that i show you to to the beginning of the of the notebook and there you can see the evolution of the the of this model so you can see how to generate the files in the upper and lower crust and how to this break up there the little tree so it is all for now thank you for your attention i hope that this tutorial will be will be useful for you and if you have any doubt about how to use mandiog or run a simulation or something you can find us in the slack channel or in their mandoc repository so thank you wonderful uh thank you very much um i can see that it was a fantastic and i loved that you started you started at the perfect level for me anyway i hope others found the same thing um i could follow along i love these animations really really cool makes me wish i had a cluster or something to to run those things on but i might give might give one a go um so anyway thank you to you and your your team i know there are a few people in the slack there and others that have been involved in writing all this code and thank you very much for open sourcing it makes a massive difference to the rest of us and to our community um see you all uh there's another tutorial in a little while later on today and then several more tomorrow so hopefully if you're watching the live stuff we'll see you soon

Channel: Software Underground

Views: 616

Keywords: geoscience, python, conference, simulation, finite element, FEM, mantle

Id: wzrIF4zpshM

Length: 72min 52sec (4372 seconds)

Published: Mon Apr 25 2022