L2 norm of Heat equation with u as vector function

mathematics
1

Dear @dokken , @nate
Thanks .

when i am trying to find out L2-norm of heat equation on the

V = VectorFunctionSpace(mesh, 'P', 1)

the error decreases and remains constant.
Here is the code

from fenics import *    
def heat(N):    
    T = 2.0            # final time
    num_steps = 4     # number of time steps
    dt = T / num_steps # time step size
    alpha = 3          # parameter alpha
    beta = 1.2         # parameter beta

    # Create mesh and define function space
    mesh = UnitSquareMesh(N, N)
    V = VectorFunctionSpace(mesh, 'P', 1)

    # Define boundary condition
    u_D = Expression(('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t','1 + x[0]*x[1]+beta*pow(t,2)'),
                 degree=2, alpha=alpha, beta=beta, t=0)

    bc = DirichletBC(V, u_D, "on_boundary")

    # Define initial value
    u_n = interpolate(u_D, V)
    #u_n = project(u_D, V)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Expression(('beta - 2 - 2*alpha', '2*beta*t'), degree=1, alpha=alpha, beta=beta, t=0 )

    F = inner(u,v)*dx + dt*inner(grad(u), grad(v))*dx - inner(u_n + dt*f,v)*dx
                                 
    a = lhs(F)
    L = rhs(F)
    
    
    # Time-stepping
    u = Function(V)
    t = 0
    for n in range(num_steps):

        # Update current time
        t += dt
        u_D.t = t
        f.t   = t
        # Compute solution
        solve(a == L, u, bc)

        # Plot solution
        #plot(u)

       
        error_L2 = errornorm(u_D, u, 'L2')
        # Update previous solution
        u_n.assign(u)
                                 

    return error_L2

for N in [2, 4, 8,16, 32, 64]:
    print("{0:d} error_L2 ={1:.2e}".format(N, heat(N)))

2 error_L2 =1.80e-01
4 error_L2 =5.12e-02
8 error_L2 =2.73e-02
16 error_L2 =2.49e-02
32 error_L2 =2.48e-02
64 error_L2 =2.48e-02

          
Kindly clarify this issue.

when I am using `V = VectorFunctionSpace(mesh, 'P', 2)`
error is coming like

2 error_L2 =2.39e-02
4 error_L2 =2.47e-02
8 error_L2 =2.48e-02
16 error_L2 =2.48e-02
32 error_L2 =2.48e-02
64 error_L2 =2.48e-02

kindly, give the necessary  suggestion to find out the order of convergence and mention if there is any issue with this code.

Waiting for a positive reply!

Thank You,
Debendra
2

Let us consider a simplified example, based on:
https://fenicsproject.org/pub/tutorial/html/._ftut1006.html
where I have changed to the manufactured solution to be beta*t*t instead of beta*t (and correspondingly the source term from (beta to 2betat) As we observe both a spatial and temporal discretization error, the key to look at convergence in one of them is by reducing the other:

from fenics import *
import numpy as np

set_log_level(LogLevel.ERROR)    
def heat(N, num_steps):    

    T = 2.0            # final time
    dt = T / num_steps # time step size
    alpha = 3          # parameter alpha
    beta = 1.2         # parameter beta

    # Create mesh and define function space
    mesh = UnitSquareMesh(N, N)
    V = FunctionSpace(mesh, 'P', 1)

    # Define boundary condition
    u_D = Expression('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t*t',
                    degree=2, alpha=alpha, beta=beta, t=0)

    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

    # Define initial value
    u_n = interpolate(u_D, V)
    #u_n = project(u_D, V)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    t = 0

    f = Constant(2*beta*t - 2 - 2*alpha)

    F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_n + dt*f)*v*dx
    a, L = lhs(F), rhs(F)

    # Time-stepping
    u = Function(V)
    error_T = 0

    for n in range(num_steps):

        # Update current time
        t += dt
        u_D.t = t
        f.assign(2*beta*t -2 - 2*alpha)
        # Compute solution
        solve(a == L, u, bc)


        # Compute error at vertices
        error = errornorm(u_D, u)
        error_T += dt*error
        # Update previous solution
        u_n.assign(u)
    return error_T

