Compute grad(u) on time dependent unknown u problem


163
views
0
4 months ago by
Dear FEniCS community, currently I am working on time dependent free convective flow.

I compute the solution for u (velocity) and T (temperature) using splitting scheme.

How can I get the grad(u) solution from the u? I need the solution for grad(u) and plotting the graph.

Thank you for your help.
Community: FEniCS Project
1
The velocity gradient is a second order tensor, how do you want to plot it?
Anyway, you can get a field for post-processing by projecting the gradient onto a suitable function space:
gradu = project(grad(u), FunctionSpace(u.ufl_function_space().mesh(), TensorElement(u.ufl_function_space().ufl_element().sub_elements()[0])))​

And write it to a file that can be read with paraview:
gradFile = File('gradu.pvd')
gradFile << gradu​

You can then view all nine components as single scalar fields in paraview.
written 4 months ago by klunkean  
from __future__ import print_function
from fenics import *

#define parameter
Tf = 1
Re = 100
Pr = 0.71
Gr = 1000000
M = 500
dt = 0.1
num_steps = 100

#create mesh
mesh = RectangleMesh(Point(0,0),Point(1,1),128,128)

#define function space
V = VectorFunctionSpace(mesh,'P', 2)
Q = FunctionSpace(mesh,'P', 1)
L = FunctionSpace(mesh,'P', 2)

#define boundaries
walls = 'x[1] < DOLFIN_EPS'
inflow = 'x[0] < DOLFIN_EPS'
outflow = 'x[0] > 1 - DOLFIN_EPS'
upper = 'x[1] > 1 - DOLFIN_EPS'

#define boundary condition

#walls
bcu_walls = DirichletBC(V, Constant((0,0)), walls)
bcT_walls = DirichletBC(L, Constant(1), walls) #T_w

#inflow
bcu_inflow = DirichletBC(V, Constant((1,0)), inflow)
bcT_inflow = DirichletBC(L, Constant(0), inflow)

#upper
bcu_upper= DirichletBC(V, Constant((1,0)), upper)
bcT_upper = DirichletBC(L, Constant(0), upper)

#outflow
#bcp_outflow = DirichletBC(Q, Constant(0), outflow)

bcu = [bcu_walls, bcu_inflow, bcu_upper]
bcp = []
bcT = [bcT_walls, bcT_inflow, bcT_upper]

#define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
T = TrialFunction(L)
s = TestFunction(L)

# Create functions
u1 = Function(V)
u0 = Function(V)
p1 = Function(Q)
p0 = Function(Q)
T1 = Function(L)
T0 = Function(L)

#define expression used in variational forms
Re = Constant (Re)
Gr = Constant (Gr)
M = Constant (M)
Pr = Constant (Pr)
j = Constant((0, 1))
k = Constant (dt)

# Step 1: Tentative velocity step
F1 = (1/k)*dot(u - u0, v)*dx \
+ dot(dot(u0,nabla_grad(u)),v)*dx \
+ pow(Re,-1)*inner(nabla_grad(u),nabla_grad(v))*dx \
+ dot(grad(p0),v)*dx \
- Gr*pow(Re,-2)*T0*dot(j,v)*dx \
+ M*dot(u,v)*dx  #(-) aiding, F=0

a1 = lhs(F1)
L1 = rhs(F1)

# Step 2:Pressure correction
a2 = inner(grad(p), grad(q))*dx
L2 = (1/k)*inner(u1- u0, grad(q))*dx

# Step 3:Pressure update
a3 = p*q*dx
L3 = dot(p1,q)*dx + dot(p0,q)*dx \
- pow(Re,-1)*div(u1)*q*dx 

# Step 4: Temperature correction
F4 = (1/k)*inner(T - T0, s)*dx \
- pow(Re,-1)*pow(Pr,-1)*dot(grad(T),grad(s))*dx \
+ inner(inner(grad(T),u1), s)*dx 

a4 = lhs(F4)
L4 = rhs(F4)

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
A4 = assemble(a4)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
[bc.apply(A3) for bc in bcp]
[bc.apply(A4) for bc in bcT]
 
# Create vtk files for visualization output
vtkfile_u = File('unsteady_free_plateD_0.002_M500/u.pvd')
vtkfile_p = File('unsteady_free_plateD_0.002_M500/p.pvd')
vtkfile_T = File('unsteady_free_plateD_0.002_M500/T.pvd')

