Modeling X with PDEs

I just said that I work in the numerics group and I typically work on PDE modeling and the finite element method. And today we're going to talk a bit about PDE models. So the most important thing about PD modeling, if you want to do PD modeling in the Wolfram language, is to look at the PD models overview page. So this page has a collection of available models that we have, and I'm going to explain a bit of how that's structured. So here you see various fields of physics like acoustic, electromagnetics, fluid dynamics, there's a few more, heat transfer, mass transport, structural mechanics, and so on. So some of these fields are fairly fleshed out and some others are we're just getting started on. So if you have a field of physics that's fleshed out, so something like acoustics, for example, you'll have a guide page and the boundary conditions page. Then there's going to be some monographs on the topic. So acoustics in the frequency domain, the time domain, and then we'll have a section on acoustic models. So those are extended application examples. And these examples typically go beyond what we have in the documentation. So they're much more. So if we look at this acoustic muffler example, they have much more content, more in depth, and typically resemble some real world type of application that you'll find in the paper or some application from industry, things like that. So this is just a bit of an overview of what that's like. The important thing is if you go to this acoustics PD boundary condition page, you'll get an overview of all the functions that we have in that specific domain. You get the links to the monographs that explain the theory behind these topics, a list again of the models that we have, and there's a small overview. So, for example, the acoustic muffler we just looked at is a 3D example. It's in a frequency domain. It uses this type of boundary condition. So that gives you a bit of a guidance of what application model you might want to look at if you have something that you want to model yourself. So this is really the important starting point if you want to do PDE modeling in the Wolfram language. Okay, so the PDE modeling started in version 12, where we started with acoustics and heat transfer. In version 12.1, we added mass transport. And we needed to do these two steps to get an understanding of how PDE modeling should work in the Wolfram language. And once we had this understanding, we started to make a framework and started to come up with basic Building blocks, which we've added in 12.2. And the acoustics, heat transfer and mass transport were then adopted to fit into that framework. In 13 oh, we added a little bit of electromagnetics. In 13 oh, we added A lot of solid mechanics. And in 13.1 we started with region symmetry and some more solid mechanics. And that's, I'm going to show you going to be showing examples about that. And in 13.2 we've added more region, region symmetry functionality. Okay, what, what's the essentials about PDE modeling? So if we look at this diffusion type of problem here, what we always have is a dependent, at least one dependent variable. So in this case it's U. We also always have an independent variable X. In this case, this is the spatial variable. We might have time or frequency variables. That's not always there very often. We also have parameters those are material parameters, like a thermal conduction coefficient or things like that. So we separate our PDE modeling into variables and parameters. Then for the different fields of physics, we try to come up with something that's called XYZ PDE component. So the heat transfer PDE component, that function gets the variables of that field and the parameters that we'd like to specify and it's going to produce the partial differential equation. We also have functions that are called something like XYZ condition. They get a predicate where they're active and again the same variables and parameters like for the components and they generate a Dirichlet condition. And we have XYZ value, again a predicate where that boundary condition is active and the same variables and parameters and they going to generate Neumann values. This is, this is kind of the basic workflow of how this works. So let's look at an example. What I have here is a heat equation that I'd like to solve over a cylinder that has, that's hollow, that has something taken out in the middle. So in the middle, the red part that has a heat flux. So the energy is going out through that part. And then the top and the bottom we have a fixed temperature. So how do we do this modeling? So we start by generating the geometry. In this case it's a region difference. So we have the outer cylinder and subtract the inner cylinder. That's the geometry of the object. We generate a mesh of this. Now we visualize the wireframe of that mesh. Now we come to the variables and parameter specification. We do a 3D problem. It's a temperature problem, a heat transfer problem. So we have T as the dependent Variables and XYZ as the independent variables. Now we need to specify the material parameters. We can do that by either saying that the material is aluminum and the solver will extract the actual data it needs from that element. Data element. Or you can, alternatively you can specify the thermal conductivity manually. So either way, if the element is available in some of the entities that Wolfram language has you good to go. You can use that. And if not, you can specify the parameters manually if you want to. Okay, now what we do is we call the heat transfer PDE component function with these variables and parameters, and that's generating us this PDE operator. One thing to note is that this operator here doesn't have an equal sign. So this is. It's not equal to anything. It's just part. It's a component of the pde. And the reason for that is we can now use these components like Lego. We can build things up from these basic components. In this case, what we're doing is we're still adding the boundary conditions. So we say a heat flux of 5,000 watts per square meter is going through this inner boundary. We specify the heat flux, boundary condition, heat flux value. So this is going to return a Neumann value, specify where it's active. So the inner boundary vars parse. And this is the specification for this specific boundary condition that we have a heat flux of this 500,000 watts per square meter. Right. The second boundary condition we wanted are the temperature on top and at the bottom. And we specify 273kelvins heat temperature conditions. So it's going to return a Dirichlet condition. We specify where this is active, bars and parameters and the specifics of this boundary condition, the reference temperature. And then we'll plug this into ndsolve. So we have this operator, set it equal to the heat flux and specify the Dirichlet condition here and we can solve this. Takes about a second to solve. I mean, from now on it's basically the standard way how we work with partial differential equations in the Wolfram language. And we can of course visualize this. So this is a 3D visualization of the solution of this specific problem. Now I'm going to try to say expand a bit. What's new in version 13.1 and what's going to be new in version 13.2? One thing we added is region symmetry, the capability to generate partial differential equations that work with axisymmetric geometries. In this case or in this, in the case we've seen just now, we have a symmetry around The Z axis. And what we're doing is we're having a truncated cylindrical coordinate system and use a 2D version. So this grayish area here is the 2D version. And then we can rotate around the z axis to get the full solution. But this is what region symmetry does. We specify the 2D geometry, that's a rectangle. In this case, we set up the variables and parameters. So now we have a 2D problem. It's R and Z are the independent variables. And the only thing we do is we take the parameters we use for the 3D example and just specify the region symmetry is axisymmetric. That's it. And now the heat transferred PD component where we plug this in, generates a new kind of a PDE. So this 2D axis symmetric PDE, basically what it does is you have the divergence gradient of the temperature temperature of R and Z, not of theta of R and Z. And then we take the cylindrical coordinates R theta and Z and again of the divergence of the divergence of R theta and Z. And that's exactly the same result what the heat transfer PD component does once you activate it. So this is kind of what's mathematically happening behind the scenes. Oops. Again, we specify the boundary conditions vars parse and the predicate where it's true, same heat flux as before. The conditions where it's true vast pass and same reference temperature as before. Same thing. We solved the same problem. The operators again, the flux and the Dirichlet condition, it's much faster now. It's about an order of magnitude faster now because it's, it's 2D only and not 3D. And we can visualize the solution. So this is the axisymmetric solution. This is a cut through the 3D solution. And we can also do a rotational plot of the solution. So the point is now that the solution in 2D, the axisymmetric case, is much faster than the full 3D solution. Okay, so where do we have this region symmetry operation? So we have it for the basic building blocks for the PDE terms like PDE diffusion, conservative convection derivative term. In 13, 2 we added to the Laplacian PDE term. We have it for the heat transfer PDE component, mass transport. And in 13.2 we've added it for the acoustics, Poisson, Helmholtz and wave PDE components. We don't yet have it for the solid mechanics PD component that's coming in 13.2. In fact, just today I've made an initial commit for this axis symmetric stuff in for solid mechanics, we also added application examples that show the usage of this in a disk break, tubular reactor, spherical capacitor, solute centrifugation. And in 13,2 we've added the Helmholtz resonator example for the axis symmetric operators. Okay, that's region symmetry. Now I wanted to talk a little bit about solid mechanics. So solid mechanics is another field that we've added. And for 13.1 we've added hyper elastic material models. So we've added the Samba Nau Kirchhoff model and a Neo Hooken model. And we have the capability or you, you can add your own material model if you would like to do so. That's also possible. And also going to talk a little bit about Rayleigh damping. All right, so what are hyper elastic materials? So hyper elastic materials are materials that can withstand extremely large deformation and are fully reversible, like rubber. So you can pull rubber and when you let it go, it goes back to its original deformation, original form. So in the example I'm going to show is we have a block of rubber that's fixed on the left hand side and we have a force on the right hand side that pulling to the right. Now we have a coupled problem. So we have the displacement in the U and then the V in the V direction. So X and Y direction. Independent variables are X and y. It's a 2D problem. And the parameter specification now is that we specify the material model is a Neo Hooken isotropic material material parameters, lame parameter shear modulus. This is a plane stress form. So it uses the plane stress formulation. And for plane stress we also need to specify a thickness of this. So the exact details of what these parameters are and how to use them are in the documentation. And there's a large monograph actually specifically on hyper elastic materials that explains the theory behind all of this. So what I've done now is I've set up small help function. Again. We have a solid mechanics PD component vars parse, a solid boundary load value on the right hand side x equals 1 vars parse and we have a force of quantity F. That's the argument to the function that goes to the right and zero force upwards or downwards and we fix the rubber part on the left hand side. That's basically the help function I've written. This is one side comment I wanted to make. I mean, why do we actually write these PD component functions? Well, it turns out that for some cases like this nonlinear hyper elastic material law, these expressions are fairly complicated and it's Much easier to have them generated than to, I don't know, derive them and type them out or things like that. So that's why we have this functionality, basically, to make it easier for customers to set up their models and less error prone. All right, so how does this work? Now, again, we'll plug it into ND solve. We have the function neo hook, and we specify that we have a force of 300 Newtons going to the right over this unit rectangle. And here is a vector displacement plot of that of that deformed rubber. So you see, we've stretched the rubber quite a bit to the right. Okay, so here's another example. Same setup, same setup. I've just increased the force, so now it's 450 Newtons to the right. And ND Solve complains that it can't find a solution. In fact, not only ND solve would complain. Pretty much any finite element solver would complain about this, because hyper elastic material laws are extremely highly nonlinear materials, and the finding solutions for them can be very challenging. Now, what can you do when you want to do when you still want to solve something like this? So there's two approaches, and I'm going to show one of them and I'm going to talk a bit about the other one in a second. So what you can do is you can set up a parametric nd sol value. So you have a parameter now, and the parameter is going to be the force that we're going to input into this material model. And we have a pfun neo hook for this thing. Okay, so what are we going to do? One approach to tackle this highly nonlinear problem. So any highly nonlinear PDE problem is by using a solver restart. What we're doing is we're specifying a maximum force 500. So it's even more than what I've shown initially, the 450. And we want to get to these 500 Newtons in 20 steps. And we start basically by a zero displacement. Now I have a do loop. So this loop in the first step, step one, uses 500 over 20 as a first force. We plug this into the new hook function. And now what's new? What we can do is you can specify an initial seeding. So this initial seeding in the first step is going to be zero. But once we get to a second step, this is going to be 500 by 20 times 2. Because with the second step, and now the parametric solver will use the initial seeding from the first step, from the previous step. So that's a way you can kind of slowly move the solver to the final 500 Newton. So and then you can get to a solution. So this is an, is an example of a solver restart. The other option you have for these type of problems is you can convert them into a time dependent problem into a pseudo time dependent problem and have the time go from zero to one and and slowly increase the force by multiplying it with a time t. So T at the beginning is zero. And then slowly while the time increases, the force ramps up to 500 Newtons. The advantage of that approach is that ND Solve has an adaptive step size mechanism. So it will take the exact steps that it needs to get to a new solution. While with the solver restart we've just specified the steps manually. The solver restart may be a bit simpler to set up. So it's a good initial try. If that fails, you can still go ahead and use the time stepping. Okay, so different, slightly different topic, we've also added Rayleigh damping. Rayleigh damping is a way to make solid mechanics model have damping. It's not a physical model, it's a phenomenological model. And what it does it basically in solid mechanics problems, you have a mass matrix and a stiffness matrix and it finds a combination of the mass and stiffness matrix C, which is a damping matrix. How to so these, these coefficients, alpha and beta can be computed. How that's done explicitly, that's explained in the documentation in more detail. So again I'm going to show a simple example here or fairly simple example. We have a time dependent problem. Now, U and V, T is the time. So this is the first time we see a time dependent example. So T is just added here in the VARS list. Now as a material we use copper, I've chosen these specific parameters. We have a plane stress example and the thickness specified. And this is the PD setup. Again, it's really simple. You specify the PD components, the load that you'd like. So this is a boundary condition again where it's active. So oh, maybe I should have mentioned. So we're solving actually this problem here, this spoon. So we're holding it fixed here at the back and we're putting a force on top of the spoon and see how it's damped when we let it go. Right. So this is where the force is applied. So we have -20 Newtons in the negative Y direction and we hold it fixed at the other end. We specify initial conditions and the Derivative of initial conditions. It's a wave equation type of PDE. Right. So we solve this with ndsolve. Takes about 60 seconds on my laptop. Again, nothing unusually. The only thing I did here is that I've added the evaluation monitor that I can keep track of when it's time integrates. And here is the displacement in the Y direction. You see, it's not entirely regular. That's because the spoon deforms a bit in itself. Right, okay, so. And here's a visualization of that still. So now you see how the spoon is actually going up and down a bit. Yeah. Okay, so that's that part as a summary. So we spoken a bit about the PD models overview and again, if you are into PDE modeling and would like to get started, the PD models overview page is the place where you should start. We've talked a bit about region symmetry that we've added. We're almost done with this. So 13.2 has a substantial new functionality there. And 13.3 will finish the region symmetry additions that we. We have for 13.1. We've added hyper elastic solid mechanics material models and Rayleigh damping. So these are two new functionalities that we've added. At the end of this talk, I'd like to take the opportunity to thank a few of our customers. Many people helped me getting this all together and I really appreciate that. I'm typically very, how to say, talkative. So if you work with PDES and like to chat with me about that, you're most welcome to send me an email and we'll see. Maybe something like collaboration works comes out of this. So I think I actually do have a little bit of time for questions and I think I would like to take one or two.