# Hold plot
for N in [8,16,32]:
    for num_steps in [4, 8, 16, 32, 64, 128]:
        E_L2 = heat(N, num_steps)
        print(N, num_steps, E_L2)

Running this, we observe the following output:

8 4 0.06539198110491778
8 8 0.042934331859209904
8 16 0.03211115627444427
8 32 0.026944719856269814
8 64 0.024468735933054046
8 128 0.023267784380047587
16 4 0.05233550568117215
16 8 0.028522261840029685
16 16 0.016710087456445907
16 32 0.010916975281129544
16 64 0.00811929850248597
16 128 0.0067813127318284375
32 4 0.04926874595797532
32 8 0.025199243181007972
32 16 0.013173630876173495
32 32 0.007176429318273071
32 64 0.004201096625917421
32 128 0.0027412031272620015

I.e. using the coarsest grid, it does not help refine the number of time steps after num_steps=32.
Similarly for the coarsest temporal discretization, it does not help refining the mesh, (the error is almost constant). However, if we look at the finest temporal discretization, and gradually increase the number of cells in the mesh, we observe convergence.

3

Dear @dokken ,
When i will change scalar function u to vector function u, V should be changed in to VectorFunctionSpace and how to assign source function f in the time loop.

Please clarify this.
In the previous problem , I am curious to know what were the wrong step and where it is , Please let me know. so that I will implement accordingly.

4

See for instance: Bitbucket

You can also just continue using an Expression for f.

As i said in the previous points, the temporal discretization (num_steps) dominate the error, thus refining the grid would not help. Therefore, using more time steps should resolve your problem.

5

Dear @dokken ,

Please see the code

from fenics import *
import numpy as np

set_log_level(LogLevel.ERROR)    
def heat(N, num_steps):    

    T = 2.0            # final time
    dt = T / num_steps # time step size
    alpha = 3          # parameter alpha
    beta = 1.2         # parameter beta

    # Create mesh and define function space
    mesh = UnitSquareMesh(N, N)
    V = VectorFunctionSpace(mesh, 'P', 1)

    # Define boundary condition
    u_D = Expression(('(pow(x[0],2)-x[0])*(pow(x[1],2)-x[1])*t','(-pow(x[0],2)-x[0])*(pow(x[1],2)-x[1])*pow(t,2)') ,
                 degree=4,  t=0)

    
    bc = DirichletBC(V, u_D, "on_boundary")

    # Define initial value
    u_n = interpolate(u_D, V)
    #u_n = project(u_D, V)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    t = 0

    f = Expression(('(pow(x[0],2)-x[0])*(pow(x[1],2)-x[1])-2*t*(pow(x[1],2)-x[1]+ pow(x[0],2)-x[0])', \
                    '-2*t*(pow(x[0],2)-x[0])*(pow(x[1],2)-x[1]) + 2*pow(t,2)*(pow(x[1],2)-x[1]+ pow(x[0],2)-x[0])'), \
                    degree=4, t=0 )

    F = inner(u,v)*dx + dt*inner(grad(u), grad(v))*dx - inner(u_n + dt*f,v)*dx
    a, L = lhs(F), rhs(F)

    # Time-stepping
    u = Function(V)
    error_T = 0

    for n in range(num_steps):

        # Update current time
        t += dt
        u_D.t = t
        f.t   = t
        #f.assign(2*beta*t -2 - 2*alpha)
        # Compute solution
        solve(a == L, u, bc)


        # Compute error at vertices
        error = errornorm(u_D, u)
        error_T += dt*error
        # Update previous solution
        u_n.assign(u)
    return error_T

# Hold plot
for N in [8,16,32]:
    for num_steps in [4, 8, 16, 32, 64, 128]:
        E_L2 = heat(N, num_steps)
        print(N, num_steps, E_L2)

8 4 0.3253144299364805
8 8 0.278399665188453
8 16 0.2562229962661123
8 32 0.2454549540116501
8 64 0.24015106572362374
8 128 0.23751917959250488
16 4 0.3318332807733057
16 8 0.2839529740324894
16 16 0.26132006926126433
16 32 0.2503306182192598
16 64 0.2449177223445922
16 128 0.2422317601892724
32 4 0.33354946841933775
32 8 0.2854143459371149
32 16 0.26266100208175563
32 32 0.2516130855383201
32 64 0.24617139904646562
32 128 0.24347115301669386

