You're very welcome 😊 I also enjoyed creating these videos a lot. Yes, the goal is to revive the fenics series in the future with some more advanced concepts and then also switch to fenicsx
You are great. Keep going. It will be very interesting to 1) show how to solve a coupled physics problem, 2) make videos discussing the anatomy of the solvers (e.g., how fe.solve works internally!). Thanks
Thanks a lot for the kind and motivating words :). Sure, there will be more advanced topics like this once the series continues with FEniCSx in the future.
Hi, I think you have a mistake in your 3rd equation. The velocity correction requires to add the split term of del pressure over density but it seems you have missed out on the density at the bottom. The equation should read: \frac{u^(t+1) - u^*}{\Delta t} = - \frac{ abla P}{ ho }. Not a major error at all since it's a constant in the projection method.
Hi, yes, we could add the 1/density term there. However, even in the case if density != 1, it would not be necessary as long the density is constant over the entire domain (which is the assumption for incompressible fluid flow). To explain why, take the equation for incompressibility projection and denote by u* the tentative velocity and by u** the velocity at the next time step. Then, we have (1) (u** - u*) / dt = - 1/rho grad(p) Taking the divergence on both sides gives - div(u*) / dt = - 1/ rho laplace(p) Rearraing gives the pressure Poisson problem laplace(p) = rho / dt div(u*) Let's define the inverse Laplace operator (which is linear) p = rho / dt * inv_laplace(div(u*)) Then, we can plug this into (1) to get (u** - u*) / dt = - 1 / rho * rho / dt * grad(inv_laplace(div(u*))) Rearranging for u** let's us cancel all the dt and rho terms, and we get u** = u* - grad(inv_laplace(div(u*))) Sorry, that this was not mentioned in the video. Hope that clears it up :).
Hi, is my understanding of the double contractor operator for (del u) : (del v) correct? del u = [du_x/dx, du_x/dy; du_y/dx, du_y/dy] and del v = [dv_x/dx, dv_x/dy; dv_y/dx, dv_y/dy] then the double contractor is: [du_x/dx * dv_x/dx + du_x/dy * dv_x/dy; du_y/dx * dv_y/dx + du_y/dy * dv_y/dy]. In otherwords the del operator on the vector creates a jacobian matrix which then double contracted becomes a single column vector in this case? And I do element wise matrix multiplication then I sum the rows (not columns) to create a new vector?
Hi, it's been some time since I uploaded the video. Do you have a timestamp to the part in the video you are referring to? Generally speaking, a double contraction of two matrices (as in your case the Jacobian in the form of the velocity gradient) reduces these two matrices to a scalar. So in the case of 2x2 matrices: [a b; c d] : [e f; g h] = a*e + b*f + c*g + d *h Another remark to what you wrote: there are two conventions to layout Jacobian matrices (what the velocity gradients of "del u" and "del v" are): en.wikipedia.org/wiki/Matrix_calculus#Layout_conventions . You can find a discussion on how they are treated in FEniCs on page 58 of the original "FEniCs" tutorial book: link.springer.com/book/10.1007/978-3-319-52462-7
Thank you for such a wonderful and informative video on this interesting topic. I have some questions regarding the pressure poisson equation (PPE). We see that the RHS of the PPE is coming from the tentative velocity-field which is not physical. This will affect the pressure or not? The second question is which type of boundary conditions are suitable for the PPE? Thank you again and waiting for more advance lecture on this topic again like the high order time discretization scemes.
Hi, thanks for the kind comment and the nice feedback 😊 Your two questions are linked and have been for debate in research papers for a long time. Check for instance this paper by Gresho: onlinelibrary.wiley.com/doi/abs/10.1002/fld.1650110509 (if you use connected papers you will find a series of related articles). Simple answers to your questions: Regarding the first: yes, there is an effect to it. Due to the operator splitting we introduce an O(delta t) error to the solution and in the way presented in the video, there is no cycling back and forth between momentum equation and pressure correction which would be necessary in order to obtain a velocity field that both satisfies the (discrete) momentum and continuity equation. Regarding the second question: That's also a debated point in the articles. In the simplest case of wall boundary conditions (i.e. the fluid cannot leave the domain), it is reasonable to have homogeneous Neumann boundary conditions on the pressure, because then there will be zero pressure gradient pointing outward of the domain, and, hence, the corrected velocities respect the walls.
Thank you for this, great content. A question : if you sum step 1 and 3 of the strong form, you see that the approximation to the time derivative of the velocity (u^(t+1) -u^t)/ Delta t is equal almost what we want (ie sum of convection term + pressure gradient + diffusion term) except that the diffusion term is calculated from u^* instead of u^t or u^(t+1). What are the consequences of that? Are there ways to also correct that diffusion term?
Hi, thanks for the comment with the kind words and the interesting question. Indeed, using operator splitting (and in particular projection) methods introduces quite some headaches. Researchers devoted many papers to their analysis, in particular for Boundary Conditions of the intermediate "artificial" PDEs, see for example the publications by Gresho in the early 90s. I must admit my understanding of projection methods is not too deep yet. One consequence of the problem you mention is that the full solution process to the Navier-Stokes equations will always be first order in time. Here, this is not an issue because we are using a first order integrator anyways (an Euler step). If you were to use a second-order method like Crank-Nicolson you would observe a first order convergence despite a second-order consistency. That's really interesting. Regarding your second question: Are there ways to correct it? I would argue that if you are fine with first order convergence in time, there is nothing to worry about. One possible remedy is to use a so-called coupled solver. In a sense, such a solver is not splitting the problem into three sub-PDEs, but rather solves the incompressible Navier-Stokes as one system. If you did the FE discretization, you would get a large sparse system matrix that is in terms of both velocity and pressure. On a first glimpse, that seems intriguing. However, this system matrix is highly ill-condition (due to a so-called Saddle Point Structure). There are a lot of mathematical resources on how to still solve such systems, but one has to dive quite deep. Bottom-line, I would say it's not too worrying that the diffusion is treated at the intermediate step in time, as long as you are fine with first order convergence. If you need better, you probably have to look into more advanced methods.
Thank you for your helpful video. I have a question, Some people or papers write the weak form (or variational form?) of the viscous term as (∇u,∇v) (inner product?). But the other person like you describes the viscous term as a ∇u:∇v. Who's right? What's different?
That's a great question 👍 Both approaches yield the same result. Let me explain: The ":" (colon) operator here refers to double-contraction. If you have only dot between the two vectors, this typically means the dot-product (so a euclidean (l2/p2) inner product), in which you elementwise multiply the two vectors and sum up the elements afterwards. With the ":" you would do sth similar for two (identically sized) matrices: elementwise multiply the entries and sum them up (resulting in one scalar value). Since u is a vector, it's gradient (the velocity gradient) is a matrix. As such we need the colon. However, we could also use the notation involving the inner product between two matrices, since this one then usually refers to the generalization of the euclidean (l2/p2) inner product, both approaches are identical. So both are a way of contracting two matrices down into a scalar value. The advantage of using inner(a,b) is that it works no matter if a&b are scalars, vectors, matrices etc. because by axiom an inner product must always map to a scalar value. (In the corresponding FEniCs videos I also use it that way). Hope that helps 😊
@@MachineLearningSimulation I'm not good at math, so this is the first time I've realized that the generalized (,) operator maps to a scalar value without regard to the target's rank. Thank you so much for your answer. It's very valuable.
Hi, thanks for the question! :) There are built-in ways to export simulation results into a format that Paraview can understand. For instance, you can check this video of mine: ruclips.net/video/ibgALqJOX-I/видео.html In general, the meshes in FEniCs are unstructured and due to the FEM representation, the dof values do not necessarily represent the values on mesh nodes (but only coefficients for basis functions). [This can coincide for simple basis functions, like linear functions]. What I would recommend to you is to create a uniform mesh over the domain you are interested in and then use the FEniCs functionality to interpolate the solution function at these coordinates. You probably have to check the documentation for this. Best of luck :).
That refers to vectors and matrices. One underscore is a vector. Two is a matrix. More refer to higher tensorical quantities. Such a notation is often used in handwritten notes for mechanics, but is seldomly used in machine learning. Personally, I like it a lot since it gives more insight into the linear algebra operations.
Thank you very much for the videos on Fluid dynamics. And please consider switching to FEniCSx for your future demonstrations...
You're very welcome 😊 I also enjoyed creating these videos a lot.
Yes, the goal is to revive the fenics series in the future with some more advanced concepts and then also switch to fenicsx
Sounds great!!😊😊
I see u have put out a lot of great FEM content. Sadly due to exams I only watched the NS videos , and I have to say excellent work!
Thanks for the kind words 😊
Hope your exams went well and enjoy the videos if you want to catch up.
@@MachineLearningSimulation I will , sadly I have 2 and a half more weeks. But when I do finish I plan to watch them. Keep it up!
You are great. Keep going. It will be very interesting to 1) show how to solve a coupled physics problem, 2) make videos discussing the anatomy of the solvers (e.g., how fe.solve works internally!). Thanks
Thanks a lot for the kind and motivating words :).
Sure, there will be more advanced topics like this once the series continues with FEniCSx in the future.
glad i found your channel. great explanation👍
Welcome aboard! 😊
You're very welcome 🤗
Hi, I think you have a mistake in your 3rd equation. The velocity correction requires to add the split term of del pressure over density but it seems you have missed out on the density at the bottom. The equation should read: \frac{u^(t+1) - u^*}{\Delta t} = - \frac{
abla P}{
ho }. Not a major error at all since it's a constant in the projection method.
Hi,
yes, we could add the 1/density term there. However, even in the case if density != 1, it would not be necessary as long the density is constant over the entire domain (which is the assumption for incompressible fluid flow).
To explain why, take the equation for incompressibility projection and denote by u* the tentative velocity and by u** the velocity at the next time step.
Then, we have
(1) (u** - u*) / dt = - 1/rho grad(p)
Taking the divergence on both sides gives
- div(u*) / dt = - 1/ rho laplace(p)
Rearraing gives the pressure Poisson problem
laplace(p) = rho / dt div(u*)
Let's define the inverse Laplace operator (which is linear)
p = rho / dt * inv_laplace(div(u*))
Then, we can plug this into (1) to get
(u** - u*) / dt = - 1 / rho * rho / dt * grad(inv_laplace(div(u*)))
Rearranging for u** let's us cancel all the dt and rho terms, and we get
u** = u* - grad(inv_laplace(div(u*)))
Sorry, that this was not mentioned in the video. Hope that clears it up :).
Well done. Thank you very much!
You're welcome 🤗
Thanks for the kind feedback
Thank you! Great video!
You're very welcome :) Thanks a lot.
Hi, is my understanding of the double contractor operator for (del u) : (del v) correct? del u = [du_x/dx, du_x/dy; du_y/dx, du_y/dy] and del v = [dv_x/dx, dv_x/dy; dv_y/dx, dv_y/dy] then the double contractor is: [du_x/dx * dv_x/dx + du_x/dy * dv_x/dy; du_y/dx * dv_y/dx + du_y/dy * dv_y/dy]. In otherwords the del operator on the vector creates a jacobian matrix which then double contracted becomes a single column vector in this case? And I do element wise matrix multiplication then I sum the rows (not columns) to create a new vector?
Hi,
it's been some time since I uploaded the video. Do you have a timestamp to the part in the video you are referring to?
Generally speaking, a double contraction of two matrices (as in your case the Jacobian in the form of the velocity gradient) reduces these two matrices to a scalar.
So in the case of 2x2 matrices:
[a b; c d] : [e f; g h] = a*e + b*f + c*g + d *h
Another remark to what you wrote: there are two conventions to layout Jacobian matrices (what the velocity gradients of "del u" and "del v" are): en.wikipedia.org/wiki/Matrix_calculus#Layout_conventions . You can find a discussion on how they are treated in FEniCs on page 58 of the original "FEniCs" tutorial book: link.springer.com/book/10.1007/978-3-319-52462-7
Thank you very much!
You're welcome! 🤗
Thank you for such a wonderful and informative video on this interesting topic.
I have some questions regarding the pressure poisson equation (PPE). We see that the RHS of the PPE is coming from the tentative velocity-field which is not physical. This will affect the pressure or not? The second question is which type of boundary conditions are suitable for the PPE?
Thank you again and waiting for more advance lecture on this topic again like the high order time discretization scemes.
Hi, thanks for the kind comment and the nice feedback 😊
Your two questions are linked and have been for debate in research papers for a long time. Check for instance this paper by Gresho: onlinelibrary.wiley.com/doi/abs/10.1002/fld.1650110509 (if you use connected papers you will find a series of related articles).
Simple answers to your questions:
Regarding the first: yes, there is an effect to it. Due to the operator splitting we introduce an O(delta t) error to the solution and in the way presented in the video, there is no cycling back and forth between momentum equation and pressure correction which would be necessary in order to obtain a velocity field that both satisfies the (discrete) momentum and continuity equation.
Regarding the second question: That's also a debated point in the articles. In the simplest case of wall boundary conditions (i.e. the fluid cannot leave the domain), it is reasonable to have homogeneous Neumann boundary conditions on the pressure, because then there will be zero pressure gradient pointing outward of the domain, and, hence, the corrected velocities respect the walls.
Thank you for this, great content. A question : if you sum step 1 and 3 of the strong form, you see that the approximation to the time derivative of the velocity (u^(t+1) -u^t)/ Delta t is equal almost what we want (ie sum of convection term + pressure gradient + diffusion term) except that the diffusion term is calculated from u^* instead of u^t or u^(t+1). What are the consequences of that? Are there ways to also correct that diffusion term?
Hi, thanks for the comment with the kind words and the interesting question. Indeed, using operator splitting (and in particular projection) methods introduces quite some headaches. Researchers devoted many papers to their analysis, in particular for Boundary Conditions of the intermediate "artificial" PDEs, see for example the publications by Gresho in the early 90s.
I must admit my understanding of projection methods is not too deep yet. One consequence of the problem you mention is that the full solution process to the Navier-Stokes equations will always be first order in time. Here, this is not an issue because we are using a first order integrator anyways (an Euler step). If you were to use a second-order method like Crank-Nicolson you would observe a first order convergence despite a second-order consistency. That's really interesting.
Regarding your second question: Are there ways to correct it? I would argue that if you are fine with first order convergence in time, there is nothing to worry about. One possible remedy is to use a so-called coupled solver. In a sense, such a solver is not splitting the problem into three sub-PDEs, but rather solves the incompressible Navier-Stokes as one system. If you did the FE discretization, you would get a large sparse system matrix that is in terms of both velocity and pressure. On a first glimpse, that seems intriguing. However, this system matrix is highly ill-condition (due to a so-called Saddle Point Structure). There are a lot of mathematical resources on how to still solve such systems, but one has to dive quite deep.
Bottom-line, I would say it's not too worrying that the diffusion is treated at the intermediate step in time, as long as you are fine with first order convergence. If you need better, you probably have to look into more advanced methods.
Thank you for your helpful video. I have a question,
Some people or papers write the weak form (or variational form?) of the viscous term as (∇u,∇v) (inner product?).
But the other person like you describes the viscous term as a ∇u:∇v.
Who's right? What's different?
That's a great question 👍
Both approaches yield the same result. Let me explain:
The ":" (colon) operator here refers to double-contraction. If you have only dot between the two vectors, this typically means the dot-product (so a euclidean (l2/p2) inner product), in which you elementwise multiply the two vectors and sum up the elements afterwards. With the ":" you would do sth similar for two (identically sized) matrices: elementwise multiply the entries and sum them up (resulting in one scalar value). Since u is a vector, it's gradient (the velocity gradient) is a matrix. As such we need the colon.
However, we could also use the notation involving the inner product between two matrices, since this one then usually refers to the generalization of the euclidean (l2/p2) inner product, both approaches are identical.
So both are a way of contracting two matrices down into a scalar value.
The advantage of using inner(a,b) is that it works no matter if a&b are scalars, vectors, matrices etc. because by axiom an inner product must always map to a scalar value. (In the corresponding FEniCs videos I also use it that way).
Hope that helps 😊
@@MachineLearningSimulation I'm not good at math, so this is the first time I've realized that the generalized (,) operator maps to a scalar value without regard to the target's rank. Thank you so much for your answer. It's very valuable.
Thanks 🙏
I'm very happy, I could help.
Good luck with Your further learning journey 😊
you are great
Thanks 🙏
How can we save the simulation results in a way that I could read using numpy to train some machine model on them (like a pinn)?
Hi, thanks for the question! :)
There are built-in ways to export simulation results into a format that Paraview can understand. For instance, you can check this video of mine: ruclips.net/video/ibgALqJOX-I/видео.html
In general, the meshes in FEniCs are unstructured and due to the FEM representation, the dof values do not necessarily represent the values on mesh nodes (but only coefficients for basis functions). [This can coincide for simple basis functions, like linear functions].
What I would recommend to you is to create a uniform mesh over the domain you are interested in and then use the FEniCs functionality to interpolate the solution function at these coordinates. You probably have to check the documentation for this. Best of luck :).
What do the underlined letters refer to?
That refers to vectors and matrices. One underscore is a vector. Two is a matrix. More refer to higher tensorical quantities.
Such a notation is often used in handwritten notes for mechanics, but is seldomly used in machine learning. Personally, I like it a lot since it gives more insight into the linear algebra operations.
@@MachineLearningSimulation ok, thanks a lot.
@@McHumaty You're welcome :) Thanks for the question.
Bro, you are the 🐑