Alright, well let's get started. So my name's Mark Oaks. I'm the head of the technical committees for the Islas Americas conference and organization, and part of that responsibility is organizing the the webinar series. So this is our kickoff webinar for the, for this year.
And we're gonna have hopefully a webinar of just about every month leading up to the Islas Americas conference, which will happen in May. Our next, our next webinar will occur in October. Alessandro month tonarrow, will give us a talk and you should be able to see the kind of slide.
That's going by and there's a QR code there, so if you want to register for that, that webinar, you can click on that or scan that with your phone. We'll also send out kind of a an email as well. So if you don't get it today, we'll advertise it the same way as we advertise this webinar. And then just kind of a plug for updating conference upcoming conferences, so Atlas Americas will be happening in.
San Francisco in May, so put that on your calendars and then next year we'll have I class, which will happen kind of later in the summer, and that will be in Miami, so the two two great conferences coming up that you can look forward to. But today we have Babyan, he's gonna give our our 1st webinar talk for the to kick us, get us kicked off. And he.
We're in this position because he were on the Martial award for, for his paper. So the Martial award is for the kind of best overall conference paper so it's quite an honor and we wanna share kind of some of that research a little broader through the webinar series.
So I was kind of reflecting on back and how I met Fabian and we, I, my PhD advisor was Olivia dejar Dawn and Fabian did a postalk with, with Olivia. And so I kind of got to know him through, through that connection and.
Early on I felt like Fabian kind of just picked up the numerical methods for multi phase flows really quickly and and then kind of the the numerical methods for multi phase flows was kind of stuck with planes and Fabian said well that's not enough, we're gonna go on to this more complicated thing which everyone up into that point thought was just kind of crazy.
But he's kind of taken that on and kind of pushed the limits of what numerical methods can do, and I think that's really why he got the Martial award. He he's kind of worked on something that many thought was kind of impossible and and made a ton of progress with it. So, he's gonna share that story I think. So I look forward to hearing that story again and and kind of learning all the little details. So with that I will stop sharing and turn this over to Fabian.
He's gonna present for 40 ish minutes, then we'll have time afterwards to ask questions. And at that time I encourage you to to unmute yourself, turn on your camera, and we can have a little conversation about things. So alright, I will stop talking and turn this over to you Fabian.
Alright, thank you. Thank you very much Mike for the, for the kind introduction. I'm sharing my screen now. I'm hoping i'm hoping you can see it properly.
Yep, looks good.
Good perfect. Alright, well I guess I would like to maybe start by thanking Ailas and I guess the community that that decided to give us this price last year. It was, it was a very nice recognition, a very nice surprise. Mark already introduced the journey that led us here quite nicely and.
It's one where we try to do something that that hadn't been done and looked a little crazy and it to be honest to me sometimes it still looks a little crazy and I'll try to, to give you a, a nice overview of how we got there and where we are now. It's obviously not.
Quite there yet, but it it does show a lot of potential. So this is this is the title of the talk that was rewarded by, by ISLA conservative to facebooked simulations using peace waste quadrantic interface reconstruction.
And it's a paper that we worked together with and and Olivia Jada at Cornell, and it's, it's really the the result of of I would say a couple years of of research that, that led to this point.
I just want to start by acknowledging a few people that were instrumental to to all of the work I'm going to present today. So Rob Cody, who I met at Cornell, worked together at Los when he was at Los Anamas and, and has now gone private, working for luminary.
Cloud, Austin Hann, who also, whom I also met at Cornell isn't at Argon. And then Ashish Dana Kota and and who were both undergrads working with me here at the University of Illinois and and I'll show some of their work in in two dimensions. I think I saw Ashish in the list of participants, so hi Ashish, and they both went to Georgia Tech for grad school. So so they're not.
So they're not working with, with me anymore on this and also parentury is a is a grad student that joined my group last year and and and really is the one that that picked up all of this work and I'll show some very preliminary things that he's done and and I'm hoping that we'll be able to present some exciting new, new results at at Iaas next year. The funding wise so that there's there's multiple agencies that, that.
Supported this work. So it all started by some funding from the German research foundation. The bulk of the work was done through a, a Marie Schooldosca Curie fellowship that that basically allowed me to come to Cornell, and very recently the NSF has has given us some some funding to, to continue this work. So we're pretty excited to to get that going. Alright, so quick outline of, of what we'll talk about today.
We'll start with some very brief motivation and challenging. I don't think I need to motivate the work to to most of the attenders. We we kind of know why we're here. And then we'll do a quick primer on, on interface capturing and and and kind of focusing on, on the method of choice for today, which will be the geometric buff. And then we'll talk about.
What's the current state of of and and and how to extend it to to high order alright, so we'll start, we'll warm up in two D and then we'll move to three D and and I'll show you where, where we stand more or less now and and what I think is going to be needed in the in the coming few months and and what we plan on doing as well. So a motivation this is obviously a Island's webinar, so I think we're all aware that liquid sprays are are pretty.
Important topic and and having predictive capabilities for them is is also pretty crucial. So that ranges from, you know, sort of traditional fuel or or liquid sprays in incompressible compressible flows to ultrasonic optimization, so that's what you see here on the right, more traditional inject printing bottom left or or shock droplet interaction which we've seen quite a lot of at is this.
Here, all these to say that essentially the the tools we're trying to build, they aim to address all of these flows at once, really. Why is it challenging to address this topic where essentially you're dealing with a fluid system that displays discontinuous properties that sharply vary across your interface. Your interface is subject to a.
Surfaceful distribution due to surfacetion which can become singular, and that also leads to a whole family of, of, of problems. On top of that, so surface tension is dependent on surface coverture, which is a differential quantity that you need to extract from the knowledge of the interface and and we'll see that the state of the outset.
Is does its best at at approximating it but it's it's not perfect. And obviously these interfaces especially in the type of flows that were interested in, you know, in the islands community, they will typically go through multiple breakup or coalescence events so that the topology constantly changes and and that changes is also important to cap.
Sure accurately, obviously. So that's it for the motivation. Now let's let's talk about the numerical methods that that allow us to, to look at this type of flows, right? And so I, I was preparing this and and I I I decided to classify through three big families and I tried to give some strength and weaknesses of it.
So obviously this, you know, I'm pretty partial to one of them and and if you, if you feel strongly about what I'm saying, please feel free to to raise your hand at the end of the, of the talk. But the 1st big family is is the front tracking methods which rely on some explicit.
Reconstruction really of the interface through some fully connected triangulated surface or or let's say computational approximation of the surface. This triangulated surface is evolved with the flow of velocity at the forms. This deformation typically means that you're gonna have to remesh or or at least.
It's correct the interface constantly so a pretty high computational complexity, but the advantage is that this gives us a full discrete domain for the interface on which additional fields or property can live. So e.g., this is pretty well suited through the transport of the fact that.
But when it comes to using this for, for itomization, there's some some big issues that that means it's it's usually not used or or there are very few instances in literature that useful tracking for itomization. One of the, the, the, the 1st issues is mass conservation, so the approach.
In itself is is not mass conservative. You need to pay careful attention on on how you do this transport and how you, you, you rememsh the interface and even then you you typically have to apply some form of correction to it. This is a fairly complex numerical method mainly related to the constant.
Adaptation and remissing you need to apply to it. On the other hand though on the plus side, this is perfect if you're interested in really accurately capturing breakup orcharisons because you have complete control on the topology of the interface, so you can decide precisely when you want it to break up orcharis. And, you know, if you.
Side for that form to be refined way below your grid size, then you can very easily represent certain skill structures as well. But for automation again, not really a method of choice. Next to it we have well I lumped everything together in what I called level set/diffuse interface methods and that includes all the recently feature on face fields.
The broad idea is that your interface will be approximated as some ISO contour of some continuous function. So in the traditional level set sense, this would be a signed distance function in the conservative level set, this would be a hyperbolic tangent function of the signed distance function and so on.
And if I use my little red, green orange classification that I started with from tracking, I I put this as as yellow for mass conservation although I I do acknowledge that, that, you know, recent version of this have some much stronger mass conservation properties which are typically related to some.
CFL constraints. The big advantage is that because of the smooth nature of this interface representation, you can pretty much use your traditional finder different schemes like off the shelf and extend pretty easily to high order and so on. So that's definitely a strength and and and you know if if you.
You're interested in running this on on sort of massively parallel compu computing resources. Well, you can make it really scalable fairly easily. On the minus side, well, because of this representation, you pretty much have zero control on numerical breakup and or or any ability to do sub scale representation and as as we'll see a bit later on, this is something that's that's really kind of growing in the field of, of.
Automation and and you know the the this is an area in which these these diffuse interface methods tend to like. And so finally we bring we're coming to the the last family, which is that of of volume of fluid or or moment of fluid methods.
And you could technically classify them as a subset of level set methods in the sense that the ISO control you are interested in now is the ISO controller of some heavy side function that identifies your fluid phases, right? So you go from zero in in one of the phases to one in the other phase.
And your interface is essentially at the transition between this, this discontinuous profile. Now one of the key advantages of the volume of method is that you can design it in a way that is probably conservative in a sense that, you know, you, you're gonna conserve mass or volume down to machine precision.
And that is regardless of the CFL that you're using for your simulation. This comes at a certain cost, so the cost is mainly in that you need to apply some pretty complex.
Computational geometry algorithms to conduct all the different intersection that that need to be done for the direction and and we'll talk about that at length later. So I won't say much more about this now. One of the like sort of I would say more recent advantages that has come across with, with.
This method is that. It does enable some degree of service care representation and and as a result, some degree of control on breakup and which is which has great potential for automization because it means ok with re relatively close resolution, we might be able to keep track of structures actually below the great size and.
And retain some, some pretty predictive capabilities for, for this type of flows. So as you can imagine that this talk is gonna be focused on, on the latter, the the volume of method, and and we'll talk about data a bit more now. So.
I, this is this is I, I took this as a webinar so so I I I, I decided to kind of go through the underlying theory and mathematical foundation that, that, you know, led to the design of those schemes. And so we're starting with with the the indicator function itself that, that.
Defines the volume of methods so the way you identify your your phases is purely with this binary indicator zero in one phase one in the other phase. Mass conservation tells you that this indicator function which is discontinuous is is subject to a very simple transport equation. It's material derivative is equal to zero. And so this allows you to discretize this equation and solve it with the relevant scheme, right? And so in a final volume sense, for instance, you would integrate.
Great in time, the integral over a computational cell of this transport equation, the transent term gives you simply that the volume of the cell is multiplied by the difference between the volume fraction in the cell at the next time level minus the volume fraction in the cell at the current time level. And then the right hand.
The 2nd term gives you this big flux integral, so integral in time of of the integral over the surface of your cell of the flux of the indicator function pretty much. And this is where essentially all of the challenge lies when it comes to the volume of how to estimate this integral accurately and in a way that guarantee is gonna.
Limit numerical dissipation and and you know preserve some, some accuracy. And so the state of this, this is none of this is, is new, right? And the state of the art in approximating this interval has been for several decades, some piece wise approximation of the interface and typically this approximation is piecewise linear, right?
And so from the knowledge of this interface position with, sort of like purely geometric routines, we can approximate what should be the flux of indicators to function through this phase. What are the methods that are used for this interface reconstruction and flux calculation? Well, the interface reconstruction.
Back in the days I would say was, was done through some modified finite difference operator, so whether it was used directly on the Voian fraction or on some regularized version of the boeing fraction fields.
High functions have also had a very important role in the estimation of these interface differential quantities. And so to these dates, you know, the the open source code that is developed by Stefan Popineer and previously Gary's have heavily relied on high functions and and.
Been very successful at using them for, for the volume of methods. But there is also another family of methods which rely on some minimization of a cross function for what would be the optimal piecewise linear reconstruction of the cell, based on some metric that is completely user defined, right? And, and the most famous of these methods is the method by P.
Which essentially makes you minimize the volume fraction error in a local stencil of cells under the constraint that the volume fraction has to be conserved in the central cell of interest, so essentially mass has to be conserved. This becomes a nonlinear constraint minimization problem, which you have NO.
But to solve iteratively essentially. So it comes with some costs but it can be shown to be second order accurate and also it is relatively straightforward to extend to non Cartesian measures. When it comes to the infection itself, well historically things were pretty much started using operator splitting, which is fairly standard in the sense that instead of.
We're Solving for that big flux integral all at once, we split the problem directionally and that relies that hinges on the fact that the mesh isian, so we can do one sweep in in in the X direction, e.g., followed by one sweep in the in the Y direction and so on.
By doing so, the flux volumes that you need to intersect with your interface reconstruction are rectangular into the your rectangular kubroids in, in three D, which allows us to use analytical expressions for the calculation of the whole infraction. And so that's, that's a huge.
Huge plus which makes it pretty efficient. When it comes to the conservation properties of the scheme, well in two D you can get away with some combination of implicit, explicit formulation of the scheme and and be probably mass conservative. In three D that's a different story and essentially you need to introduce.
The dilutation term that, that will ensure that your volume fractions remain bounded and and this is the work of of Waymooth and you from 2010. Without it, you, you really can't do three D infection with provable boundiness of your voying fraction using operator splitting. So conservation.
It's possible, it's not it's not necessary straightforward or or requires the introduction of this additional term and the extension to momentum and energy conservation is also doable but, but not straightforward either. On the plus side complexity scalability is is definitely a green because we can make use of these analytical expressions for computing the the fluxes but you can't extend this to instructed message, right? Because this, this, the method itself is built on the fact that you.
Can do this directional splitting, so your, your, your mesh is aligned with these directions. But it is, it is a method that is used quite widely in the community and and it has been used for simulating automation. So what you see here at the bottom is the simulation of a of a pulse liquid jet that was published recently in in GFM which usesk.
As opposed to operator split invection, obviously there's another familiar method which we call unsplit methods or like schemes, which essentially rely on some approximation of the.
Transport of your cell in a lagrangent manner. And I'll I'll I'll maybe formalize that this a bit better in a in a couple of slides but you essentially solve for the invection all at once now, which means that your flex volumes become much more complex.
They essentially take an arbitrary shape, which is you know dependent on the topology of your, of your cellfaces. But now you can guarantee boundiness of the volume fraction is as a result conservation by making sure that the volume of these fluxes matches the solinoidol fluxes.
Is that are essentially so the result of your flow solder, whether it's a pressure correction or or or else. And so as a result, the, the scheme is by construction probably bounded and it's not subject to an ECFR constraint for the conservation. So that's definitely a plus you can make it momentum conservative and as we've seen with the the the paper from Zongaozu and and Olivia.
At Highlas this year, you can also make this energy conservative which is definitely a very appealing feature. This comes at a cost of of geometric complexity, so you have to deal with arbitrary polygonal intersection now with a plane if you consider the piecewise linear standard and and as you will see in a few slides with a power bottoid if you wanna do this with with high order interface reconstructions. A big plus as well is that none of this relies on the nature of the grid. So I've I've illustrated it on the right with the cartisian grid, but you could apply the exact.
And there's really not much more additional complexity. And again, this has also been used quite extensively for automization simulations. You have an illustration year of a simulation running in in Olivia's group at at Cornell a couple of years back of a core actual itomizer. So let's talk about this method in a bit more detail and and what is.
What is, what are the maths essentially behind it? So the idea is that if you have a given computational cell and and you want to figure out what is the volume fraction in that cell at the next time level of your flow server, where you're gonna assume that this cell is essentially identical to some like rangeant volume that is moved by the flow at the time TK plus one, right? The time at which you want to predict the volume fractions.
This like rangent volume is is following the floor, right? So at the time TK for the previous time level it would be located in this position. And if you apply the renals transport theory into this volume, you end up with this relation which tells you that the rate of change of your.
Indicate the function integrated inside this volume is equal to zero in the case of incompressible flow or is equal to a term that is proportional to the velocity divergence in the case of compressible flow, right? So if we stick to incompressible flow essentially what this tells you is that the volume of liquid that you have in that.
Like rangent volume, and by the way, the instance that Tk at at time TK we typically call the pre image is conserved, so we know exactly how much volume should be in the in the cell at the next time interval. We don't know where in that cell the volume is located, which means that we we are going to need an additional interface reconstruction step, but at least we know how much is there.
Now, numerically, once you try to do this, you're you're going to introduce a number of different errors, right? The 1st one comes from the, the, the transport itself or the corners of your cell along the path lines of the flow. You're typically gonna do this with some numerical integration scheme, and so you're gonna introduce a.
The numerical error resulting as a result of that. 2nd error comes from the fact that as you transport this cell backward in time, well, it's its faces or edges get get curved by the flow typically and and if you don't want to have to deal with arbitrary curves as the edges of your cell, you will approximate them and and the standard is typically to do this with a piecewise linear approximation. And we call this the representation error, which relates to the the representation of the flex volumes or or the pre image essentially.
And finally, we do not have exact knowledge of the interface at time TK, but instead we have some approximate knowledge of the reconstruction that comes from its reconstruction. And this reconstruction can be piecewise linear and we'll show examples later on where we try to extend this to higher order.
If you add up all these errors, well, you make an error obviously in the estimation of the volume fraction. Although you can still guarantee bandoness if you construct your flux orium appropriately. And if you test this on a very simple infection test case, well, you realize that this error of the transport of your volume fraction is essentially second order.
Right and it relates to the fact that everything that every approximation that we've made is is piecewise linear essentially. So this is a case where I transport this very cool interface in a periodic mesh which is unstructed, so it's a it's a volunoid mesh essentially and the the interface reconstruction is piecewise linear, and if you.
Should the error in the volume fraction field after one flow through the box, then you find that your shape error which is essentially just aL1 norm of the volume fraction error converges with second order accuracy. The, the symbols here, so the squares are for and the pentangles are for the instructed mesh I show on the left. As a result of 2nd order accuracy of the volume fraction transport, the curvature that you will estimate from the.
Reconstructed interface will only be zero zero order accurate, and that comes from the fact that curvature essentially requires a a differentiating your, your interface quantities twice, and so you essentially lose two orders of accuracy, so the the curvature error stagnates, meaning that you just simply.
We can't throw infinite computational resources at this problem to solve it because the the curveature error simply won't won't go further down. So that's a that's, you know, clearly a a limitation of the scheme. It doesn't prevent it from being used successfully, but you know oh sorry but as as a as a numerical modeler and and you know as someone with a apply mathematics background, obviously.
Not great to have a method that doesn't converge with mesh refinement. So one of the things we started looking at a few years ago is, ok, what do we need to do to improve this and and essentially what what that means is what do we need to do to increase the order of accuracy of this scheme so that we can get a convergence estimation of the coverage. And so that's what we're gonna talk about now and and we'll start by by looking at the problem simply in two D.
So ok if we want to increase the order of accuracy of the method, we need to reduce the leading order of all these different error terms that I've that i've lead, that I've listed earlier. The 1st term which relates to the transport of your cell corners backward in time is is relatively.
Easy to handle because you just need if you use high order interpolation, high order time integration, you can pretty much reduce the leading order of that term fairly easily. The main challenge comes from the other two terms, so the the error resulting from the piecewise linear representation of the fluxes and the error coming from the piecewise linear represent reconstruction of the interface. And so the the what we propose to do is to well.
Switch from a piece waste linear representation or reconstruction of the interface to something that is piece waste quadratic so we could essentially bumps the leading order of all these terms by one. How do we do this? Well essentially if I go back to that sketch I had of my like original.
Production scheme, what I need to do now is that I need to replace everything that I had piecewise linear by something that is not piece was quadratic. And so that's what we do here in two D the edges of my pre image become quadratic arcs and the piecewise reconstruction of my interface becomes quadratic. In that.
In that sketch, what you see here as the points are the control points of these quadratic arcs, meaning the intersection of the tangents at the endpoints of these quadratic arcs. Ok, so let's let's talk about how we would do this starting with the transport of that premium is backward in time. So we need to be able to transport the edges of our computational cell backward in time. The way we do this is by transporting both the.
Points and the tangential directions of the edge at those endpoints. You might end up in a configuration like the one pictured here where your tangents directly intersect and that allows you to directly construct a quadratic art. Sometimes these tangents won't necessarily intersect, in which case what we do is that we split the original edge.
In half transport the tangents once again, hoping that they intersect again, if that's not the case, we keep on doing this recursively until we have something that that is satisfactory. You end up with this reconstruction of the edge, and as I said if you want to guarantee boundiness of the transport of your fraction, you need to make sure that the, this hashed area here exactly matches the discrete flux that your flow solver returns that is divergence three, and so we do this by.
Essentially moving the control points accordingly so that we can exactly match the the volume of the flux. So this is for the, the, the transport of my cell backward in time, so essentially the construction of the pre image. If you want to transport the volume fraction, what you need to do next is intersecting this pre image with the quadratic reconstruction of your interface, right? And so in two D the way we do this is that well each.
Arc is a is a quadratic arc, so essentially an arc of a power bottle, we parameterize it with a busier quadratic busier parameterization, same for the interface, and so finding the position of the intersection between my interface and my quadratic pre image essentially boils down to solving a qualtic polynomial, so polynomial of all the four, which is not that all challenging. Once you do this, you can parameterize the the region that encompasses.
Is this this hashed area here and you can just simply integrate its area using the divergence theory and and it just requires you to integrate this this guy which is also a polynomial. So nothing, nothing too fancy. You could do the same thing for higher moments if you wanted to, but.
We won't do it here. Finally interface reconstruction is also a key ingredient of the, of the approach. The way we do this for this simple to the example is that we just essentially take the the Elvira method from linear interfaces, and we just kind of like naturally extend it to.
To quadratic reconstruction, so we minimize the volume fraction error in my stencil under the constraint that volume fraction has to be conserved in the central cell. And on top of that we had this, we had a small penalty curve, so this guy is only here and only really becomes active when the curvature of your interface the radius of curve.
Goes below the grid size because we at this point in time we we we don't want to focus on on cyberscary construction, we want to avoid having anything that that goes below the grid size in terms of that, that raises a whole bunch of of additional problems. So on the right here you just have an illustration of of the results on the top is.
The sort of classical piece was linear virus on the bottom you see the the quadratic version of it and and the reconstruction of a circle from volume fractions. Once you put all of this together and you look at the rigid rotation of the circle in a2D domain, you can start looking at the the behavior of the scheme with different flavors of the of of quadratic representation, right? So.
What you see here on the right is the volume fraction error, which is which I already plotted for the three D case earlier. In gray is the source state of the art, so piecewise linear flux is piecewise linear interface reconstruction. In orange I use piecewise linear flux representation, but this was quadratic interface reconstruction, and in blue everything is piece was.
Quadratic. And so what you see here, this is a case where the CFL number is, is, is fairly high 2.5. Essentially the the choice of flux representation dictates the order of accuracy. So even if I have a quadratic reconstruction for the orange case, it essentially almost exactly overlaps with the the gray line.
Meaning that the flux representation error completely dominates the error. If I start reducing the CFL number now, well, I see that the the behavior slightly changes, right? And, and in particular, I can see here pretty clearly the point at which my flux representation error starts dominating the interface reconstruction error and my orange line goes from essentially what's very high order accuracy to some.
I should, yeah, I didn't mention it earlier, but the the blue line is always 4th order accurate, right? Which is kind of one of the messages of this. And if I keep reducing the CFL number, so I go all the way to .05, which is very, very low, but I I see here that both my blue and orange lines converge with 4th order accuracy meaning that I'm I'm in a case where the flux representation error is much smaller than the flux the interface reconstruction error and so I can still recover high order accuracy even though I have only the inf.
The face reconstruction that is quadratic and not the flux representation.
So this is a small two D kind of academic exercise that we did together with with Ashish and and Richek because we wanted to, you know, verify the assumption that we might be able to get away with not having a piece ratical instruction or.
Presentation in three D so that we can only purely focus on having a a quadratic interface for construction because that one that's already complex enough and two, just the thought of intersecting curve surfaces in three D and completing the volume of that intersection is, is, is.
Is really something that that that would tend to give me a nightmares so we wanted to keep it as tractable as possible. Okay, so moving on to three D what is what is what is the strategy? Well, we've done our sanity check in two D meaning that we verify that as long as we keep the CFL number low enough, we, we can get away with having.
Linear representation of the fluxes or pre image and only focus on having a piecewise quadratic reconstruction of the interface. And so that's what we're gonna try to do and so that requires introducing what we call PPX or piecewise power botic interface construction with the hypothesis that this should enable higher accuracy in in three D provided we keep the CFL low. Another byproduct of this is that this should enable the transport of high order geometrical moments.
And I'll tell more about that in, in a sec. I want to show this this little graph here because I'll talk a lot about what I call the forward problem and and backward problem in in the following slides. The 4th problem is the problem where I have knowledge of my powerbotic reconstruction and I wanna compute the geometric moment.
That rely from that result from the clipping of a polyhidrun with this public reconstruction. The backward problem is that in which I know the geometrical moments and I want to find what is the public reconstruction that.
Matches these geometrical moments in some optimal way, ok? And as you advect your whole infraction in time, you essentially keep going back and forth between these two problems, ok? You need to be able to solve both. So let's start with the with the forward problem and and we'll, we'll define it mathematically and then we'll go through it. So I was.
I was going through the slides last night and I was a little hesitant in, in going too much into the math and then at some point I was like, you know what, if there is, if there is one event in which I can probably, you know, show the community what's under the hood with all these, this geometry cutting operation, that's webinar, so let's let's buckle up and and and do it.
So the 4th problem solution is essentially the following from the knowledge of some peace waste quadratic interface approximation. I want to calculate what is the geometric moments of the volume that results from the clipping of a polyhedral cell or flux volume with this power bulloit, ok? And that's what I've illustrated on the right, what I call S in the is the is the power bluet surface. P is the.
And so I want to know the moments of which is the intersection of my polyhedron with a queue which is the the volume that resides below that powered surface. The implicit definition of my surface is given here. It's essentially ah function fee of x, which is equal to zero and and that function fee I give in its most minimal possible form, alpha x square plus beta y square plus Z is equal to zero. And the reason why we do this is we want to have as.
Few parameters as possible in the mix for calculating the, the, the moments and anybodroid that lives in the three D space can be expressed in the following form after some combination of translation and rotation, right? So we really cover all possible options by by considering this this limit case.
Okay, so, what is the quantity we're interested in calculating is the integral over my clip volume of a bunch of monomials, right? So the, the, the 1st one is is just one. So that's the integral of the volume. Then we have xyz, that would be the integral of the components of my barry center. And I could add more.
So if I'm interested in the moments of inertia, I could add X square, X Y x Z et cetera, and, and integrate 2nd moments. But for for the the the presentation here, we'll just keep the the zero and 1st moments only. The 1st step, if you want to compute the the moments of this guy is to convert this surface, this volume integral sorry into a surface interval using the divergence theory and that's what i've written here. There's multiple multiple ways in which you could formulate this, which is to formulate it such that essentially we're integrating along the Z direction for this.
And so this is the surface integral and 2nd line is the the splitting of this integral into a contribution of the curve surface, what I call this S tilled here, and the contribution of the F hat surfaces which are just the.
Flat faces that have been keep by the surface. Second step, I'm going to convert this surface integral into another surface integral, but this time on the projection of this boundary onto the XY plane.
And this is essentially a variable substitution, right? So I want to substitute the DA element I had in the previous interval into something that looks like DX times DY. And what's nice about this transformation is that actually the.
Integrats that result from this are just simple polynominals, so they're kind of trivial to integrate. What's not trivial is the Alex on which we need to integrate sorry the the the areas on which we, we need to integrate this and that's why there is some extra step. The next step is to convert this surface integrals into line integrals, so we.
Apply the divergence theorum another time. This is what the expression looks like at the end. It's become a bit more convoluted in terms of of mathematical notation, but essentially what it means is that we're integrating along the edges of these projected phases, ok? The last step.
I insert some parameterization of these edges in, into the mix, and so I can convert this into essentially a double sum, so some over the faces, some over the edges of my faces, of some integral for the parameter T which is the parameter of the parameterization of each one of those edges.
Good. So here we are, this is the the expression we need to calculate to be able to to solve this forward problem and and on the left here you see an example an example of the integration domains, right? So the integration domain is this oriented curve that you see here. It's a combination of trade.
Lines and and so the parameterization of these straight lines is trivial and a combination of konic section arcs and and that's slightly less obvious what to do there. The instinct would be to parameterize these connic sections using trigometric function because you know that's the sort of natural way to do it. The problem is that this could lead to singular parameterization, particularly in cases where your face is your face normal is off.
To the Z axis of your and and that's what I illustrate here with this curve is the same axis of my chronic section parameterization if I rotate my plane around the X axis, which means that if you were in the case where that angle theater is close to zero or pie, then.
Essentially the same access would would go to to infinity and so you would have to do arithmetic operation with with quantities that go to infinity. Your computer would essentially run out of digits and and give you a very poor estimate. What we do instead is we parameterize these connect section arcs using rational prenominals, so this is also known as a rational basic curve parameterization. This is unique and this is importantly well posed, ok? So there is NO none of these coefficients essentially go to infinity even even for those limit cases that are highlighted earlier.
Alright, we're almost done, so you could go ahead and and complete this interval, but you will still be played with some random errors and so you have to be careful in the way you do it. The way we solve the issue is that we split it into three contributions, so essentially this this integration domain that we have at the top here, we split into these three integration domains.
The reason for this again is is is error cancellation so that we don't have any singular configuration. Once you do this, you have a set of closetal expressions which are a non singular, well posed and essentially not really subject to to significant random errors.
But that's not the end because you can still have issues for very specific cases in which you end up with what I call ambiguous intersection topologies. These are cases in which for instance, your powerboroid directly intersects with one of the vertices of your polyhiton or your power boardwid is directly tangent to one of the favor.
This is of your polyhid. What we do to handle this is that we have sensors that detect those cases when the sensors are triggered, we switch to quadruple precision and we start shaking this this interface configuration essentially up until we find a suitable and and well posed.
Intersection topology. And that seems to do the job. So we've tested this extensively. Yeah, this is just a slide where I show all the formulas we have in the paper it's just to show you the the what they look like. It it's a bit bulky, but it's essentially mostly just a bunch of of polynominal expression. We have some hyperbolic.
Sorry some trigonometric function we need to to pay attention to, but otherwise fairly easy to implement and I would say most importantly very fast to calculate, right? Because there's just a bunch of addition and and multiplication really. Okay, how did we test this thing? We considered a bunch of different shapes that are representative to the type of of cells we would encounter in simulations. We have some cases where we had non complex.
This is because we wanted to test that and and this stands for Bunny has as a hundred thousand triangles, so we wanted to see if we were able to scale up. We did 200 million random test cases. We also generated specific cases where we have these ambiguous topologies to test the robustness of the scheme. These are examples of these ambiguous topologies and and we made.
Measured the moment calculation error and the integration time and and the conclusion is that essentially we're machine accurate, so the average error is on the order of the machine ten two -16. The maximum error for these 200 million cases is on the order of ten two -14, so pretty happy with that. It's robust, the algorithm doesn't fail on all.
All these these these test cases we've considered and it's relatively fast, and what I mean by that is it's it's only six times the cost of the piece was linear cutting, which is significant, but I think it's still within the range of of what we could afford for the, for these type of simulations. The code is open source, so feel free to go take a look and and and re.
Chat if you if you're interested in using it. Okay, so that's for the solution of the forward problem. The forward problem is essentially pretty much solved and it gives us a result that is machine accurate, unique, and robust. The backward problem is a different story and and and the the main difference between the backward and the forward problem is that.
The backward problem has essentially an infinite amount of solution in the sense that depending on the way you define what is an optimal public reconstruction based on the knowledge of the moments, well, you can have as many different type of optimal reconstructions, right? And so I'll list some of the different approaches you can consider, we've tested them and and.
I'll show some results. This is actually something that is, that has been done quite a lot in the literature already because people were using this to compute interface curveature for self astansion, right? So you could fit your power bottoid two points, that's a linear regression, pretty cheap, doesn't work too well. Second option you could try to minimize the volume difference between your click and your borderwide in a way that is formulated as a linearly square problem. So that's pretty fast.
And were actually very fast and actually works pretty well too. So very good method. You could also go brute force and just solve the nonlinear minimization for minimizing the volume fraction or the moment error in your neighborhood. That works pretty well, that converges with mesh refinement, but much more expensive.
And and we've compared actually these three approaches in the paper with with Austin hand and and Olivia that was published in IGMF last year. The main conclusion here is that the best bang for your back is actually the method so the 2nd option. But we'll see later on that it's it's, it's not necessarily the the solution.
That's that's gonna enable us to do a PP green. You could go further with this brute force approach, you could keep adding moments to your minimization, you could start adding constraints, so you could do sort of type I've already introduced that. You could also start adding higher moments, so you could add your 1st moment, your second moments, and so on if you want to. The idea being that the more data you have, hopefully the the more accurate.
The reconstruction should be. What I'll the results I'll show today is with move to reconstruction, so essentially we solve the nonlinear minimization in which we try to minimize the 1st moment and second moment error in the local neighborhoods.
Under the constraint that volume has to be conserved in the central cell. We again have some penalty terms that prevent the formation of of very highly curved power bottoids and and what you see here at the bottom is is some illustration of this reconstruction compared to click on the left, so click is on the left and this move to method is on the.
Right at different mesh resolution going from four cells across the diameter to 32 cells across the diameter. So I've mentioned higher moments, obviously if you want to be able to use higher moments for your reconstruction, you also need to transport them, right? And and so for the 1st moment, this would be the sort of traditional moment of read method of the and and in which you apply the transport theorum this time to X times the indicator function, it gives you an.
Make an equation for the rate of change of the the center in your cell, but you can also extend this to higher the moments, right? And so for instance, you can extend this to 2nd moments, the moments of inacture of your liquid if you, if you want. And this is the expression you end up with at the bottom for the rate of change of 2nd moments as a function of the sort of moments transported moments with the.
The truncated series expansion of your velocity field inside your.
Once you put all of this together, you can go back to your direction test case and, and start testing a few things, right? And, and I'm starting with a case where I'm using a piecewise linear reconstruction, and I get this sort of traditional result that my volume fraction error is second order. My 1st moment error is 1st order.
And my second moment error, which you don't really see much in the literature but which would use this this equation that I just showed in your previous slide is about the zero, right? So essentially, every time you consider a moment of one order higher, you lose one order of accuracy, which makes sense. If you now go the same test case, but you replace your in.
Interface reconstruction by this move to method, you essentially bump the order of accuracy of everything by, by one. So your volume fraction is 3rd order and all the way down to your second moment, which is 1st order. What you see here on the on the bottom left is the comparison between in orange the initial reconstruction and gray the final reconstruction after one revolution. So on the left it would click and on the right is with this reconstruction.
This is a slightly different test case which we designed so that the pre images are not piecewise linear because actually in the case of a rigid rotation, your, your pre image is piecewise linear so you so that that kind of.
Hides some of the potential errors you could make with your scheme, and the conclusion is is you know fairly similar to what we had in two D, which is that if your CFL number is low enough, you have this high order convergence behavior. If your CFL number increases, at some point you'll revert back to the sort of traditional secondary.
The accuracy of infraction and that's expected because that's where your flux representation error would become dominant. This is another well known test case of the community, so the the formation case, and this is the snapshots of what it or the video of what it looks like with this move two methods. On the right you see a comparison with the final results using peace was linear reconstruction. The key takeaway here.
There is that because we go, so this is a ranger resolution for which the interface thickness reaches the, is below the great size at some point, we essentially revert back to second order accuracy because our scheme in this form doesn't have service reconstruction capabilities, right? So we're essentially limited to second order accuracy.
Because we have this machine use breakup at some point. I wanna stay on this case a little bit because there's some interesting discussion to be done here. So, what's really noticeable in this case is that the the reconstruction method is able to quite nicely I would say break up the interface once it reaches a subgrid sized thickness without catastrophically failing at doing so and and you will see in a cycle in a in a, in a second where.
What I mean by this? And that's actually a pretty important feature of such schemes, right? Because when you look at atomizing gas liquid flows, you will have multiple, a lot of breakup instances and if your scheme is not capable of of breaking this up cleanly or if you don't have any cyber scale model that allow you to keep track of the interfaces below the grid size, this could really put.
With your interface your your simulation heavily and and render it completely pointless. And I'll show you an example of what this means. So on the left here is that jibin method that I presented earlier worked fantastically well if the interface is well resolved, but as soon as you have cyber scale.
A thickness involved, the the algorithm kind of really catastrophically fails and start throwing errors everywhere which propagating the domain and and and your simulation really almost, well essentially in that case stops because I I just have way too many interfaces everywhere. So it kind of hints at the fact that.
There's really a balance to be find between cheap, affordable reconstruction method and which are robust to, this type of of features. On the right here you see some some of the work that that Parin has been doing in my group, where we use some form of interface blending, so we have different sensors for figuring out where the.
The interface is is behaving robustly or not and and and merging all those together. And so this is a result of this approach which does a very nice numerical breakup of the, of the interface once once you go below the grid size, you might tell me, well, this looks exactly the same as this move to method that I just had in the previous slide. Well, the, the key here is that this, this interface blending approach is.
There's one to maybe even two orders of magnitude faster than than the mother method, so it's it's really something that would actually be pretty viable to use in in in actual gas liquid for simulation. So this is still early days. We we're still exploring the limits of this, but but we're hoping we'll be able to present some, some very nice results at, at the next meeting with the.
Just aware that one of the nice features with this reconstruction is that, well, because the reconstruction is curved, we directly have access to so we don't need an extra step to compute it. And in fact, we can directly integrate it on the on the clip surface using this this integral here.
If you couple this with the flow server, one of the 1st test cases you do is a laplastic librium, so this is essentially a droppet translating in a periodic domain, and you measure the magnitude of the paracetic flow currents that are induced by the transport and curvature calculation errors, and so on the.
Left here we have the the sort of traditional click method on the right is the PIPEC. The only thing that changed between these two simulation is the type of interface reconstruction, everything else is otherwise adventical. And you see that essentially the the magnitude of your paracetic for occurrence is is pretty much decreased by more than an order of magnitude by switching from the from this piece was linear to piece was.
Quadratic interface reconstruction representation, so I'm pretty happy with this. So I, I I realized I've spoken for a while now, so I'll I'll try to wrap this up. What's next? So where do we stand now? Well, we have proved now that, you know, high order of is enabled by these piecewise quadratic interface reconstruction, but, but we're not quite there yet, right? And so in particular.
This hinges on the fact that the CFL number has to be low enough, so I guess one of the main thing we need to do is quantify what this means for the CFL to be low enough, maybe formalize it in a way that that could represent equate some form of CFR constraints. We're still making progress on on solving the backward problem in a way that is robust and and cheap.
So so hopefully we'll have some nice results to present. Another current limitation is that none of this is GPU friendly and that is due to the the data structure and the and the algorithmic implementation that we have. So I think we'll have some, some thinking to do on how to to restructure this so that we could extend it to, to, to GP architectures and and make use of of massively powered resources. But this has huge potential so.
Obviously the convergence estimation of curature is is a big plus, but on top of that, one thing I haven't really mentioned much is the capability for service reconstruction, which is really enhanced with these higher representations and so in particular, we're currently working on extending this to, to be able to capture service ligaments. On the left here you see a simulation from Austin which.
I'm sure lots of IDLAS people would have seen before where where he represents cyberscale films with with two planes. Well, we could do this with high older surfaces and we could also do ligaments. And actually, I just want to show you some some brief results on that. So you could stretch to represent cyber scale ligaments but that is subject to some pretty significant random errors cause you need very high coverages. So instead what we did is that we extended the approach to quadratic cylinders, and what I mean by this is extruding ellipsoids sorry extreding ellipses or extreding hyperbolas in three D.
And this was work done by my my Studental Basquel, who was here last year, and we've essentially reached more or less the same results as we did for PPEC, so machine accurate robust, it's fast ish. So I I managed to sell it to you earlier that six times flick was was cheap. I'm not gonna try to say that 30 times flick is cheap, but at least, you know, this is something we want to limit to specific cases in which we have service ligaments, so it's it's hopefully something you would only use in, in very specific parts of your sort of.
So that should be a cost we might be able to afford. And this is the type of of cyber scale cylinders we're able to represent now with this quadratic cylinder intersection routines. So, I guess this brings me to my conclusion. I'm already, I think almost at the end so I'll I'll just put those up and and and leave the.
For question. Thank you for, for attending.
All right, thanks so much Fabian. I really enjoyed your talk. We are a little after three, so I think let's just all thank Fabian right now and then feel free to stick around and and ask questions if you're able to. Does anybody want to have a question? I'm gonna jump get us started.
Okay.
Somebody under, I should have wrote down some questions I had some as we were going through there, but I didn't run down, I can't remember on the spot.
That's ok. Save them for islands next year in, in San Francisco. Yeah I'm sure I'm sure we'll we'll talk more about this.
Can you hear me now? Yeah, we can hear you. Oh, ok, yeah. So I was trying to unmute but this is a website issue, yeah. So, this is Jelly. Yeah, Nice talk Fabian.
So.
Yeah, I, I look at this and you know one thing I noticed that probably this is a trivial issue for this message just in terms of how to impose some kind of, wall boundary condition like a contact angle all that. I assume it should be pretty much the same as the previous linear construction.
So that that's, that's an excellent question. I I, I wouldn't go as far as saying it's gonna be trivial, but it definitely has potential in the sense that, you know, with, with, with click, the way you would enforce contact angles is is by having.
Some ghost layers and and try to match the the curvature that would lead to some contact angle, right? Here was quite nice is that we could technically enforce the contact angle as a constraint of the reconstruction algorithm, right? Or the reconstruction minimization. It's something I haven't tried yet, but it's very much on the to do list because I I'm curious to see if if that leads to anything meaningful or powerful.
Okay, that's good. Thank you.
I guess one question that I had Fabian, you, you kind of presented this forward problem and backward problem, and it's like we do the forward problem and we get the moments in the in the news in the solid time TM plus one.
And we do a reconstruction and we throw those moments away, and then we kind of do that again. Do you think there's any value in keeping a history of the moments? Would you transport the moments, you do a reconstruction, but you still have those moments? Could you, could you do anything with those?
To kind of transport the moments and and continue on. So your index reconstruction has like multiple layers of moments that you've transported. Some are the reconstructed moments and some are kind of the older moments. Are there any is the I guess the question is like is the reconstruction introducing enough complex errors or complexity that it's worth like trying to skip that step in some degree?
That's that's a very nice suggestion and I must say I I've never really tried or or thought much about it. I guess it would be very so you're trying to get me to add even more variables in my minimization problem, right? As if I don't have enough already it it could well be a good addition to the the data we use for the reconstruction, right? But.
I guess maybe another way to to approach this would be through enforcing some form of realizability of the transported moments, right? Because I guess if your moments are realizable, then you might be able to match exactly what what you transport and if they're not, then there's.
There's not even a a chance that that you could do it. And so realizable volume is is trivial, it has to be dounded. Realizable by center, ok, it has to be in the cell. Realizable moments of inertia, we're still we're still trying to figure out with with my student what really that means and and and how this this would look like in terms of numerics.
But yeah, that that's that's something we have in mind and yeah maybe maybe maybe we just do this and on top of that we keep some, we keep some history of the moments. That's, that's a that's a good idea.
Could the history of the moments help reconstructing the new ones faster?
I'm.
So if you transport them without without taking too much care of their, like that they're super accurate, just accurate enough to give you a faster trajectory to the solution of to the reconstruction of the new moments.
I have the feeling if anything this would increase the cost, but you you you might be right in the sense that it might make the convergence of that non linear non linear minimization faster, and so your your your resulting reconstruction net cost would be lower.
The thing is, so we're very much trying in, in, in my group at the moment to move away from iterative minimization because it's just, it's just way too expensive. So, yeah, if we could find some heuristical way of using this, this history to improve.
The heuristic reconstruction or or let's say linear regression that we use to reconstruct the interface that that could be beneficial, yeah.
Okay.
All right, I guess if anyone else has a another question, please jump in. If not, let's thank Babyan one last time. I will clap for everyone that's muted.
And, yeah, thanks so much. I really enjoyed your talk. Good background information, good detail, so thanks thanks for kicking off the HIGHLAS webinar series and yeah everyone that's on the call look for the emails about the October webinar and get registered for that. We will have a very different talk on kind of experimental methods and engine characteristic kind of efficiency of engine combustion and things like that, so.
Perfect. Thanks. Thanks a lot Mark. And yeah, see you at the next webinar.
All right. Sounds good. Bye everybody.