 like Brian at the spot and we're awake now that's good we have to move on I think we'll have a chance to revisit the same questions a little bit later we're gonna move on to the numerical part of the problem aspects of the problem the next speaker is going to be Louie Moreci of Melbourne University telling us about the quagmire algorithm parallel landscape evolution models or coupling to large-scale tectonics thank you so yeah I'm I'm going to get I'll take the the obvious statement that you're a geodynamicist probably think I moved out of geodynamics a while ago so I'm not quite sure where I fit in any of the communities but I have followed I have been following a lot of the discussion about trying to make trying to make large-scale tectonic models interact with some form of surface process algorithms so what's the the inspiration for the sorts of things that I became interested in or think about a questions on a very very large scale I think and those those questions relate to the long tectonic evolution of the whole plate-scale processes and thinking about what the surface is telling us about the evolution of the and it's telling us a couple of things obviously most of most of what we heard about is is topography that that raises and lowers as its primary direction if you like and that we can describe everything in terms of the latitude and the vertical direction but it's also worth you know if you look at it I made this map say well if you look at look at what's happening in some typical piece of topography around the world you can say well there are places where everything more or less an equilibrium in the sense of things flowing downhill and very broad sense and there are there are other places where things take very fortuitous to get where they're going and obviously before some thing about the history I'm gonna take tonic modeling point of view these these sort of continental conditions have been the bread and butter of many many modeling groups over the years and trying to understand how to interact with the continental the continental with the continental strain rates and so we can build quite sophisticated models in which we understand how the how these how these collisional processes happen what happens about what happens from the underlying process the question what more can we bring to that if we do some better understanding of what the surface flows are doing and the information carrying to us about this process so here is here's a movie of sort of a typical collisional experiment that we would do in the conic modeling world of some the production zone evolving plates moving plates mostly here laterally and and as they and in this case we have a an indentor going into the continental block here which is creating enormous amount of lateral movement consequence there's there's vertical motion but it seems to me that one of the one of the things that really important to try to be able to do is to fill a couple with codes like that is not necessarily assumed that the motions are parameterized everything that way and that's the starting point that I took when I was a problem so in a way there are these sort of two different things that we're always trying to look at as we as we kind of approach these problems and I guess I introduce what I think of as the kinds of models that I hear a lot about actually which is trying to match patterns of in this case rifting models this is Roman Boucher's work and and he used the underworld code for doing this is very defamation and coupled that with the badlands comes out of the university group and the idea there is really trying to forward model stratigraphy in this case it's a it's a large-scale project that's interesting in understanding base and fill and how the hinterland contributes to that it's a very nice work on that but it's not obvious to me that the same techniques that are used to study those detailed processes will be the ones that we might look at if we were very interested in that original very large tectonic length and time scales and we might want to go back some more generic type of models very simple stripped back kind of models and just really understand what's going on but on the other hand we might want to consider we're starting again in terms of how we think about the problem what's the best way to build algorithms that can sit nicely with the kinds of codes that I'm talking about in the large scale electronics which have very large deformations in both the horizontal direction and vertical very vertical movements and tectonic erosion of topography by sort of folding and buckling and stretching and other processes so from the point of view of how you build algorithms like that the idea was to try to come up with some way to look at to build a sort of a flexible approach doesn't really depend on running with one particular code or in particular running with one particular set of equations or one particular problem in mind so of course I'm so thinking about that and it's taking forever to make much progress so it is interesting to look at to think about some of the aspects of this problem and some of the reasons that the problem can be a bit hard to solve in general and so you know one of them of course come up many times is that there is some form let's be abstract about it some form of downstream aggregation across a whole large natural extent of the system there is this accumulation of perhaps it perhaps it would be water flow in livers but it could be ice agree any of these things where there is an accumulation and you don't you know the more and more accumulation you have the more potential you have to greater time so in this case from this little kids experiment you can do you can play with this thing in Houston and if you connect sufficient numbers of the pathways together you can drive a little water wheel here some exciting fast rate that causes everyone to get splashed so that's kind of the principle of the operation it doesn't really need to be much more complicated to develop the computational side of that problem you know it's very well behaved and then you can you can look at these that you can develop that pathway as a tree as a sort of as a graph which has some very nice properties that's not cyclic everything is aggregating past parts will be very much independent of one another that's really that's really quite a nice property and it's exploited extensively in making efficient algorithms and we can use this kind of all this characteristic is the is the party sort of the area term that appears in most of these stream power formulations and the matter how you vary those points important to be able to look at some of these upstream type processes where aggregation is occurring from the point of view of from the point of view of the other side of this tectonic problem at a large scale what I would just mention you know Australia sort of and tectonics aren't really the sort of things you think about together at least not a present day but but that's not really the case right it's just very slow and very subtle and so in a very flat landscape and one in which there's very little but anything going on no we're not really not very much flow and not very much erosion then it's might maybe possible to record a lot more subtle signals right then you can record or then then thinking about the very very large-scale variations in the in the tectonic environments that have any collisions on so just a few meters shift in topography on the very very longest wavelength also the system and the fact that nothing get nothing you know that that is then recorded for a very long period of time so there's actually a possibility here that you can start to see through the see through this system into the deep lithosphere and the potentially the instabilities that occur below that are very subtle and you don't normally see so there's a sort of again it's a very long wavelength very large-scale process that we're interested in trying to understand but it's almost the exact opposite of the monogamous and there are quite some modeling challenges in this example including the fact that doesn't look anything like the graph that I just drew so many of the river systems in Australia don't really go anywhere they sort of go into dust right I mean look this one's backwards that comes out of here and then this and actually it turns out that so the characteristics of the problem I don't think I really need to slide the point of view of planning anything but just really the point out that if you think about the way that these there are there are two sort of sets of processes one one a local one having to go to do this aggregation they really need to be treated differently and in and of course they always are treated but very much from the sort of development of a parallel algorithm that couple with other codes be nice to be able to treat those in a way that is sort of not particularly special you don't have to develop some something that is not compatible with other kinds of codes and in particular one of the things that many of the tectonics codes do they're just very much parallel and distributed distributed across many many cores of large-scale computers and under those circumstances you know we normally have a pretty standard approach to solving that problem make the equations render everything in a bunch of things find suitable solvers for each chunk of the problem put it all together in a way that and that process of going through the decomposition is sort of the work of that is put on the library and so something we sort of forget about that but if you look at this sort of pathological example of taking a model in which information being aggregated then tells you some of the hard parts of doing right which is that if you start so I just point out the flow of information here and I numbered these boxes representing how you distributed this domain across a bunch of processes somehow or other I think I shifted that a bit so the ones that get numbers but anyway the information sort of starts out up here moves around we're needing to track this information as we go around and as you can see this is this is cyclic in the sense that the information passing from one place and algorithms which are actually explicitly worrying about the passage information place to place actually have to worry about you know have I finished this one can I carry on so there is this dependency which is non-trivial and although as I pointed out this is a contrived example it's also not completely the problem itself is not available you actually build an algorithm that is based on detecting whether or not you completed something for moving on to the next process and then that cyclic dependency becomes much more hard work so I so there's this bug me for quite some time and and then I decided to have a bit of a think about it and I realized that the answer was to go back to the simple things which is to think about stuff in terms of let's see if we could do the same thing render all of our equations in matrix form and then not have to worry anymore you know give all the work over to the linear algebra linear algebra packages design the dealing with and it turns out of course that there is a way to do that everyone already knew except me apparently but that you just render everything that in any graph you can draw your matrix that tells you how you passed along the graph so these simple directed graphs which have very simple structure to them actually lend themselves naturally to the process and you can you can look and see how it works that if you have a whole bunch of information abcdef nodes and you have some arrows which you reproduce the passageways and the toy problem then a simple matrix multiplication with let's say this matrix times the original nodes produces this distribution of information on the notes in other words we've taken a and b and put them together at C and we've taken e and move it to D and that's one iteration if you like of the information across the circle you can see that you keep cycling through this you will propagate information on wherever it is on the domain out towards so we can think of that as a simple as a matrix there is a matrix you can think of a matrix B as something which is an operator to move information one node in a downhill one node across graph and where D then it moves to find and D itself is incredibly fast and the total the total operation that we're trying to do is all of them one again and again imagine we'd like to accumulate so each of them that's not a sparse but in times but if we were to do that then well parallelism is quite straightforward because we could build a matrix and we use parallel matrix library to the actual implementation you would probably not do by building that matrix doing a lot of matrix multiplication because there's a recursion that exists if you keep multiplying and and summing this problem over and over again then you will achieve the results if you do it any times you will achieve of building the matrix and using it once approach so they are equivalent the two representations are formally the same and they look the same domain decomposition is in principle pretty straightforward and you can you can sort of convince yourself that as long as you have some overlap of the two domains in which the structure of the this case triangulation so the structure of the the neighbor relationships in the points would be reproduced then you can guarantee that you can actually build all of these matrices piece by piece without having to construct the whole whole so in principle this is not a this problem and here's an example of doing something like that I can show you how to do this one but essentially this is a sort of I mean there's some bits and pieces here and there but I took this I took I took the top one the photography model I did that because it's relative to geoid so at least as a starting point it should be straightforward and just took a million or ten million points or something just a triangulation and then these domains here are the whatever whatever just let Pepsi in fact do some decomposition and then just showed how the flow paths can be calculated and so these these here are just pathways not actual which is the flow pathways if you like across the landscape and they flow naturally across the boundary there's no special boundaries are not special so we can do that there are some actual there's some nice things so for example if we want to go backwards if you want to start up the alpha points and propagate information upstream then we just use the transposes that original matrix and that's sort of one-to-one mapping anymore that's what you want right you want to take each one every neighbor that it has and so that's a trivial operation you can find the transpose of the matrix immediately and you can multiply it out using exactly the same routines that I just talked about you might not want to do a summation but you might well will do transfer information you know that's just that's I'll show you what you can do with that the second thing you can another thing you can do very easily is that multiple pathways for flow a trivial to accommodate because we can build sort of a downhill flow matrix that goes one neighboring node or two just and then you can do one to a second neighboring node and we just have them together with some weighting function which is a diagonal diagonal matrix and w here can be whatever you wanted to be slow to be based on the delta function if you want to do things differently however you want to do it because because I've learned very much that whoever is in the audience today will probably want to do it differently because next time so people who worry about sorry people who worry about matrix parallel worry about doing this efficiently will worry that D1 and D2 are quite fast so I should I should point out these ones of these ones are both sparse of course together they're a bit less sparse but this matrix here is now equivalent to the one I would use if you want to go back up hill using these multiple paths it might make sense if you want to do it then you just do the transverse of this one and again you can do you have this concept of multiple actually by the time take these ones with multiple pathways and take the full matrix then it's in it becomes pretty dense and the purpose of doing well one of the purposes of doing that there are there are two purposes I can think of one of them is sort of actually physically represent the fact that multiple pathways will be very useful and low relief setting the other reason it might be interesting is to look at this which is which is this is what happens if you build a matrix out of the first downhill neighbor direction on a perfectly smooth of this part of the surface that is very defined by the cylindrical function and then these these points here and the interesting thing is to look and see how much of this totally uniform distribution ends up uniform and these are the catchments that you might form out of that work and of course these are not uniform and that's entirely to do with mesh discretization if you introduce making two neighbors or three neighbors then this kind of eliminates some of the effect of the mesh is I think another thing that you obviously need to do is if you just take a dm and randomly do this process to that is that there's a discretization error associated with sampling that if you don't sample the center of the then you'll find high points you might sample on the side in one place base somewhere else so there's a the low parts of the surface very artificially interrupted particularly in this case this is a model that's just a couple of kilometers in scale and so most of it is but but that means that we need some approach to saying well there are in places where actually the original original surface is not really where the flow paths go so here I took a landscape and we shot gun to it blasts a load of holes into it and then ran an algorithm over this in which we use that hatch with extra flood-filling process of reverse reversing the flow pathways to detect where the still ways are and then just fill them in well and that's what we get these holes are if you rain on that here's a slightly simpler example to follow but it shows that it works in parallel because here's a sort of egg box that's tilted and and then all of these sort of dark patches and light patches of the parallel domains and then we apply that filling algorithm potentially agnostic about the parallelism and we contain something that looks like that waffle waffle filling algorithm we're going to apply that to that landscape I just showed you which is what Macquarie if anyone's interested and and this is the sort of small changes in height that are needed to correct the landscape I won't say corrected but develop where the pathways should be so in other words if if we take the raw landscape and start filling it up to find out what which pieces which is not which is not to say that there's something doing to it just talk about the one do some analysis of the pathways that then redeveloped and you can look on Google Earth and see which ones are actually real and which bits of this something out there something else in the problem doesn't have but it doesn't show up correctly in the DM and we can do the same on the very large-scale and well they're not we're sorry they're not reverse these are the flow paths in Australia most of these don't connect up this is dry barren spots which don't get any rain to make these sort of connections occur and in fact over here of course really does get stuck so we implemented all of it I mean you can try this so we called it quagmire no particular reason but it's Python based but it's based also on Etsy and Pepsi for pie so it makes it reasonably fast it's not optimized for fast but parallel fast and parallelism more or less trivial the Etsy aficionados that we use we use the appropriate data managing structures within Etsy which gives us a whole bunch of interesting stuff we also had to go back and build a bunch of matching classes for this well they don't play with those and as you can see it's implemented mostly in terms of operating for those of you like working in Jupyter notebooks that's probably way easy so you can run it in the notebooks parallelism that's probably not really enough but you can just export a Python script maybe take out some of the interactive parts and then you can run it with MPI run the interactivity though it's quite nice because a lot of these kinds of algorithms they do I mean you they're not supposed to be bulletproof you just play with them and do what you want to do and we want to have a look and see what you're doing so the interactivity is not a useful scripting is also there for doing very large-scale things so the philosophy is very similar to that for underworld but we converted underworld to be very similar to this just a very big bag of relatively efficiently implemented Python classes that you can throw together to do whatever you want but within the limitations of the things that we provide you do it it's an example to make that work and the idea is that the science questions then dictate what the codes do rather than the other way around which is what we've been trying to do in the tectonic side much more than you know much more than the other way around you don't really want to be doing certain problems just because that's what code we make the code very very flexible the disadvantage of that is that so anybody in the community then can define their vision their vision their version of the problem and and this code would support that choice of course it's open-source hopefully an open community as well we leverage as much as possible on things like Petsy and produces strategy then for you to interoperate with other Python codes and installation would be through PIP or SPAC or one of those kinds of package managers I haven't uploaded the code for you to do that yeah but you can see they're fairly soon it's not shrimp rat rat to do anything other than show a few examples in that case it's kind of a typical Python project is a bunch of stuff you can do maybe it may be just people take it and do some of the algorithms I don't know if you want to try then the most robust way is to go to use the Docker for a limited time only you can go to this you can go to this URL and run it and if you do that so I was running it on my phone this morning so that's running in the cloud and you can log in through this port there's only one machine behind that so I guess if everyone does that I'm actually last very long but you can more than welcome to try and I'll restart it if it's a github project that's that's where we're at you can go if you go into this URL you'll see sort of a front page like this really just kind of a nice landing point after which makes a lot of learning time for questions both course University of Mines thank you very much for a very cool talk and I just want to reiterate for the people in the audience that don't get catch that the major problem for doing lithospheric scale models essentially we must use MPI parallel codes there's no way to do 3d without this and this has been the main hampering of coupling this kind of code because almost all landscape evolution models have been running either on a shared memory machine or on one processor so I think what you are showing here is an absolutely crucial stab in making progress so I do have a few more remarks regarding the you were showing that you implemented the matrix form and also the graph form what is more efficient for which not well optimized the one optimization you do actually sort of pre-multiply a bunch of combined but basically I would get my guess and the second question would be you were showing a little bit about the match sensitivity and we haven't heard very much about that for the last few days but my sensitivity is kind of a big issue and the landscape evolution models is in my experience but does that get rid of that by connecting it two or three neighbors I guess what you see is one gap and the other thing that you can look at the other mesh sensitivity is sort of an obvious one is that if you're doing something like this where you do the pre-processing or the other part of the mesh sensitivity I can think of is if you mesh once mesh twice and there's some randomness in it then you might get a particular example poorly poorly chosen example because that's a problem and there's multiple channels that we connect so depending on how you try to try to mesh this sometimes one way or sometimes the other way so all of these things could be parameterized as possible channel and sense that one mesh they exist in a different connection so there's a mesh sensitivity to the meshing so it would not be about multiple paths represented with the paths of the same mesh but it's multiple paths of multiple meshing instances so I think you could you could certainly construct different you could certainly would give them problems let me say if this could be plagued by the difficulty around somewhere around thanks for the very cool talk I um interoperability with Underworld do you mean you can truly do um two-way interaction or that quagmire just whizzing in the Underworld yeah so what we had in mind uh so for this at the moment the only impediment I guess to simply firing the two things up within the same Etsy environment is that the Underworld Etsy environment's been out of date dates back a few years when this was not quite so straightforward but you should in principle be able to use the same collection of processes and contribute across so that your data sharing is done across just the surface processes of the processes that's computing the surface of Underworld right but on the other hand it's I mean I don't know if that's that equally may not be the most efficient on the same machine and you talk through some other mechanism I'm not really I'm not really sure about that part of the problem hi uh Ralph Alto um that was a very interesting talk thanks very much um quick question you provided some great examples of uh Australia uh Australian rivers don't necessarily follow steepest descent flow and you made some comments about the model stability using steepest descent I was a little unclear on whether or not this code is steepest descent or multiple multiple directions it's both so that was the point you can you can construct actually I mean it's by default it uh constructs constructs the first two descent directions and your choice of which one is your choice which one to use the reason I have two directions is for the filling algorithm that is much more straightforward if you have two neighbors to know when you're a saddle point but to find to find when you're straddling to when you're when you're straddling a sort of a ridge the easiest way is to find out if you're two two or two of your directions fall into different categories and that's sort of how the process will be worked but it constructs both it's your choice of whether you use two or one or two more and use that it's just briefly coming back to Boris's question partially given that the the catchment areas will change you're going to have to reassemble this d matrix all the time and can you exploit parallelism for that in a very straightforward way so the catchments are not in not d doesn't contain explicitly what that information is so if you want to know what the cash information is you go backwards from the output points it's a formal construction to calculate what they are so the d matrix you do you build it on the fly so you do need to build it each time you remesh because it's basically a connectivity make it tells you what the connectivity mesh is sort of the connectivity across the topography so topography changes you should change it so you are always rebuilding it but it's totally local as long as you've got shadow information that goes as far as your stencil and of course you would you would do that right so so the shadow information stencil is is the nearest neighbor's bus a few extras to give you a bit of extra choice to find things there's an alumna Adina Pusuk, SIO, UCSD actually something related to Torsten's question so you assembly the matrices on the fly so dynamic assembly how does that affect the memory load? we know those matrices have a maximum non-zero structure and so for any given one that we make we know we know a priori what the memory allocation is for a given distribution of nodes we're given decomposition of the nodes across the system we know we know how many non-zeros effectively there are in the matrix maximum and there'll be a few of those will actually be zero after all because you'll be at the boundary or something so in terms of memory allocation that's already known and so you just use this filling to make sure actually so you pre-allocate all the memory it's not really efficient. Hi Alison Anders I had a question about the pit filling and external drainage so you force you have pit filling but you don't end up with Australia becoming all externally drained can you talk about how that works there's there's one yeah so the one thing is that it finds lake air a lot of things work their way down to lake air which is below sea level so it just decides that that's one of the places to go and then I guess just like the real Murray River it struggles a bit to find the way out without filling it and it'll fill up so I mean I just terminated that I ran it a bunch of times figured out that there was a low point in the system around just in land of Adelaide I guess and you couldn't buy it you didn't find a way out without continually cycling again but I actually thought it was kind of interesting I didn't yeah I didn't completely fill right out to the top so I mean I just I just put in the termination. Great talk Laurie and I know you and others have been working on this for many many years and I'm really glad to see that it's coming to this. Question you you use this tree structure which others you said have been using also to build the flow path and the accumulation of flow but as you may know we have also used that to solve the equations and it helps because you can actually in some cases it leads to an implicit formulation of the integration of the tree. Can you do this here as well? So I mean the thing is I don't build the tree I built a tree for illustration purposes to make sure they were the same and that was the last time I did it right so I did I sort of demonstrated for equivalence and I think that was my response to Boris's question that was that we I'm guessing that once you formulate things in the matrix form there'll be a way to look at it and go ah I see how to do this and the the abstraction that you so that you could work with the you could work with the full matrix representation and look to see if there's some way to speed that up um but but I guess the we then built it as a virtual matrix which is that we use the recursion to do an iteration that shouldn't necessarily be a problem if you have the right structure that's what we've been doing in Stokes flow for years right we never build the whole the whole thing that we want we just leave it as a virtual virtual matrix and solve that um so I I've thought about that there's a lot of I do think I figured I would just release this let's end on that note and thank both of our speakers again I've got one more thing but there's one more thing I should say no one announced the obvious and um that's good because that means that no one was like sitting there trying to do that but if you do want to have a password is back my up all right thanks a lot Louis for making this accessible and just a reminder we have the picture coming up there's the coffee break and then yes and please if you're on the round table please see Louise who's in the back