when I run my code

from fenics import *

def heat(N):    
    T = 2.0            # final time
    num_steps = 4     # number of time steps
    dt = T / num_steps # time step size
    alpha = 3          # parameter alpha
    beta = 1.2         # parameter beta

    # Create mesh and define function space
    mesh = UnitSquareMesh(N, N)
    V = VectorFunctionSpace(mesh, 'P', 1)

    # Define boundary condition
    u_D = Expression(('(pow(x[0],2)-x[0])*(pow(x[1],2)-x[1])*t','(-pow(x[0],2)-x[0])*(pow(x[1],2)-x[1])*pow(t,2)') ,
                 degree=4,  t=0)

    bc = DirichletBC(V, u_D, "on_boundary")

    # Define initial value
    u_n = interpolate(u_D, V)
    #u_n = project(u_D, V)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
    f = Expression(('(pow(x[0],2)-x[0])*(pow(x[1],2)-x[1])-2*t*(pow(x[1],2)-x[1]+ pow(x[0],2)-x[0])', \
                    '-2*t*(pow(x[0],2)-x[0])*(pow(x[1],2)-x[1]) + 2*pow(t,2)*(pow(x[1],2)-x[1]+ pow(x[0],2)-x[0])'), \
                    degree=4, t=0 )

    F = inner(u,v)*dx + dt*inner(grad(u), grad(v))*dx - inner(u_n + dt*f,v)*dx
                                 
    a = lhs(F)
    L = rhs(F)
    
    
    # Time-stepping
    u = Function(V)
    t = 0
    for n in range(num_steps):

        # Update current time
        t += dt
        u_D.t = t
        f.t   = t
        # Compute solution
        solve(a == L, u, bc)

        # Plot solution
        #plot(u)

       
        error_L2 = errornorm(u_D, u, 'L2')
        # Update previous solution
        u_n.assign(u)
                                 

    return error_L2

for N in [2, 4, 8,16, 32, 64]:
    print("{0:d} error_L2 ={1:.2e}".format(N, heat(N)))

the out put is

2 error_L2 =3.17e-01
4 error_L2 =3.21e-01
8 error_L2 =3.43e-01
16 error_L2 =3.50e-01
32 error_L2 =3.52e-01
64 error_L2 =3.52e-01

when You are changing u_D in time loop as u_D.t=t, and f.assign() when f is a scalar constant
in u as vector function problem u_D.t=t and f(?)

Kindly let me know how to assign f in the time loof in u as vector function and kindly inform difference between u_D.t=t, and f.assign() .

Thank You,
Debendra

6

You keep on changing what problem you are solving. I do not have time to sit down and go through your now quite complex u_D and f. I would suggest you start with the scalar example I supplied you, and gradually extend it to a vector problem.

7

Dear @dokken ,
Please see the code for vector problem, which is working

from fenics import *
import numpy as np

set_log_level(LogLevel.ERROR)    
def heat(N, num_steps):    

    T = 2.0            # final time
    dt = T / num_steps # time step size
    alpha = 3          # parameter alpha
    beta = 1.2         # parameter beta

    # Create mesh and define function space
    mesh = UnitSquareMesh(N, N)
    V =VectorFunctionSpace(mesh, 'P', 1)

    # Define boundary condition
    u_D = Expression(('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t*t','1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t*t'), 
                      degree=2, alpha=alpha, beta=beta, t=0)

    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

    # Define initial value
    u_n = interpolate(u_D, V)
    #u_n = project(u_D, V)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
   
    t = 0

    f = Constant((2*beta*t - 2 - 2*alpha, 2*beta*t - 2 - 2*alpha ))

    F = inner(u,v)*dx + dt*inner(grad(u), grad(v))*dx - inner((u_n + dt*f),v)*dx
    a, L = lhs(F), rhs(F)

    # Time-stepping
    u = Function(V)
    error_T = 0
   
    for n in range(num_steps):

        # Update current time
        t += dt
        u_D.t = t
        
        f.assign(Constant([2*beta*t -2 - 2*alpha, 2*beta*t - 2 - 2*alpha]))
        # Compute solution
        solve(a == L, u, bc)


        # Compute error at vertices
        error = errornorm(u_D, u)
        error_T += dt*error
        # Update previous solution
        u_n.assign(u)
    return error_T

