Hello everyone! I try to study the sensitivity to base-flow modifications on the solutions of Navier-Stokes equations.
The expression is
\begin{equation}
\nabla \sigma = -(grad(u))^{H}\cdot u^ (grad(u^))\cdot u^{*}
\end{equation}
where u^ is the adjoint solution.
I need to calculate adjoint(grad(u)), and I believe I have to convert grad(u) to numpy array, but I don’t know is it the best way.
Thank you
PD: How can I write latex in the forum?
To write latex, use $-encapsulation, for instance
$\frac{a}{b}$ yields
\frac{a}{b}.
To convert grad(u) to a numpy array you need to project it to an appropriate function space, and evaluate it at the coordinates of interest. Anorher option is to get the values at the degrees of freedom, procjeted_func.vector().get_local() But this would Give you a 1D array. Using appropriate subspaces along with this function, you can build the matrix.
Note that there are several implementations of adjoint based sensitivity methods for dolfin, such as dolfin-adjoint. This add-on to dolfin can automatically compute sensivities with respect to for instance boundary conditions or the shape of the domain.
Mmmm. procjeted_func.vector().get_local() doesn’t work for me. If I try these and works, I can do
procjeted_func.vector().get_local()=adjoint(procjeted_func.vector().get_local()) and
dot(projected_fun,u)=(\nabla u)^{\dagger}\cdot u if project_fun = \nabla u?
This is the code :
#Solution of LNS, velocity and pressure
P2 = VectorElement(“P”, mesh.ufl_cell(), 2)
P1 = FiniteElement(“P”, mesh.ufl_cell(), 1)
TH = MixedElement([P2, P1])
W = FunctionSpace(mesh, TH)
wauxa=Function(W)
waux=Function(W)
r,c,rx,cx = eigensolver.get_eigenpair(indexea)
ra,ca,rax,cax = Aeigensolver.get_eigenpair(indexea)
waux.vector()[:]=rx
wauxa.vector()[:]=rax
uaux,paux= waux.split(deepcopy=True)
uauxa,pauxa= wauxa.split(deepcopy=True)
Q = FunctionSpace(mesh, ‘P’, 2)
V = VectorFunctionSpace(mesh, ‘P’, 2)
uv=project(uaux,V)
auv=project(uauxa,V)
T = TensorFunctionSpace(mesh, “P”, 2)
grad_u=project(grad(uv),T)
auxiliary=grad_u.vector().get_local()
and I have another question, but I think it would be appropriate to put it in another post.
It’s about how obtain the norm of Vector Function in each element.
This is my way:
uvx,uvy=split(uv)
auvx,auvy=split(auv)
p_uvx=project(uvx,Q)
p_uvy=project(uvy,Q)
np_uvx=p_uvx.vector()[:]
np_uvy=p_uvy.vector()[:]
p_auvx=project(auvx,Q)
p_auvy=project(auvy,Q)
np_auvx=p_auvx.vector()[:]
np_auvy=p_auvy.vector()[:]
np_norm=np_uvxnp_uvx + np_uvynp_uvy
np_norm=np.sqrt(np.array(np_norm))
np_anorm=np_auvxnp_auvx + np_auvynp_auvy
np_anorm=np.sqrt(np.array(np_anorm))
norm_uv=Function(Q)
norm_uv.vector()[:]=np_norm
norm_auv=Function(Q)
norm_auv.vector()[:]=np_anorm
I think this works, but there is another quicker way?
Thank you!
Please specify what you mean by it doesn’t work for you? Does it throw an error?
Also please format your code using 3x` encapsulation, such that it is formatted properly.
Finally make the problem a minimal working code example, as described in:
https://fenicsproject.discourse.group/t/read-before-posting-how-do-i-get-my-question-answered/21/3
For your particular example, the following is a minimal working code example:
from dolfin import *
mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, "CG", 1)
u = project(Expression(("x[0]","x[1]"),degree=1), V)
# Note that the gradient is in a DG space of one degree lower than the original space
T = TensorFunctionSpace(mesh, "DG", 0)
gradu = project(grad(u), T)
print(gradu.vector().get_local())
Sorry… I had another error in my code and thought that was it. These work!
I think that this code should wokrs
grad_u=-(skew(grad_u)-sym(grad_u))
Transpose doesn’t work: Error name ‘transpose’ is not defined, but here it is
Again, you need to supply a minimal code example, and not just one line of a longer code. Please also spend some time investigating the cause of the error.
Here is my minimal working version:
from dolfin import *
mesh = UnitSquareMesh(10, 10)
V = VectorFunctionSpace(mesh, "CG", 1)
u = project(Expression(("x[0]","x[1]"),degree=1), V)
T = TensorFunctionSpace(mesh, "DG", 0)
gradu = project(grad(u), T)
print(gradu.vector().get_local())
gradu=-(skew(gradu)-sym(gradu))
Dokken, thank you so much for your time! Your help is greatly appreciated. I had a issue with Jupiter lab.
All works perfectly! I’m sorry how I wrote the messages. The next post will be better.
Hi! I was wondering how we do this in dolfinx, since there is no project function anymore.
See my comment at: How to compute the strain & stress matrices? - #7 by dokken
You can also use DOLFINx Expression for this purpose.
In particular, use: dolfinx/function.py at 9548c27e45bd0ca78d351964820cc4a4a7c0957d · FEniCS/dolfinx · GitHub
dolfinx/test_expression.py at 9548c27e45bd0ca78d351964820cc4a4a7c0957d · FEniCS/dolfinx · GitHub
You could also have a look at
https://jsdokken.com/fenics22-tutorial/heat_eq.html#post-processing-without-projections
Thank you so much for your prompt reply!