# Create progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

# Update current time
    t += dt
    print('current t is',t)

    # Step 1
    begin("Computing tentative velocity")
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u1.vector(), b1,'gmres', 'ilu')
    end()

    # Step 2
    begin("Computing pressure correction")
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p1.vector(), b2,'gmres', 'hypre_amg')
    end()

    # Step 3
    begin("Computing pressure update")
    b3 = assemble(L3)
    [bc.apply(b3) for bc in bcp]
    solve(A3, p1.vector(), b3,'gmres', 'amg')
    end()

    # Step 4
    begin("Computing temperature correction")
    b4 = assemble(L4)
    [bc.apply(b4) for bc in bcT]
    solve(A4, T1.vector(), b4,'gmres', 'amg')
    end()

    # Plot solution
    plot(u1, title='Velocity')
    plot(p1, title='Pressure')
    plot(T1, title='Temperature')

    # Save solution to file (vtk)
    vtkfile_u << (u1, t)
    vtkfile_p << (p1, t)
    vtkfile_T << (T1, t)

    # Update previous solution
    u0.assign(u1)
    p0.assign(p1)
    T0.assign(T1)

    # Update progress bar
    progress.update(t / Tf)
    

# Hold plot
interactive()


This is my coding. I dont know where to write the coding on grad(u). I already tried but failed. Thank you for your help.

written 4 months ago by raihan asimoni  
Hi raihan,

do you have any code that you have tried? is a 1D or 2D problem? do you want to graph some sort of vector field?
written 4 months ago by Ruben Gonzalez  
Yes, below I attach my coding on klunkean comments. I am doing 2D on vertical plate. Could you help me? Thank you.
written 4 months ago by raihan asimoni  

1 Answer


3
4 months ago by

I made three changes: I defined a new FunctionSpace below your FunctionSpace definitions (lines 21-23), I created another pvd file below your files (line 120) and performed the projection which then is written to file (lines 172-173). This Code works for me:

from __future__ import print_function
from fenics import *

#define parameter
Tf = 1
Re = 100
Pr = 0.71
Gr = 1000000
M = 500
dt = 0.1
num_steps = 100

#create mesh
mesh = RectangleMesh(Point(0,0),Point(1,1),128,128)

#define function space
V = VectorFunctionSpace(mesh,'P', 2)
Q = FunctionSpace(mesh,'P', 1)
L = FunctionSpace(mesh,'P', 2)

FE = FiniteElement('DG', mesh.ufl_cell(), 1)
TFE = TensorElement(FE)
TV = FunctionSpace(mesh, TFE)

#define boundaries
walls = 'x[1] < DOLFIN_EPS'
inflow = 'x[0] < DOLFIN_EPS'
outflow = 'x[0] > 1 - DOLFIN_EPS'
upper = 'x[1] > 1 - DOLFIN_EPS'

#define boundary condition

#walls
bcu_walls = DirichletBC(V, Constant((0,0)), walls)
bcT_walls = DirichletBC(L, Constant(1), walls) #T_w

#inflow
bcu_inflow = DirichletBC(V, Constant((1,0)), inflow)
bcT_inflow = DirichletBC(L, Constant(0), inflow)

#upper
bcu_upper= DirichletBC(V, Constant((1,0)), upper)
bcT_upper = DirichletBC(L, Constant(0), upper)

#outflow
#bcp_outflow = DirichletBC(Q, Constant(0), outflow)

bcu = [bcu_walls, bcu_inflow, bcu_upper]
bcp = []
bcT = [bcT_walls, bcT_inflow, bcT_upper]

#define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)
T = TrialFunction(L)
s = TestFunction(L)

# Create functions
u1 = Function(V)
u0 = Function(V)
p1 = Function(Q)
p0 = Function(Q)
T1 = Function(L)
T0 = Function(L)

#define expression used in variational forms
Re = Constant (Re)
Gr = Constant (Gr)
M = Constant (M)
Pr = Constant (Pr)
j = Constant((0, 1))
k = Constant (dt)

# Step 1: Tentative velocity step
F1 = (1/k)*dot(u - u0, v)*dx \
+ dot(dot(u0,nabla_grad(u)),v)*dx \
+ pow(Re,-1)*inner(nabla_grad(u),nabla_grad(v))*dx \
+ dot(grad(p0),v)*dx \
- Gr*pow(Re,-2)*T0*dot(j,v)*dx \
+ M*dot(u,v)*dx  #(-) aiding, F=0

