Martin J. Gander: Multigrid and Domain Decomposition: Similarities and Differences

I'm very happy to be here it's not my first time here I've been several times I enjoyed this workshop very much I think it's on the top of my list of this summer's conferences I went to the tropics are extremely interesting very stimulating I really really appreciate this organization my topic is multi grating domain decomposition similarities and differences so it's the main solvers that I want to talk about and I want to start with similarities between these methods so the first similarity is that both are iterative methods for PD's this you all know the next similarity is both have dedicated conference series to the topic so there is a multi-grade conference it was the 19th copper mountain conference that I went to this year there is also the European multi grade conference it's now called the International multi grade conference and this year it was in China in Kunming there is the 26th American position conference gonna be in Hong Kong so the figure is a lot of things happening in China so that's another point in common both are optimal methods and optimal has here a mathematical definition so this is a copy from the book of tous le and Wheatland definition 1.2 optimality an iterative method for the solution of a linear system is said to be optimal if it's rate of convergence to the exact solution is independent of the size of the system this is something you can check you can prove four methods but it doesn't mean you have a good method you can have an optimal method that is a hundred thousand times slower than another optimal method they just do not depend on the size of the system so optimality doesn't mean it's good so if you have a paper claims there is an optimal method in it when it uses this definition you still should very carefully read the paper and study if it really is a fast method or a slow method so these two classes are both optimal that's what I'm coming to know so this optimality result is based on the laplacian basically the generalization of this method to non symmetric problem advection diffusion for example is if not trivial and to wave propagation problem the generalization is hard for both types of methods for domain decomposition we've seen there is a quite some progress for multigrid it's basically stuck there is no good multigrid method for a wave propagation that I know of there is not a direct relation that I'd like to mention here it's a beautiful relation that goes back to change I'll show in this seminal paper iterative methods by space decomposition and subspace correction in SCI ref in 1992 and he explained to me the idea that he had when he tried to link these methods in domain decomposition of subdomains and in general have many unknowns in each subdomain now imagine that these subdomains become smaller and you have less and less unknowns and at some point there is just one unknown left then if it was a Schwarz method you have now appoint Jacobi method so there's a direct relation between a point Jacobi iterative smoother and the Schwarz method that's the link he exploited in this paper so it's a direct relation Kingdom's method but to explain more about these methods I'd like to start by explaining the methods and for me I had in my title it's a personal view on this subject for me multigrid methods are essentially discrete methods they have already grid in the name they're based on a discretization they're soft they might to solve a linear system au equals F and the system comes from a PD it comes from a discretization of a PD but it's a discrete system and here is a multi grid method the classical multigrid method it's not complicated to describe you first start with some initial guess you zero you apply some smoothing steps here it's new one smoothing steps this can be a point Jacobi with the relaxation parameter otherwise it's not a smooth or it can't be gauss-seidel that doesn't need a parameter to be smoother you do a few steps you get a new approximation then you calculate the residual F minus L um you restrict this residual to a coarse grid you do a coarse salt you prolong gate that solution that you obtain to the fine grid and you add the correction and then you might do some post smoothing new to post smoothing steps like you had three smoothing steps here we get a new approximation you iterate on this prolongation is naturally done by interpolation you want to go from a fine from a course to a fine grid you interpolate the notes in between restriction can be the transpose of prolongation that's a natural choice there also other restrictions and the coarse matrix you can choose a go working approach to use your prolongation and your restriction from the big matrix to get the small one or you can do a discretize that's a standard multigrid method these are all components that are needed how does it work here is the example that I will use throughout my talk it's an academic example it's on a unit square but it represents a small room a student room so you have zero temperature on the walls so it's pretty cold if not well insulated but you have a door luckily there is a door and the door goes into a warm hallway so here the door is warm and this warms goes also a little bit into the room and also this wall here is connected to a hallway so it's also a bit warmer than the outside part and he has a heater in the middle so he can heat the room a little bit so that's the distribution of this temperature in the room it's a solution of a Poisson equation so I apply multigrid to this it's a V cycle I go down to one unknown I use the Jacobi smooth ER I use as damping parameter two-thirds one pre and one post smoothing step this is the first iteration second iteration third iteration fourth iteration fifth iteration now we see with Brazil converged the fast solve for this problem what results do we know about multigrid the first result goes back to federico in 1964 and it's a result that I like very much it's based on fury analysis it gives a lot of insight into the functioning of this method I can see the full analysis like the Swiss pocket knife the Swiss Army knife for numeric analysis it's a small tool you get a lot of information very precise you might not be able to hold a big truck with it but a lot of information this full analysis then in 1979 Hawk Bush came up with a very general convergence estimate for multigrid and he writes the multi great iteration in one step you see there's a UN plus one of my previous algorithm related to un with a matrix G and here you see the components that I had before there is the smoothing the residual calculation the coarse grid correction and then added as a correction so this is exactly the algorithm I had on the previous slide except it's written in one line and I don't do post smooth in here because Hawk Bush was in his original paper didn't do post smoothing he just considered the pre smoothing step and it has a contraction estimate for this method and the contraction estimate says for any contraction rate you want for example you want the air to be multiplied by half in each iteration poked by 1/10 or by 100 any contraction you want there exists a number of spooning steps new such that you get it this iteration contracts in the norm at the rate that you want and Mike agree plan his PhD thesis in 1990 calculated how many floating-point operations you need for this V cycle independent of the size of the problem 1923 now I remember this from a talk a long time ago but I didn't have the formula the precise number so I send an email to him two months ago and he answered and he told me this is a route from my thesis and I said can I get your thesis and then he said I'm gonna try to dig out a copy for you and he sent me a copy and it's very nice because this is for the five point stand so finite difference blood loss and he also asked for other stencils so the numbers are slightly different might be different this is for number one smooth except the full D cycle 1b cycle for unknown it's 23 floating-point operations yeah the smoothing step goes into a factor I think it's two in this so this is a truly optimal result somehow because if you know that this contracts independent of the size then you get a very very good solver for these systems it's very hard to beat the solver like this it's a very very good song now the domain the composition what's domain decomposition the main decomposition for me is essentially a continuous method I don't need a great to define it domain decomposition method all the methods I know I can write at the continuous level at the PD level so here is the Schwarz method the Schwarz method was a tool to prove result in analysis to close the gap in Riemann's PhD thesis at the end on the riemann mapping urine invented by schwartz so here is the original example we've seen this once already in the conference there is a domain that's pick it's the composed of a disc and of the rectangle and schwartz wanted to solve a lot of provocation on this with given boundary conditions and this alternating schwartz method it just looks at the domain as two overlapping subdomains of simple shape disc and the rectangle and then you solve on the disc you don't know what the value of your solution is along gamma one you put something zero for example and you find a solution on the disc then you trace this solution on the boundary of the rectangle you impose it as a boundary condition on the rectangle you solve you have a solution on the right the blue domain it's not going to be the solution of the original problem because you started with the wrong data but you can iterate on this so you see this is the alternating Schwarz method no discretization is all at the continuous level Schwarz proves convergence of this using the maximum principle it's not a parallel method it's a sequential method you solve left right left right but Leo's about a hundred years later introduced the idea to just start on both domains with some data zero and you solve on both the domains then you exchange on both domains you exchange and then it's parallel and you can imagine if you have many subdomains you can use many processors that was the idea of Leon's to introduce this change so this is the parallel shorts method of Leon's this whereas the other method was the alternating method of words so here is an example the same example as before with the temperature in the room I'm first going to do a alternating shorts method so I have two subdomains I have a subdomain on the left left part of the rule subdomain and the right part of 0f overlap between point four and point six and a trait so here is the solution on the left of the main I start with zero like Schwartz then solve on the right solve on the left solve on the right solve on the left it's all from the right so on the left so on the right and so on and you see it converges now do the parallel method so I solve on both subdomains I exchange data solve exchange solve exchange solve exchange solve exchange you see it's a bit slower I need about twice the number of iterations this is because the parallel method calculates to alternating sequences which are just embedded you can see it just goes like this so for two subdomains is the direct relation for many subdomains it's not so directly relation what convergence results are known about these Schwartz methods the first result I told you about it's due to Schwartz it's done by maximum principle arguments in 1869 we also had a quote 1870 in this workshop the talk was given in 69 and the publication appeared in 70 of words then Leon's about 100 years later came up with a completely new tool variational analysis based on projections it's a very very interesting approach to study convergence of these methods then we have China & Rosasco who introduced fully analysis to study methods of this type and then Leon's in 99 also gave a maximum principle proof so you can wonder why why order to maximum principle proofs because there is a an inaccuracy in the proof schwartz there is a study where he was not very careful it's just claimed something is true and Leon's actually shows exactly why it's true so it's in a more accurate estimate and then there's a very very interesting paper by bramble pashya Quang and shoe in 1991 where they have both continuous and discrete the orthogonal projection estimates so they go using techniques of Leon's and all the results in these papers give contraction estimates exactly like for multigrid these are iterative solvers and they contract at a certain rate and the contraction rate is estimated they're all at the continuous level as well except there is also in Brammer porsche gwangju some discrete estimates as well now there's another method that was invented around the same time and this is the additive Schwarz method this was invented by 300 London in 1989 it's a discreet method and they have a two-level condition number estimate so it's a two-level method and they have a condition number estimate that says the condition number of the precondition matrix is bounded by a constant plus 1 divided by capital H over Delta where Delta is the overlap of the subdomains and capital H is the size of the subdomains and usually in an application this overlap will be a few mesh sizes you could also have very general overlap we've seen methods like this when this is proportional to capital H so for additive Schwarz there is a condition number estimate there is no contraction estimate that's the first interesting difference now this estimate is a discrete estimate but I told you this method domain decomposition method can't be defined at the continuous level so does this hold it does you can prove that this estimate also holds if you don't discretize there's no fine mesh you get the same result the estimates can be done they're technically not so easy but they can be done and we have a result in this from 2014 so this estimate also holds at the continuous line then dr. Neumann Lyman and 50 methods we've seen a lot about the 50 methods the Neyman diamond method goes back to a publication by Boogaloo in school italic and Velasco in 1989 and I'd like to give you the method exactly like they describe it it's a method that has two major steps and I just described it for two subdomains like I did for the Schwartzman so how does this method work you have to start with some guests at the interface between the subdomains and there is no overlap so it's not like shorts there's no overlap you take some guests for the solution like in Schwartz zero for example which I call H and then you solve to directly problems with this boundary condition because now you know what the solution is at the boundary so you can solve on both subdomains so you obtain approximate solutions UI and then you follow with the correction computation what is the correction computation doing you take the traces the normal derivative traces of these two solutions you calculate it and you add them and you impose this as an ointment condition to the equations in the subdomains and no source here it's a correction so there's no source and you calculate tube side so each step takes two direct results and to nine months offs and once you have these solutions of side you take their typical it races you saw them and you add the relaxation parameter you can calculate the new difficult race at the interface so you see exactly how it works you start with some data you solve directly you do now man you solve Neyman you have a relaxation of a new dealer theta that's the Neyman diamond method we've heard a very nice talk already by constructs a via today about the Fetty method the fakie method was invented in a quite different setting from the way we've seen it it was invented by imposing the problem as a minimization problem it comes from a variational formulation so you don't look at the Euler Lagrange equations which are the Laplace equations you look at the functional to minimize you cut the domain into two you say we want to minimize only in the subdomains but then it's not consistent at the interface so introduce Lagrange multipliers to make it consistent that's what how the method was invented but this method is essentially the same as this method but in the other order you start with the Neyman data and Neyman source and then you do direct Corrections like it's done here for Lyman and you have a relaxation step so these methods are very very much related [Music] so how does it work on my example I just showed the Neyman name and the 50 would be very very similar I start with the non-overlapping decomposition of smaller sub domain and the bigger sub domain and I start with a zero approximation and I solve the Didache resolves that's what I get now I take the normal derivative at ease I sum them I do now and corrections and I add them and to Adira call solve again so that's the next step Wow you see that looks almost like the solution so I do another one so we saw it overshot a little bit are you another one not much changes it's very fast like so what happens if i take subdomains that are exactly the same size so I cut exactly in the middle and I remind you this is a Poisson equation same operator same subdomain size same boundary conditions theoretically all around so I do one correction and I have the solution this is a direct solver for this configuration the direct solver if I do once more it doesn't change that's nine one nine one fifty has exactly the same properties what is known about these methods what is known about nine one nine month there is a result by dark and italic in 1991 they have a continuous analysis for this method there is a continuous analysis that shows that the method contracts but there are no cross points in this analysis and with cross points there is discrete analysis and then there is again a fundamental condition number estimate that goes back to Mandal and Brezina and you see it looks like the Schwarz estimate the condition number of the precondition system is bounded by a constant times 1 plus a log term of capital H of little H this is a discrete estimate and pick domain size and the mesh size coming in this estimate and there is a similar result also by our model it was very fundamental contributions here just a little later 450 a very similar condition number estimate and here it's really appropriate to mention that excel has worked tremendously on moving the theory of these methods forward he was instrumental their papers with all of papers with the Rhine box paper with with Max Trier there is a lot of contribution that we heard axel talk about the 30 methods today as well there's a lot of important analysis for these methods again I'd like to emphasize these are condition number estimates whereas only in the first paper without cross points there was a contraction estimate so we've seen two differences so far between multigrid and the main decomposition for me multi great as I told you is defined that the discrete levels the discreet method TV methods are defined at the continuous method to continuous level their continuous methods that's for me a fundamental difference the second difference is that we have contraction estimates for multigrid that all contraction estimates we have contraction estimates for some domain decomposition methods but there are a few of the very popular important methods that do not have contraction estimates they have condition number estimates and additive Schwarz is one of them no I'm on I'm on wait cross points is one of them and also fakey with cross points is one of them so why why do we have condition number estimates for these methods well as the methods as I've shown you they just contract their PD methods they just contract why is that I will get to that yes yes so let's look at additive shorts for the temperature in my room so I do additive Swart with the same two subdomains I've shown you before alternating and parallel ones are actually during first iteration second third fourth fifth six seven eight nine you see it's not converging there's not a convergent method additive Schwarz is an iterative solver is not a convergent method it was not advertised as a convergent method is was advertised as a pre-conditioner but I think this is important to know so this is very different from the solutions I've shown you before what happens if I use more subdomains let's say I use 25 subdomains I solve the same problem and now I indicate the subdomains here so you can see I have 25 subdomains 5x5 and I use additive Schwarz in this example first iteration of additive Schwarz you see already here at my door there are over shirts in the overlap I have over shoots here and here I start heating now if you look at the iteration of additive Schwarz as an iterative solver you see these cross points are extremely bad for a dish once it's not converging but we have a condition number estimate so there is a bound on how bad this goes and a crowd of method can take care of this but as an iterative solver it cannot be used like this see AC divergence is very violent 9/10 I get to exactly right I just wanted to say now I learned from Danielle shield about 15 or 16 years ago that there is a linear algebra community but this is well known and you can introduce damping in the method and if you dump the method then you can make it converge so I damp it so I introduce the damping which is strong enough to count these modes so now this is the damped additive Schwarz iteration right up and now you see these oscillations still there but they're controlled by damping so ie trait you see it works I mean 20 iterations I get to about the height of the last red line now what happens if I use the parallel Schwarz method of Leon's on this example because that one does work so I just do it this is now the parallel Schwarz method of Rio's we just iterate no problem there is a proof that this converges by maximum principle and you see in 20 traditions it converges to a much more accurate solution because the damping that you need to damp these divergent modes it also harms the convergence of the method but it's very important adesh what was never advertised to be used like this it was advertised to be used as a precondition but I think it's important to know what happens if you use Krylov acceleration here is a comparison here is the number of iterations here is the error and here are two curves there is a curve which are called wrath and there is another curve which I call as Erika Schwarz or restricted at before so you see the Krylov method needs to do some work to remove these modes that additive schwartz produces but it converges it's a good method but if you use wrasse it takes less steps so it's beneficial to have a method that converges as an iterative method in the situation for the number of iterations so again it's very important additive Schwartz is a pre-conditioner it was never meant to be a stationary solver but if you have a good stationery show but like wrasse you need less iterations now I should tell you about Russ Russ is the correct discretization of the Leon's parallel force method if you apply brass you get exactly the Schwartz method at Leon's design so this is the method that you've seen on the previous slide it Rast respected Irish Wars if you want to read about this it took quite a while to understand this this is a paper which is very easy to read just two or three subdomains and this is a paper that proves it in full generality and I'm very happy these papers are outside that more and more I think people want to know about this they want to know why rice is better than Ali - words now let's go to nine one nine one I also told you now I'm on line one has a condition number estimate no contraction estimate as soon as there is a cross point and here a shower result that's due to my PhD student faisal actually implemented this myself this summer as well and it was just before coming to this conference i wanted to add the course correction that I want to show later I did not manage to I'd finish so I showed the calculations that he did this is now not my room but it's also on a square it's cut into four sub squares and I run the Norman Norman iteration with a cross point and I start with this initial guess it's a smooth initial guess and the solution is the pump that we want to calculate and I ran it as an iteration stationary iteration exactly like I described it now you should look at the scale here as well why I iterate this is the first iteration second iteration the reiteration look at what happens at the cross point fourth iteration we're at the scale 10 now 53 or the scale 30 this is also not converging that's why there is a condition number estimate there is no contraction estimate but here the situation is much much more severe because this this method is not well posed in any function space if you look at it at the continuous level you cannot find a functional setting such that this method is well posed and the problem is that you have moment traces which gives you new data collected races which are discontinuous at the cross points which you will solve and then take Norman data from a discontinuous theory or trace on the boundary this is going to have a scene rarity for which you impose a directly trace again after you've averaged this method has no chance to work as an iterative method but again this growth is controlled this gives you the long-term in the condition of resting so that's why these methods are preconditions and are not stationary iterative methods but I like these methods at the PD level so I wanted to write this method at the PD level and it's possible but you have to add a correction and here what we do is we do a course correction not for scalability it's to fix the methods and then you get a method which is well-posed h2 this is a very recent result even if you have a cross point and if you run it iteratively exactly the same example but now with this law is a local course correction it doesn't do global communication then you get a method which converges again you can wonder is this important the kraut method fixes it anyway and it's true it fixes it this is without this local course correction and this is with the local course correction so it does help we take less crowd of iterations you also pay something for the local course correction it's not free but it's not an expensive correction it's not a global one so now there is a method which I have not shown you yet but it's a method that has been mentioned very often here in this meeting it's the optimized towards method and it's a method that goes back to Leon's as well it's amazing what Leon's had in these three papers in the early days of the didi conferences just amazing the idea is a very simple one you remember much words did he just took tarantella data at the interfaces and what this Leo's proposal is you take a combination of tearing out no human data at the interface that's the only difference you can do it in parallel like I showed you for Leos you can do it alternatingly like the alternating Schwartz method it's a new class of methods there's an early result by frederik notice that shows that if you have two subdomains and if you use for this parameter here at the interface which act on the delegate race the Terracotta Norman operator this is an expensive non-local operator you get convergence in two steps it's a direct order this is like the name alignment method if you're exactly in the middle but it works for any position of the interface and it works with and without overlap so this is a very general early result and in addition he shows that if the decomposition is a strip decomposition into J subdomains then you have to iterate J times with this parallel method and you still have a direct order of the generations you have the solution only more recently we realized that in this case when you use the derivative Norman operator is the parameter and if you do an alternating method so you go from left to right and then from right to left it also is a direct solver and in that case there is a very nice relation of this method to linear algebra because you can show it's a the exact block Lu decomposition of the matrix but written at the continuous level so these there are two norman operators if you wrote your matrix you did a block Lu decomposition you would see sure compliments appear and these are the short compliments which are here to do her to normal operator so how does this work on my example here is first an example with two subdomains I do an optimized towards method so I just have a number and we have formulas what a good choice is for these numbers for many many equations we know many equations so this is a Poisson equation so we know what the formula is this is the first iteration it's a parallel version second iteration third iteration you see this is a very fast method cost per subdomain is the same it's just the robin condition instead of a dirichlet condition now I'd like to show you what happens at the cross point with this method so here I show you the error so you can see and there is a decomposition into for subdomains and there is a cross point and I've shown you 30 and Neyman Norman as an iterative method they have a problem at the cross point so here is an optimized Schwarz method we iterate you see these are contracting methods as a stationary iteration and they're well posed in the function space now this idea of optimized transmission conditions has been reinvented over the past years by many many researchers and scientists here is the list the last one I discovered it's also an optimized Schwarz method is an algorithm by bank and GMAC in 2001 that has an analysis in 2008 by bank and Vasilevsky then there is an elliptical earth complete Lu this is where a relation with this block Lu decomposition was made but I didn't yet see it's exactly the block and the composition in the sweeping sense then they were to sweeping preconditions by engquist in being they're optimized Schwarz method where is a PML to approximate the direction to normal operating there is a single layer potential method by stock in 2013 there is a source transfer method by Chiang and Jana and when I met Jimmy Chang he's from the Academy of Sciences of in China I asked him why they invented this method because I didn't understand the method yet and then he said we couldn't understand the sweeping precondition of inclusive doing and so they came up with their own method called its source transfer domain decomposition there's the method of polarized traces from the MIT school and it takes a mind-boggling effort to understand these methods because they're invented in different settings for example this one here is invented the source transfer is in a DD setting but with greens functions formulations the sweeping pre-conditioner is invented in the discrete setting not having subdomains but just grid lines it's very difficult to link them and it took us about five years to show that all these methods are based on the same mathematical principle at the continuous level it's this optimal Schwarz sweeping method and at the discrete level is the plot Lu decomposition the exact plot how do you get the composition is the paper that just appeared now in Siam review a simpler introduction is given in 2006 but there are no clue yet about all these were the relations these are we do not bear it so let's look at my room what happens if I do one of these sweeping types so I need more subdomains now so suppose I have four subdomains in this direction and I sweep in the room and I use the optimal Schwartz method the one that has the Dirac Latin human map so I start on the left sub domain and I go and then I come back I must have the solution and that is the solution so it's just one forward solve and one backward solve that's an optimal Schwarz method now I can use an approximation for example I use a robin condition with the P that I know what its value is so I go forward and backward it's pretty good but this is not the solution it's not the exact plot coming the composition so I go once more and you see it does correct a little icon once more it still corrects a little and now we're converged so these are very fast methods and the most of these methods were invented for Helmholtz but it has nothing to do with hormones the fundamental structure is the block tell you the composition of the optimal Schwartz method PD is is a completely irrelevant let's look at the classical Schwartz method alternating sweeping what is not very fast you can see but I remind you this is also an optimal method if you fix the overlap and you'll find the mesh the number of iterations is not depending on the number of unknowns but it's a much much slower method much much slower I think I stopped at some point I hope might be a day eater so yes I did okay so I've shown you two main differences so far I've shown a discrete multigrid versus continuous DV and I've shown you contraction estimate versus condition embracement there is a third major difference for me it's course spaces course bias Corrections and I want to explain this in the remaining time so what's the course correction that's needed in a multi grid method in a multi grid method the course correction is based on what's called an error equation and here you see the error equation I have a you and I add and subtract my approximation um and the error is U minus UN so I get a times the error is a times my approximation equals F and so I can put the egg you win on the right and I find what's called the error equation a times the error rates depend equals the residual and the course correction is based on this error equation in multigrid so what do you need to do with this error equation you could just solve it if you solve it and if you add this error to the solution to the approximation you have the solution you're done that's the best you could do but that's as expensive of solving a problem so you need to get a cheaper way to do this and what you do is you looked at some chekhov is smoothing steps which you do before you do the course correction and you look what's left in the error this is what I plotted here for the room example so if I do with three Jacobi smoothing steps with the damping parameter 2/3 that's the best damping on the high frequency for to copy then I get this thing here and you see it looks rather smooth I started with a random initial guess so there's a lot of oscillations initially it looks rather smooth and here is the residual and so the idea in multigrid is because its a lapras equation you have this property from Jacobi or gauss seidel and you approximate this error on a coarse grid and you solve this problem on a coarse grid that's a good course correction for multigrid now how does this look like in the main decomposition so what does one classically do as a core space there is a quote from the book of the salient Midland as well the subspace V zero is usually related to a cores problem often build on a coarse mesh that's what we usually do like in multi grid we build the coarse mesh with a coarse solve and the residual equation and then we interpolate but look at what happens in the error that I obtain if I do one parallel fourth iteration with a random u0 you see it also looks smooth but it's not smooth in the same sense multigrid you see in the subdomains it's extremely smooth these functions are harmonic in the sub domain but routine subdomains there are jumps and if you look at the residual it looks completely different from the multi grid residual this is the residue the residual is zero within every subdomain why is this because domain decomposition methods solve the problem in the sub domain there is a zero residual left the only place where you can have residual is what the domains are connected so you see what the shape of the error is and you also see what the shape of the residual is from this experiment so I'd like to show you what the error equation is for domain decomposition like I've shown you for multigrid so we know that any domain decomposition this is 15 I'm on line munch Ward's optimized Wars no matter which method you take it solves on the subdomain the problem it's Lu J M equals F that's what these methods do and so if you define the error like I did it's for multigrid it's the solution minus the approximation then you see if you apply the operator independent of the PD if you apply the operator to this expression then you get L times u which is f and you get L times UJ n which is f so you get zero that explains why the error is always harmonic that's why it looks so smooth in each subdomain here it's a harmonic function which means it solves your homogeneous PD no matter what your PD is it solves the homogeneous PD so from this you can immediately define what an optimal course P is not optimal here in the sense of scalability or in the sense of number of iterations optimal in the sense there is no better space this is the best space you can choose so what's that space you just need all piecewise harmonic functions in each subdomain because that's where your error lies in no matter what your DB method is now what's the equation what's the error equation that one has to solve from the approximations that we calculated we can calculate the directly jump and we can calculate the Neyman jump because we know the UI and UJ on each interface and so what does the error satisfy the error that we've seen in the previous slide the error satisfies the homogeneous equation in each subdomain and at each interface it satisfies a jump relation this is called the transmission problem we've seen transmissions problems earlier in these constants as well so if you solve this problem the solution lies in this space and then you have the solution of the DD method after the next step this is the exact error equation like in multigrid that's the one if you solve this one it's a direct method yes it's a very very good question so this works for any method overlapping or non overlapping but the correction in the core space is done on a non overlapping setting so use the subdomains solutions in a non overlapping decomposition and the jumps are posed at the non overlapping interfaces now this problem can be simplified if with an ion lineman method you know the derailer jump is zero if it's a xxx method you know the Neyman jump is zero so you get a reduction in the size of the problem but this is the exact equation that the core space should solve so how does such a core space look in two dimensions here is an example there are nine subdomains 3x3 I construct overlapping ones but like Rob said I also put the lines here of the non overlapping composition because that's what the course space correction is applied and so the optimal course space in this case it's a very big core space because you need all the traces to the left of the sub domain boundary and to the right of the subdomain boundary and then you need for each of these degrees of freedom the harmonically extended function to have all discrete harmonic functions and in that space you find the correction in one go it's a direct order now that looks very expensive so I'd like to first show you an experiment this is again the example I had it's an optimized Schwarz method I could do this with any other method as well I have one Cross point and I showed me only the error again we've seen this already but now look carefully what's left when I enter it look at these functions we know they're harmonic because it's the error that means they satisfy the equation but they have a very typical shape their q1 functions why do I know these are q1 functions because all these Schwartz theorem I'm alignment methods are smoothers this is somehow the bootstrapping process you've seen in algebraic multigrid you see I've run the master three four times this is what's left it's a q1 function I know this is a priori for this Laplace equation so I know what it is so I just put the q1 function in my course grid so I put for q1 functions that I correct in each iteration that's the first correction you see these q1 functions are gone but you see a second function appearing look carefully at the interface here here and here I started with random so there is noise but this is a smoother so the second iteration the iteration these are sine functions it's the first sign on this edge the first sign on this edge the first sign on this edge we've seen in the talk by Axel the sine is the solution of an eigen value problem along this interface for the laplacian it is a sine function so this space could be enriched I can put the sine function then this one is gone as well but then what's left is the second sine function if I put the second one as well and reach the core space there's a third sine function it's just a spectral basis of my very rich space that I've shown you which has too many points so if I compare various q1 spaces no but I could have a q1 core space that has a degree of freedom in the middle of every domain I could have a q1 space that has a each cross point a degree of freedom or I could have four unknowns in each subdomain Sir Richard space or I could put four points around each cross point we've seen which one approximates the best space it is this one because I think if I put all the nodes that's what I would need that's why the q1 was a good function for the La Crosse and this is see immediately if you run this as an iterative solver so this is an experiment now we're sixteen by sixteen subdomains it's the parallel Schwarz method of Rio's I implemented it as Russ here is the number of iterations here is the maximum error and you see if you put these q1 functions you get a much faster convergence and if you put any other option of a q1 core space this is really an approximation of this optimal core space here is a comparison that was calculated by search on creaking in who is a researcher at Idris and I'd like to mention how I got to know him I had never met him when he sent me an email and he said don't you think one should do a careful comparison of these methods that one person implements in one software package and test oh yeah great that would be fantastic and he said this post allows him to do this on the side so we would like to start and then we start to have an email exchange and he went over here and after the year he had comparisons of all these methods and I realized the guy is very good he's a very clever guy he understands the methods he can call them and so I invited him for a seminar and we discussed the results and we're working in the meantime on non symmetric problems because all these ideas apply to non symmetric problems but here our results for a symmetric problem so this is all results that were done in pets if I search on creaking in from Idris and this is a scaling experiment we have 200 400 800 up to a thousand subdomains which means process it's a petty implementation and these are iterative methods no crowler for the moment so I cannot show additive Schwartz on this so this is just the parallel Schwartz method and the top line is the parallel Schwartz method with a classical course correction which means one unknown person terrain in the middle then the next line in blue is the classical parallel Swartz method but now with an optimized course pace which is the q1 course base which has four nodes around it's just on a square and then here is an optimized Schwartz method with a Robin transmission condition but the classical course pace and here is with optimized transmission conditions and optimized core space so you see all these methods are optimal but some are more optimal than others you see but multigrid is still faster the number of iterations multigrid is faster in number of iterations on this problem yes and this Laplace problem you can actually do it's very simple but for a general problem is much harder yes yes yes yes this is not no Crowell of Normandy first station or iteration so what happens if you use Krylov everybody becomes faster because Carla finds a better polynomial and now I can also show additive Schwartz so in red there are lines added for active Schwartz and additive Schwartz is symmetric I use counter the gradient whereas for the other methods are non symmetric I need to use generous but you see number of iterations even though you have a better cry love my food additive Schwartz is paying for these cross point divergence modes so number of iterations is slower than the parallel Schwarz method of Leon's optimized Schwartz is still faster but here I use Krylov as well so it's still multigrid Aziz's much lower iteration count so let's look at CPU time because now does the same guy who is implementing it's the same guy who's testing your CPU time same experiment you see the parallel ports MA there's an iterative solver with the classical course correction it's a by far the slowest method if you put the peddle course credits already much better method if you use optimized transmission conditions but the wrong course grid it's still better and with the best course credits better than all of them and now you see that multigrid is paying in communication time that's a well-known problem of hyper it's happens also with the geometric multigrid but I think one could tune there if I use Krylov on top then all the Schwartz methods come down you see krylov's really helps a lot now they're all better than the multigrid method that we tested in wall clock time and if you blow it up you see the same hierarchy basically remains so you have the optimized Schwarz method with an optimized course pretty much the fastest whereas classical course click classical transmission conditions are the slowest here I just run the solver in the package it's it's it's the no we have the geometric as well it's just not on this ball but I have it as well it's very similar yeah no no this is really this the best that I can find but but I would probably need to work on tuning it's the standard package yes yes yes so I did both geometric and algebraic yeah I don't want to say multi grid is slow this is behind by no means what I want to say here I mean this is I told you at the beginning for me multi great for Poisson is very difficult to beat this is really good solvent this is just how it came out in this comparison and it was meat both packages yes previous flight okay it's just a q1 here no no other enrichment what could further enrich please just the q1 I've shown you know you don't work harder at all if you look here it goes from this blue line to this blue line if you have an a synchronous method you cannot use crown off oh yeah poor optimized Schwartz it's much less this is classic reports yes yes you could but it's not more complicated it depends where you put your notes this is really not complicated you can choose your notes in the middle or you can choose your notes at the cross point this is really not more complicated it's the same code same right it's really not complicated but the important thing is you could enrich it as well so I showed you results with enrichment on a more difficult problem so this is now a problem which has high contrast like axel had it you see there is high contrast that crosses interfaces this is an example where we still have rectangles but you can imagine now that domain can have an arbitrary shape you will still get not a q1 function for the function that's one at the cross point going linearly down to zero if it was a laplacian or here it's a multiscale element that you get when you run and look at the error it's just a multiscale element so it's a multiscale element that Ivan Graham has introduced as a course face already which can be used or you can enrich it so with one function per interface two functions three functions four functions and you observe if you go to two functions things improve in iteration numbers as contrast goes up but it's not robust if you go to three functions certainly its robust its robust because you have at most three inclusions that cross the interface so you need three functions for three inclusions but you can do it adapt if we like axel said if there is just one thing crossing one function is enough and so it also robust if you do adaptively one function on one three functions if there are three this is an example for a multi scale core space we call this Shems back to the harmonically enriched multi scale cost core space and it's exactly what I described for the laplacian it's q one functions plus sine functions and here it's multi scale functions plus solutions of local interface dimensional eigenvalue problems so to conclude I would like to remind you of the three differences that I tried to show the first difference was discrete versus continuous the second was contraction estimates versus continued condition number estimates and the third one was coarse faces you can find preprints on my web page and I'd like to advertise it's now the end of the summer is probably the last time I advertise this at the beginning of summer I complained with Siam that we only got 5 copies to give away of this my little PD book it's a little PE book under 50 pages and then a little bit Greenspan said yeah it's ok but if anybody wants a book they just send me an email so this is my little PD book I like it very much it has four main components it describes the four main discretization methods for PDS there is finite differences there is finite volume the respective methods and there is finite element methods for each of the methods there is the precise historical way how the method was found then there is a complete convergence proof they stripped bare down to the bones so you can basically just do it on the board there's no result needed from any analysis book it's worked out in the simple small setting but it contains all the pieces that are needed in a general setting as well it's just simplified and for each method there is a little code that you can try to play play with and if you send a little bit Greenspan an email it seems she sends you a free copy so thank you very much [Applause]

Channel: Centre International de Rencontres Mathématiques

Views: 1,100

Rating: 5 out of 5

Keywords: Cirm, CNRS, SMF, Mathematics, mathématiques, Marseille, Luminy, Centre international de rencontres mathématiques, aix-marseille, université, AMU, aix, marseille, university

Id: 8fcFI6CgY-4

Length: 56min 56sec (3416 seconds)

Published: Tue Oct 08 2019