# Calculating the jacobian of norm and square root terms in the Finite Element Method

In the code that my group is writing (Lethe) we use a stabilized approach to solve the Navier-Stokes equation. The GLS stabilized method we use has a stabilisation term which contains a stabilization parameter such as: $$\tau =\left[\left( \frac{2\rho |\mathbf{u}|}{h} \right)^2+\left( \frac{4 \mu}{\gamma h^2} \right) \right]^{0.5}$$

where $$\gamma$$ is a constant that depends on the type of element, $$\rho$$ is the density, $$h$$ is the element size and $$\mu$$ is the dynamic viscosity.

This stabilization parameter does not depend on the flow variable except for the first term which contains the absolute value of the velocity. I find myself unable to calculate the jacobian of that term, so right now we assemble an analytical jacobian for everything but this term.

Is there anyway to proceed to establish the jacobian of this term analytically? I understand that the $$u^2$$ term can be treated straightforwardly, but I am unsure on how to calculate the Frechet derivative of the square root.

• The term is squared. You could just think of it as having a $u^2$ term instead, which is amenable to analytical derivatives.
– Paul
Nov 16 '19 at 23:27
• Yeah that make sense. I am more confused about the square root part. I am under the impression that it would lead to an atrocious number of terms :(
– BlaB
Nov 16 '19 at 23:33
• @Paul do you mean $\tau^2$? or I totally misunderstood both the question and your comment? Nov 17 '19 at 19:46
• @AntonMenshov what I am actually wondering is how to calculate the jacobian of $\tau$ with respect to $\mathbf{u}$. So what you are suggesting is that I derive $\tau^2$ with respect to $\mathbf{u}$, which is actually very doable, and use that to establish my jacobian by dividing the result by $\tau$?
– BlaB
Nov 17 '19 at 22:01
• A good sanity check for any solution you do come up with is that the quantity you describe is convex with respect to $u$. Nov 21 '19 at 18:44

You can arrive at the Jacobian analytically it just takes a few steps

So assuming we have our typical FE field values:

$$u_i = \sum_j \phi^j u_i^j$$

Where $$i$$ represents coordinate direction, $$\phi$$ our shape function, and $$j$$ index to DOF.

$$\tau = \sqrt{\alpha \lVert u \rVert + \beta}$$

Define some functions for ease of writing derivatives

$$g(\mathbf{u}) = \alpha f(\mathbf{u}) + \beta$$ $$f(\mathbf{u}) = \lVert u \rVert = \sqrt{u_x^2 + u_y^2 + u_z^2}$$

So now stabilization parameter can be written as

$$\tau = \sqrt{g(\mathbf{u})}$$

And our derivatives will be in the form

$$\frac{\partial\tau}{\partial u_i} = \frac{\partial g / \partial u_i}{2 \tau}$$

And so we find derivative of $$g$$

$$\frac{\partial g}{ \partial u_i} = \alpha \frac{\partial f}{\partial u_i}$$

and finally derivative of $$f$$

$$\frac{\partial f}{ \partial u_i} = \frac{u_i \phi^j}{f(\mathbf{u})}$$

Putting it all together:

$$\frac{\partial\tau}{\partial u_i} = \frac{\alpha u_i \phi^j}{2 \tau \lVert u \rVert}$$

You may want to double-check my math but it should give you the idea of how to do it.

• Explained this way it is perfectly clear. Thanks :)!
– BlaB
Nov 22 '19 at 7:44