# Hold plot
for N in [8,16,32]:
    for num_steps in [4, 8, 16, 32, 64, 128]:
        E_L2 = heat(N, num_steps)
        print(N, num_steps, E_L2)

the output is

8 4 0.09247822654901844
8 8 0.0607183144067209
8 16 0.04541203270680158
8 32 0.03810558825508058
8 64 0.034604018210649644
8 128 0.03290561623663456
16 4 0.07401358192796759
16 8 0.04033656952372397
16 16 0.023631632309345764
16 32 0.015438934502666543
16 64 0.01148242205916962
16 128 0.009590224436037588
32 4 0.06967652873487848
32 8 0.035637111468114664
32 16 0.018630327450786882
32 32 0.01014900367132535
32 64 0.005941247825202164
32 128 0.003876646639798319

it is working.

my doubt is when source function is of the form in Vector problem, for instance,

f =Expression(('x[0]*x[1]*t', '2*x[0]*t*t', , degree=2, t=0))

in the time loop I tried

f.assign(Expression([x[0]*x[1]*t, 2*x[0]*t*t]))

it says that t is defind in this scope.

Kindly suggest me ,in this scenario what to do ?

Thank You,
Debendra

8

Why dont you simpy do

f =Expression(('x[0]*x[1]*t', '2*x[0]*t*t', , degree=2, t=0))
# ....
for step in range(num_steps):
    t += dt
   f.t = t
9

Dear @dokken ,
when I tried for

f = Constant((2*beta*t - 2 - 2*alpha, 2*beta*t - 2 - 2*alpha ))

    for n in range(num_steps):

        # Update current time
        t += dt
        f.t = t

in the previous problem , i.e

from fenics import *
import numpy as np

set_log_level(LogLevel.ERROR)    
def heat(N, num_steps):    

    T = 2.0            # final time
    dt = T / num_steps # time step size
    alpha = 3          # parameter alpha
    beta = 1.2         # parameter beta

    # Create mesh and define function space
    mesh = UnitSquareMesh(N, N)
    V =VectorFunctionSpace(mesh, 'P', 1)

    # Define boundary condition
    u_D = Expression(('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t*t','1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t*t'), 
                      degree=2, alpha=alpha, beta=beta, t=0)

    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

    # Define initial value
    u_n = interpolate(u_D, V)
    #u_n = project(u_D, V)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
   
    t = 0

    f = Constant((2*beta*t - 2 - 2*alpha, 2*beta*t - 2 - 2*alpha ))

    F = inner(u,v)*dx + dt*inner(grad(u), grad(v))*dx - inner((u_n + dt*f),v)*dx
    a, L = lhs(F), rhs(F)

    # Time-stepping
    u = Function(V)
    error_T = 0
   
    for n in range(num_steps):

        # Update current time
        t += dt
        u_D.t = t
        
        f.t=t
        # Compute solution
        solve(a == L, u, bc)


        # Compute error at vertices
        error = errornorm(u_D, u)
        error_T += dt*error
        # Update previous solution
        u_n.assign(u)
    return error_T

# Hold plot
for N in [8,16,32]:
    for num_steps in [4, 8, 16, 32, 64, 128]:
        E_L2 = heat(N, num_steps)
        print(N, num_steps, E_L2)

the out put shows

8 4 0.23487372975997028
8 8 0.23523062490311686
8 16 0.23565645690141793
8 32 0.235725475260137
8 64 0.23575858789865342
8 128 0.2357745663902976
16 4 0.25938489860171965
16 8 0.258545925999806
16 16 0.2581682885143698
16 32 0.25808261069208344
16 64 0.2580333321270201
16 128 0.2580089778972719
32 4 0.26592361102628886
32 8 0.2650489259387714
32 16 0.2646124704658123
32 32 0.26439659229880275
32 64 0.2642990034035616
32 128 0.26425537078655464

the ereror is still increasing,
the same problem shows for my problem.

Thank You,
Debendra

10

Here you haven’t defined f as an Expression as I suggested.

11

Dear @dokken ,

Thanks, it is working ,
Kindly corect me if the error for L2 norm is true. If I go in this way,
Kindly see the code ,

from fenics import *
import numpy as np

set_log_level(LogLevel.ERROR)    
def heat(N, num_steps):    

    T = 2.0            # final time
    dt = T / num_steps # time step size
    alpha = 3          # parameter alpha
    beta = 1.2         # parameter beta

    # Create mesh and define function space
    mesh = UnitSquareMesh(N, N)
    V =VectorFunctionSpace(mesh, 'P', 1)

    # Define boundary condition
    u_D = Expression(('1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t*t','1 + x[0]*x[0] + alpha*x[1]*x[1] + beta*t*t'), 
                      degree=2, alpha=alpha, beta=beta, t=0)

    def boundary(x, on_boundary):
        return on_boundary

    bc = DirichletBC(V, u_D, boundary)

    # Define initial value
    u_n = interpolate(u_D, V)
    #u_n = project(u_D, V)

    # Define variational problem
    u = TrialFunction(V)
    v = TestFunction(V)
   
    t = 0

    f = Expression(('2*beta*t - 2 - 2*alpha', '2*beta*t - 2 - 2*alpha'), degree=0, alpha=alpha, beta=beta, t=0)

    F = inner(u,v)*dx + dt*inner(grad(u), grad(v))*dx - inner((u_n + dt*f),v)*dx
    a, L = lhs(F), rhs(F)

    # Time-stepping
    u = Function(V)
    error_T = 0
   
    for n in range(num_steps):

        # Update current time
        t += dt
        u_D.t = t
        
        f.t=t
        # Compute solution
        solve(a == L, u, bc)


        # Compute error at vertices
        error = errornorm(u_D, u , 'L2')
        #error_T += dt*error
        # Update previous solution
        u_n.assign(u)
    #return error_T
    return error

# Hold plot
for N in [8,16,32]:
    for num_steps in [4, 8, 16, 32, 64, 128]:
        E_L2 = heat(N, num_steps)
        print(N, num_steps, E_L2)

Out put is ,

8 4 0.047027897959540356
8 8 0.030740075568045148
8 16 0.022886472490036138
8 32 0.019138354440597413
8 64 0.01734319360591817
8 128 0.01647294543724075
16 4 0.03786306256872677
16 8 0.02059403675056871
16 16 0.012024321622422125
16 32 0.007819586443231547
16 64 0.005788648904869915
16 128 0.0048176909914040425
32 4 0.0357104343077472
32 8 0.01825519729141756
32 16 0.009532617676186328
32 32 0.005182098906917493
32 64 0.0030231423494833484
32 128 0.0019635176300912332

For the Source function the degree depends always on spatial variable?

Thank you,
Debendra

12

The degree of the source function should reflect which finite element function space the Expression is in. I.e. if you have a quadratic term, it should be at least 2.

1 Like
13

Thanks @dokken ,
In the Previous problem for Poisson ,
set_log_level(LogLevel.ERROR) is not applied ,

Here You are using, what is the role of
set_log_level(LogLevel.ERROR)
for this problem.

14

It controls how much output you get from dolfin.
The different levels are covered here: How to disable message "Solving linear variational problem." - FEniCS Q&A

1 Like
15

Hello dokken

I am trying to get a graph of the magnitudes of my solution u of the navier stokes problem, where u is a vector function. I would also like to make an animation of the solution like as gif, however, I have not found the answer, any suggestions?

import fenics as fe
from tqdm import tqdm # optional
import matplotlib.pyplot as plt
N_POINTS_P_AXIS = 41
TIME_STEP_LENGTH = 0.01
N_TIME_STEPS = 10000
KINEMATIC_VISCOSITY = 0.01 # -> Re = 100
def main():
    mesh = fe.UnitSquareMesh(N_POINTS_P_AXIS, N_POINTS_P_AXIS)
    velocity_function_space = fe.VectorFunctionSpace(mesh, "Lagrange", 2)
    pressure_function_space = fe.FunctionSpace(mesh, "Lagrange", 1)
    u_trial = fe.TrialFunction(velocity_function_space)
    p_trial = fe.TrialFunction(pressure_function_space)
    v_test = fe.TestFunction(velocity_function_space)
    q_test = fe.TestFunction(pressure_function_space)

    # Define the Boundary Condition
    stationary_wall_boundary_condition = fe.DirichletBC(
        velocity_function_space,
    (0.0, 0.0),
    """
    on_boundary && (x[0] < DOLFIN_EPS || x[1] < DOLFIN_EPS || x[0] > (1.0 - DOLFIN_EPS))
    """
    )
    moving_wall_boundary_condition = fe.DirichletBC(
    velocity_function_space,
    (1.0, 0.0),
    """
    on_boundary && (x[1] > (1.0 - DOLFIN_EPS))
    """
    )
    velocity_boundary_conditions = [stationary_wall_boundary_condition, moving_wall_boundary_condition]

    # Define the solution fields involved
    u_prev = fe.Function(velocity_function_space)
    u_tent = fe.Function(velocity_function_space)
    u_next = fe.Function(velocity_function_space)
    p_next = fe.Function(pressure_function_space)

    # Weak form of the momentum equation
    momentum_weak_form_residuum = (
    1.0 / TIME_STEP_LENGTH * fe.inner(u_trial - u_prev, v_test) * fe.dx
    +
    fe.inner(fe.grad(u_prev) * u_prev, v_test) * fe.dx
    +
    KINEMATIC_VISCOSITY * fe.inner(fe.grad(u_trial), fe.grad(v_test)) * fe.dx
    )
    momentum_weak_form_lhs = fe.lhs(momentum_weak_form_residuum)
    momentum_weak_form_rhs = fe.rhs(momentum_weak_form_residuum)

    # Weak form of the pressure poisson problem
    pressure_poisson_weak_form_lhs = fe.inner(fe.grad(p_trial), fe.grad(q_test)) * fe.dx
    pressure_poisson_weak_form_rhs = - 1.0 / TIME_STEP_LENGTH * fe.div(u_tent) * q_test * fe.dx

    # Weak form of the velocity update equation
    velocity_update_weak_form_lhs = fe.inner(u_trial, v_test) * fe.dx
    velocity_update_weak_form_rhs = (
    fe.inner(u_tent, v_test) * fe.dx
    -
    TIME_STEP_LENGTH * fe.inner(fe.grad(p_next), v_test) * fe.dx
    )

    # Pre-Compute the system matrices (because they do not greatly change)
    momentum_assembled_system_matrix = fe.assemble(momentum_weak_form_lhs)
    pressure_poisson_assembled_system_matrix = fe.assemble(pressure_poisson_weak_form_lhs)
    velocity_update_assembled_system_matrix = fe.assemble(velocity_update_weak_form_lhs)

    for t in tqdm(range(N_TIME_STEPS)):
        # (1) Solve for tentative velocity
        momentum_assembled_rhs = fe.assemble(momentum_weak_form_rhs)
        [bc.apply(momentum_assembled_system_matrix, momentum_assembled_rhs) for bc in velocity_boundary_conditions]
        fe.solve(
        momentum_assembled_system_matrix,
        u_tent.vector(),
        momentum_assembled_rhs,
        "gmres",
        "ilu",
        )

        # (2) Solve for the pressure
        pressure_poisson_assembled_rhs = fe.assemble(pressure_poisson_weak_form_rhs)
        fe.solve(
        pressure_poisson_assembled_system_matrix,
        p_next.vector(),
        pressure_poisson_assembled_rhs,
        "gmres",
        "amg",
        )

        # (3) Correct the velocities to be incompressible
        [bc.apply(momentum_assembled_system_matrix, momentum_assembled_rhs) for bc in velocity_boundary_conditions]
        velocity_update_assembled_rhs = fe.assemble(velocity_update_weak_form_rhs)
        fe.solve(
        velocity_update_assembled_system_matrix,
        u_next.vector(),
        velocity_update_assembled_rhs,
        "gmres",
        "ilu",
        )

        # Advance in time
        u_prev.assign(u_next)

        # Visualize interactively
        fe.plot(u_next)
        plt.draw()
        plt.pause(0.02)
        plt.clf()
16

I would suggest writing the solution to file and use paraview for such visualizations

17

Okay, in the case of animation. How can I export it to use it in a latex .tex?

18

These questions are unrelated to FEniCS. I suggest searching for LaTeX / beamer tutorials.