Great talk. So I do research in biophysics and computational biology and very big domain of PDEs that's often used in computational biology is reaction diffusion equations. So like it would be. I, I don't see it at the moment in the PD like the link that you shared, but it would be nice in future for researchers like me if, if that class is incorporated and it, it, it is often related to like things like Turing patterns and you can see patterns evolve in space and time. So it's, it's really like, interesting. So yeah, that's just kind of a suggestion.

Right, question. So I think, I mean reaction diffusion equations currently are available, so we have them and we just need to show a bit where this is, where this is to be found in the documentation. So I mean, what you can do is you can look at the mass transport section that will have some reaction diffusion equations or you can look at the PDE terms, the basic building blocks. There is going to be some, some examples for reaction diffusion equations. Again, maybe what, what I'll take from this question is what a good thing would be is to actually add a specific section, the PD models overview page for reaction diffusion type equations. Yeah, we, we can do them and they are there. You're right. Then they're probably not right there, there, there in the way that you immediately see them. But I can add that. I'll. Yeah, I'll. I'll make a note and maybe I can still get that info 13 too.

Okay, thank you so much.

Yeah, sure.

I was wondering, have you. Are you familiar with the package that helps with Maxwell's equations? I'm specifically been interested in an antenna which is a oscillating EM field with. In doing near and far field analysis. And I started and I didn't get very far. I didn't spend much time on it.

So the package you mentioned, is that something. I mean, do you know the name of the package?

No, I was hoping for a source of whether you've had any work in that area or maybe reference to somebody that has.

Okay, so let's start by what the state of our electromagnetic things are. So we have. Maybe I can. Where is my zoom? Disappeared. So we have a section on electromagnetic models in the PD models overview page that has I think three application examples. We don't yet have a theory and background information for these PDE models and they're probably not going to be there super soon. So the next field of physics I'm going to tackle probably is fluid dynamics and after that we'll tackle electromagnetics. And the reason for that is it's a bit technical. The finite element method like we have it now works well for fluid dynamics and probably also electromagnetics. But in electromagnetics often so called vector elements are used and I haven't implemented those yet. So I do the fluid dynamics and then I do the electromagnetics things. That doesn't mean that you cannot do electromagnetics now. The only thing is you'd have to type the equations yourself. There is just no way to generate these equations quite yet. So you can do them, but it's not as automated as for the other fields of physics yet. And let me think, I don't know of a finite element package for electromagnetics. So no, I don't have a name for you there for a reference except the things that we have so far? Do we have more questions? Okay, so no immediate questions. Great. Thanks a lot.