# Ploting the magnitude of a vectorField

**URL:** <https://fenicsproject.discourse.group/t/ploting-the-magnitude-of-a-vectorfield/12363>\
**Category:** General\
**Created:** [October 3, 2023, 9:23am UTC](https://fenicsproject.discourse.group/t/ploting-the-magnitude-of-a-vectorfield/12363 "2023-10-03T09:23:47Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![greg\_perez](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/greg_perez/32/5617_2.png) [@greg\_perez](https://fenicsproject.discourse.group/u/greg_perez)\
**Post date:** [October 3, 2023, 9:23am UTC](https://fenicsproject.discourse.group/t/ploting-the-magnitude-of-a-vectorfield/12363/1 "2023-10-03T09:23:47Z")

</div>

Hello everyone, Im currently trying to solve the Navier Stokes equations, however whenever I plot the velocity results, I get a vector field (as expected). I was wondering if there is any way to plot the magnitude of the velocities instead of the vectorfield itsefl?

My code:

from dolfin import \*

#Creación de un mallado. Se va a refinar en los bordes inferiores y superiores  
n = 50  
mesh = UnitSquareMesh(n,n)

#La malla puede refinarse si es necesario

#Define Taylor-Hood elements (For now, they work best for Navier-Stokes problems)  
V = VectorFunctionSpace(mesh, ‘P’, 2) #Function space that solves for velocity  
Q = FunctionSpace(mesh, ‘P’, 1) #Function space that solves for pressure

#Define Two sets of Test/Trial functions, as we are working with two unknowns.  
u = TrialFunction(V)  
v = TestFunction(V)  
p = TrialFunction(Q)  
q = TestFunction(Q)

def up\_boundary(x, on\_boundary):  
tol = 1E-8  
return(on\_boundary and abs(x[1]-1)\< tol)

def down\_boundary(x, on\_boundary):  
tol = 1E-8  
return(on\_boundary and abs(x[1])\< tol)

def left\_boundary(x, on\_boundary):  
tol = 1E-8  
return(on\_boundary and abs(x[0])\< tol)

def right\_boundary(x, on\_boundary):  
tol = 1E-8  
return(on\_boundary and abs(x[0]-1)\< tol)

DBup = DirichletBC(V, Constant((0,0)) , up\_boundary) #Note that constants must be two dimentional  
DBdown = DirichletBC(V, Constant((0,0)) , down\_boundary)  
DBleft = DirichletBC(V, Constant((1,0)) , left\_boundary)  
DBright = DirichletBC(V, Constant((0,0)) , right\_boundary)

bcs = [DBup, DBdown, DBleft, DBright] #Add all Dirichlet BC into a list

# Define functions for solutions at previous and current time steps

u\_n = Function(V) #Es importante definirlas puesto que voy a iterar en el tiempo  
u\_ = Function(V)  
p\_n = Function(Q)  
p\_ = Function(Q)

#Formulation of the variational form of NavierStokes

#Parameter Creation

T = 1 # final time  
num\_steps = 100 # number of time steps  
dt = T / num\_steps # time step size  
mu = 10 # dynamic viscosity  
rho = 1 # density

# Define expressions used in variational forms

U = 0.5\*(u\_n + u)  
n = FacetNormal(mesh) #Esto me crea un vector normal a la superficie planteada  
f = Constant((0,0)) #Fuerza constante, la idea es ir cambiando este parametro a ver que ocurre  
k = Constant(dt)  
mu = Constant(mu)  
rho = Constant(rho) #Integers must be changed into constant type

# Define symmetric gradient

def epsilon(u):  
return sym(nabla\_grad(u))

# Define stress tensor

def sigma(u, p):  
return 2_mu_epsilon(u) - p\*Identity(len(u))

#Variational form due to Chorin’s Projection.

#Velocity without pressure calculation  
F1 = rho\*dot((u - u\_n) / k, v)_dx + rho_dot(dot(u\_n, nabla\_grad(u\_n)), v)\*dx + inner(sigma(U, p\_n), epsilon(v))_dx + dot(p\_n_n, v)_ds - dot(mu_nabla\_grad(U)\*n, v)\*ds - dot(f, v)\*dx  
a1 = lhs(F1)  
L1 = rhs(F1)

#Pressure Correction  
a2 = dot(nabla\_grad(p), nabla\_grad(q))\*dx  
L2 = dot(nabla\_grad(p\_n), nabla\_grad(q))\*dx - (1/k)\*div(u\_)_q_dx

#Velocity correction due to pressure  
a3 = dot(u, v)\*dx  
L3 = dot(u\_, v)_dx - k_dot(nabla\_grad(p\_ - p\_n), v)\*dx

# Assemble matrices

A1 = assemble(a1)  
A2 = assemble(a2)  
A3 = assemble(a3)

[bc.apply(A1) for bc in bcs] #Apply boundary conditions to A

# Create XDMF files for visualization output

xdmffile\_u = XDMFFile(‘navier\_stokes\_cylinder/velocity.xdmf’)  
xdmffile\_p = XDMFFile(‘navier\_stokes\_cylinder/pressure.xdmf’)

# Save mesh to file (for use in reaction\_system.py)

File(‘navier\_stokes\_cylinder/cylinder.xml.gz’) \<\< mesh

# Time-stepping

import matplotlib.pyplot as plt

t = 0  
for n in range(num\_steps):

```
# Update current time
t += dt

# Step 1: Tentative velocity step
b1 = assemble(L1)
[bc.apply(b1) for bc in bcs]
solve(A1, u_.vector(), b1 ,"gmres", "ilu",) #Preconditioned settings for linear solving methods that accelerate process

# Step 2: Pressure correction step
b2 = assemble(L2)
#[bc.apply(b2) for bc in bc2]
solve(A2, p_.vector(), b2, "gmres", "amg",)

# Step 3: Velocity correction step
b3 = assemble(L3)
solve(A3, u_.vector(), b3, "gmres", "ilu",)

a = plot(u_ , title = "Trial simulation")
plt.colorbar(a)
#plt.draw()
plt.savefig("Instance" + str(round(t,2)) + ".jpg") #save as jpg
#plt.pause(0.02)
plt.clf()

# Update previous solution
u_n.assign(u_)
p_n.assign(p_)

# Update progress bar

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [October 3, 2023, 9:48am UTC](https://fenicsproject.discourse.group/t/ploting-the-magnitude-of-a-vectorfield/12363/2 "2023-10-03T09:48:31Z")

</div>

First of all, you should use 3x` encapsulation to make sure your code is rendered properly.

Secondly, I would advice usage of paraview for such plots.

---

<div class="post-metadata">

**Author:** ![greg\_perez](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/greg_perez/32/5617_2.png) [@greg\_perez](https://fenicsproject.discourse.group/u/greg_perez)\
**Post date:** [October 4, 2023, 1:07pm UTC](https://fenicsproject.discourse.group/t/ploting-the-magnitude-of-a-vectorfield/12363/3 "2023-10-04T13:07:27Z")

</div>

> [@dokken](#):
>
> 3x` encapsulation

Thanks a lot, that is exactly what I did