a1 = lhs(F1)
L1 = rhs(F1)

# Step 2:Pressure correction
a2 = inner(grad(p), grad(q))*dx
L2 = (1/k)*inner(u1- u0, grad(q))*dx

# Step 3:Pressure update
a3 = p*q*dx
L3 = dot(p1,q)*dx + dot(p0,q)*dx \
- pow(Re,-1)*div(u1)*q*dx 

# Step 4: Temperature correction
F4 = (1/k)*inner(T - T0, s)*dx \
- pow(Re,-1)*pow(Pr,-1)*dot(grad(T),grad(s))*dx \
+ inner(inner(grad(T),u1), s)*dx 

a4 = lhs(F4)
L4 = rhs(F4)

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)
A4 = assemble(a4)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]
[bc.apply(A3) for bc in bcp]
[bc.apply(A4) for bc in bcT]
 
# Create vtk files for visualization output
vtkfile_u = File('unsteady_free_plateD_0.002_M500/u.pvd')
vtkfile_p = File('unsteady_free_plateD_0.002_M500/p.pvd')
vtkfile_T = File('unsteady_free_plateD_0.002_M500/T.pvd')
vtkfile_gradu = File('unsteady_free_plateD_0.002_M500/gradu.pvd')

# Create progress bar
progress = Progress('Time-stepping')
set_log_level(PROGRESS)

# Time-stepping
t = 0
for n in range(num_steps):

# Update current time
    t += dt
    print('current t is',t)

    # Step 1
    begin("Computing tentative velocity")
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u1.vector(), b1,'gmres', 'ilu')
    end()

    # Step 2
    begin("Computing pressure correction")
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p1.vector(), b2,'gmres', 'hypre_amg')
    end()

    # Step 3
    begin("Computing pressure update")
    b3 = assemble(L3)
    [bc.apply(b3) for bc in bcp]
    solve(A3, p1.vector(), b3,'gmres', 'amg')
    end()

    # Step 4
    begin("Computing temperature correction")
    b4 = assemble(L4)
    [bc.apply(b4) for bc in bcT]
    solve(A4, T1.vector(), b4,'gmres', 'amg')
    end()

    # Plot solution
    plot(u1, title='Velocity')
    plot(p1, title='Pressure')
    plot(T1, title='Temperature')

    # Save solution to file (vtk)
    vtkfile_u << (u1, t)
    vtkfile_p << (p1, t)
    vtkfile_T << (T1, t)
    
    vel_grad = project(nabla_grad(u1), TV)
    vtkfile_gradu << (vel_grad, t)
    
    # Update previous solution
    u0.assign(u1)
    p0.assign(p1)
    T0.assign(T1)

    # Update progress bar
    progress.update(t / Tf)
    

# Hold plot
interactive()
1
Nice!! it works. The only issue is that Fenics assigns a different variable name at each time step, that is why Paraview crashes after the first time step. It is easily fixed by adding this line after the line  171

vel_grad = project(nabla_grad(u1), TV) #line 171
vel_grad.rename('gradu','gradu') #rename to the same name every time step
vtkfile_gradu << (vel_grad, t) #save into file​

Also it is tricky to check all the components in Paraview as it only shows you one variable, in the dropdown menu next to the eye button, in this case it would be 'gradu' but in the next dropdown box you will find all the components of the tensor.

This lead me to these question, what is the meaning of all the components 0-8? in this example most of them are zero but (0 and 3) I thought it was the position in the tensor 0 = xx but 4 = yy then why 3 is nonzero?

Last thing, could you please explain a little bit more about your selection of FiniteElement? why 'DG'?

Thank you very much for your help!!! awesome contribution
written 4 months ago by Ruben Gonzalez  
I open in paraview but failed to run. Do I need to plot the vel_grad?
written 4 months ago by raihan asimoni  
I can open the gradu.pvd file in paraview and plot the components.
What exactly fails to run?
written 4 months ago by klunkean  
I can open, but when I start to run the time, paraview automatically closed. I dont know whats the problem.
written 4 months ago by raihan asimoni  
1
Hi raihan, please read my answer above.
written 4 months ago by Ruben Gonzalez  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »