Getting same results at timesteps 0 and 1

dolfinx
1

Hello dear FEniCSX users,

I wanted to solve the following PDE:

\dfrac{\partial T}{\partial t} = \Delta T ,
T(x=1, y, t) = 1 ,
T(x, y=1, t) = 1 ,
T(x, y, 0) = 0 .

I have used the tutorials available in the FEniCS discourse when writing my code.

Here is my code:


#Below are my references
#https://fenicsproject.discourse.group/t/upadate-solutions-in-time-dependent-non-linear-problem/15719/4
#https://fenicsproject.discourse.group/t/upadate-solutions-in-time-dependent-non-linear-problem/15719
#https://fenicsproject.discourse.group/t/time-dependent-diffusion-constant-in-dolfinx/8137/5


#code to solve 2D Heat Transfer

from mpi4py import MPI
import dolfinx
from dolfinx import default_scalar_type as dst
from dolfinx import mesh, fem, io, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver


import ufl
import basix

import numpy as np
import pathlib

cache_dir = pathlib.Path('cache_2HeatTransfer2024')



Tf = 67.         #Final time
num_step = 100 
dt = Tf/num_step

T_r = 1.0    #boundary value for the dimensionless temperature at right boundary
T_t = 1.0     #boundary value for the dimensionless temperature at top boundary
T_init = 0.0  ##initial condition for dimensionless temperature

L = 1.0  #total length


# Create domain and element

# Define mesh
nx, ny = 100, 100
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0.0, 0.0]), np.array([L, L])],
                               [nx, ny], mesh.CellType.triangle)

# Function space
V = fem.functionspace(domain, ("Lagrange", 1))



# Get boundary dofs
tdim = domain.topology.dim
fdim = tdim - 1


def right_boundary(x):
    return np.isclose(x[0], L)


right_boundary_facets = mesh.locate_entities_boundary(domain, fdim,
                                                     right_boundary)
right_boundary_dofs = fem.locate_dofs_topological(V, fdim,
                                                 right_boundary_facets)



def top_boundary(x):
    return np.isclose(x[1], L)


top_boundary_facets = mesh.locate_entities_boundary(domain, fdim,
                                                      top_boundary)
top_boundary_dofs = fem.locate_dofs_topological(V, fdim,
                                                  top_boundary_facets)


# Create functions
T = fem.Function(V)  # current solution
v = ufl.TestFunction(V)

T_n = fem.Function(V)  # previous solution

# Initial condition
T_init_expr = fem.Expression(fem.Constant(domain,
                                          dst(T_init)),
                             V.element.interpolation_points(),
                             jit_options={'cache_dir': cache_dir})
T_n.interpolate(T_init_expr)
T_n.x.scatter_forward()






dt = fem.Constant(domain, dst(dt))

A = ufl.dot(ufl.grad(T), ufl.grad(v)) * ufl.dx    
lf = ((T - T_n) / dt) * v * ufl.dx

# BC
bc_0 = fem.dirichletbc(dst(T_r), right_boundary_dofs, V)
bc_L = fem.dirichletbc(dst(T_t), top_boundary_dofs, V)



# Problem
F = A + lf  

problem = NonlinearProblem(F, T, [bc_0, bc_L],
                           jit_options={'cache_dir': cache_dir})

# Solver
solver = NewtonSolver(MPI.COMM_WORLD, problem)


info = solver.solve(T)
print(info)

solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.max_it = 100
solver.report = True

'''
with dolfinx.io.VTXWriter(domain.comm, "C_moisture.bp", V, engine="BP4") as vtx:
    vtx.write(0.0)

'''

# Output
T_out = io.XDMFFile(domain.comm,
                    '2DHeatTransfer_Dec2024.xdmf',
                    'w')
T_out.write_mesh(domain)

# Time loop
t = 0
while t < Tf:
    if t == 0:
        T_out.write_function(T, t)

    t += dst(dt)
    print(f'Time {t}')
    n, converged = solver.solve(T)
    print(f'converged: {converged}, n_iter: {n}')
    T_n.interpolate(T)  # update prev. sol
    T_n.x.scatter_forward()

    # Save
    T_out.write_function(T, t)

T_out.close()

The code works well and converges in two iterations. But when I visualize the results in Paraview, the results at timesteps 0 and 1 are the same.

I am not sure if my initial condition is well implemented. I am a quite loss with this problem, and any suggestions would be appreciated.

Thank you very much in advance!

2

If you remove your line info = solver.solve(T), then it works as expected.

Right now, your first export is exporting the T from that solve. Then you do not update T_n and re-run that same computation. That gives you the exact same T, which you export again. So that’s twice the same solution in the vtk file.

1 Like
3

Thank you @Stein for your answer. When I removed the line info = solver.solve(T), I see that my solution now starts at 0 at timestep 0 and evolves over time until it reaches 1.

But the visualisation does not show the same.
Please can you upload your screenshot at timestep 1 because I get this plot (see screenshot at timestep 1).

Thank you

FEniCSX_picture_mamadou_stein_Mon9Dec2024
FEniCSX_picture_mamadou_stein_Mon9Dec20241551×798 15.3 KB

for your help.

4

Ah, you’re getting confused about the paraview colorbar settings.

Every timestep, the solution values change, but paraview doesn’t automatically update the colorbar range. On the left of the top settingsbar, there are buttons for rescaling the colorbar. If you click “rescale to visible data range” (the button with the eye symbol), then you will see the color variation in your solution at the first timestep.

You can also click ‘rescale to range over all timesteps’ to keep it constant and usable for all timesteps, but this is often also not great.

1 Like
5

Thank you again @Stein for your time and patience. I really learned something new today.

1 Like
6

Dear @Stein,
I wanted to change the problem to

\dfrac{\partial T}{\partial t} = \Delta T ,
T(x=1, y, t) = 0 ,
T(x, y =1, t) = 0 ,
T(x, y, t=0) = 1 .

But after visualisation with Paraview, the solution at the initial condition is zero instead of 1.

7

Could you post your updated code? I suspect you’re (correctly) setting T_n but are outputting T.

8

Of course, here it is

#Below are my references
#https://fenicsproject.discourse.group/t/upadate-solutions-in-time-dependent-non-linear-problem/15719/4
#https://fenicsproject.discourse.group/t/upadate-solutions-in-time-dependent-non-linear-problem/15719
#https://fenicsproject.discourse.group/t/time-dependent-diffusion-constant-in-dolfinx/8137/5


#code to solve 2D Heat Transfer

from mpi4py import MPI
import dolfinx
from dolfinx import default_scalar_type as dst
from dolfinx import mesh, fem, io, log
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver


import ufl
import basix

import numpy as np
import pathlib

cache_dir = pathlib.Path('cache_2HeatTransfer2024')



Tf = 67.         #Final time
num_step = 100 
dt = Tf/num_step

T_r = 0.0    #boundary value for the dimensionless temperature at right boundary
T_t = 0.0     #boundary value for the dimensionless temperature at top boundary
T_init = 1.0  ##initial condition for dimensionless temperature

L = 1.0  #total length


# Create domain and element

# Define mesh
nx, ny = 100, 100
domain = mesh.create_rectangle(MPI.COMM_WORLD, [np.array([0.0, 0.0]), np.array([L, L])],
                               [nx, ny], mesh.CellType.triangle)

# Function space
V = fem.functionspace(domain, ("Lagrange", 1))



# Get boundary dofs
tdim = domain.topology.dim
fdim = tdim - 1


def right_boundary(x):
    return np.isclose(x[0], L)


right_boundary_facets = mesh.locate_entities_boundary(domain, fdim,
                                                     right_boundary)
right_boundary_dofs = fem.locate_dofs_topological(V, fdim,
                                                 right_boundary_facets)



def top_boundary(x):
    return np.isclose(x[1], L)


top_boundary_facets = mesh.locate_entities_boundary(domain, fdim,
                                                      top_boundary)
top_boundary_dofs = fem.locate_dofs_topological(V, fdim,
                                                  top_boundary_facets)


# Create functions
T = fem.Function(V)  # current solution
v = ufl.TestFunction(V)

T_n = fem.Function(V)  # previous solution

# Initial condition
T_init_expr = fem.Expression(fem.Constant(domain,
                                          dst(T_init)),
                             V.element.interpolation_points(),
                             jit_options={'cache_dir': cache_dir})
T_n.interpolate(T_init_expr)
T_n.x.scatter_forward()






dt = fem.Constant(domain, dst(dt))

A = ufl.dot(ufl.grad(T), ufl.grad(v)) * ufl.dx    
lf = ((T - T_n) / dt) * v * ufl.dx

# BC
bc_0 = fem.dirichletbc(dst(T_r), right_boundary_dofs, V)
bc_L = fem.dirichletbc(dst(T_t), top_boundary_dofs, V)



# Problem
F = A + lf  

problem = NonlinearProblem(F, T, [bc_0, bc_L],
                           jit_options={'cache_dir': cache_dir})

# Solver
solver = NewtonSolver(MPI.COMM_WORLD, problem)


solver.convergence_criterion = "incremental"
solver.rtol = 1e-6
solver.max_it = 100
solver.report = True



# Output
T_out = io.XDMFFile(domain.comm,
                    '2DHeatTransfer_Dec2024.xdmf',
                    'w')
T_out.write_mesh(domain)

# Time loop
t = 0
while t < Tf:
    if t == 0:
        T_out.write_function(T, t)

    t += dst(dt)
    print(f'Time {t}')
    n, converged = solver.solve(T)
    print(f'converged: {converged}, n_iter: {n}')
    T_n.interpolate(T)  # update prev. sol
    T_n.x.scatter_forward()

    # Save
    T_out.write_function(T, t)

T_out.close()

I just changed

T_r = 1.0    #boundary value for the dimensionless temperature at right boundary
T_t = 1.0     #boundary value for the dimensionless temperature at top boundary
T_init = 0.0  ##initial condition for dimensionless temperature

to

T_r = 0.0    #boundary value for the dimensionless temperature at right boundary
T_t = 0.0     #boundary value for the dimensionless temperature at top boundary
T_init = 1.0  ##initial condition for dimensionless temperature
9

Right, and then you interpolate that initial condition expression on T_n but you’re exporting T. So that first export line should export T_n

10

Thank once again @Stein, this has solved the problem.