Results appear wrong for implementation of symmetric boundary conditions and initial conditions

[appana_3]

Hi everyone,

I am running FeNICSx on Windows using docker containers. I am implementing a transient heat conduction equation with square geometry, triangular mesh, uniform initial condition across the mesh, symmetric boundary conditions across the centroidal X-axis. The boundary conditions are only Neumann and Robin. The code below details various sections of implementation.

# Square geometry with triangular mesh
nx, ny = 50, 50
mesh = create_rectangle(MPI.COMM_WORLD, [np.array([0,0]),np.array([5,5])], [nx, ny], CellType.triangle)
V = FunctionSpace(mesh, ("Lagrange",1))

# Constant initial condition across the mesh
x = SpatialCoordinate(mesh)
T_n = Function(V)
T_n.interpolate(lambda x: np.full((x.shape[1],), 298))

# Boundary definitions
boundaries = [(1, lambda x: np.isclose(x[0],0)),
              (2, lambda x: np.isclose(x[1],5)),
              (3, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],4),np.less(x[1],5))),
              (4, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],3),np.less(x[1],4))),
              (5, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],2),np.less(x[1],3))),
              (6, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],1),np.less(x[1],2))),
              (7, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],0),np.less(x[1],1))),
              (8, lambda x: np.isclose(x[1],0))]

# Variational formulation
# T_inf, Q_gen, q_tabs, rho, C_p, k, h are all constant across the mesh at all timestamps
u, v = TrialFunction(V), TestFunction(V)
a = form(rho * C_p * u * v * dx + dt * inner(k * grad(u), grad(v)) * dx + 
     dt * h * u * v * (ds(1) + ds(2) + ds(3) + ds(5) + ds(7) + ds(8)))
L = form((rho * C_p * T_n + dt * Q_gen) * v * dx - dt * q_tabs * v * (ds(4) + ds(6)) +
     dt * h * T_inf * v * (ds(1) + ds(2) + ds(3) + ds(5) + ds(7) + ds(8)))

I am using a linear solver and I have simulated it for 50 seconds. Ideally, the results should be symmetric about x[1]=2.5

But I am not getting it that way. Below is the result I got. Can someone help point out what I am doing wrong?

Thank you

[dokken]

Your code is not complete, as it is missing several imports and definitions of your solver, ds etc.
Please set up a minimal reproducible example that can illustrate this by looking at the solution after 1 time step.

[appana_3]

Got it, @dokken . Find below a working example which shows the same.

# Imports
from dolfinx import default_scalar_type
import dolfinx
from dolfinx.fem import (Constant,  Function, FunctionSpace, assemble_scalar, form, locate_dofs_topological)
from dolfinx.fem.petsc import (LinearProblem, assemble_vector, assemble_matrix, create_vector,
                               apply_lifting, set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_rectangle, locate_entities, locate_entities_boundary, meshtags, CellType
from dolfinx.plot import vtk_mesh

from mpi4py import MPI
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, div, dot, dx, grad, 
                 inner, lhs, rhs)

from petsc4py import PETSc
from mpi4py import MPI

import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pyvista

# Mesh
nx, ny = 50, 50
mesh = create_rectangle(MPI.COMM_WORLD, [np.array([0,0]),np.array([5,5])], [nx, ny], CellType.triangle)
V = FunctionSpace(mesh, ("Lagrange",1))

# Temporal definitions
t = 0     # Start time
t_end = 1     # End time
num_steps = 1    # Number of time steps
dt = (t_end-t)/num_steps   # Time step size

# Initialization
x = SpatialCoordinate(mesh)
T_n = Function(V)
T_n.interpolate(lambda x: np.full((x.shape[1],), 298))

T_h = Function(V)
T_h.interpolate(lambda x: np.full((x.shape[1],), 298))

T_inf = Constant(mesh,default_scalar_type(298))    # Ambient temperature
Q_gen = Constant(mesh,default_scalar_type(1000))    # Heat generation inside CV
q_tabs = Constant(mesh,default_scalar_type(500))    # Heat flux to tabs
rho = Constant(mesh,default_scalar_type(25))     # Density
C_p = Constant(mesh,default_scalar_type(90))      # Specific heat capacity of the cell
k = Constant(mesh,default_scalar_type(25))         # Thermal conductivity
h = Constant(mesh,default_scalar_type(10))          # Convection coefficient

# Creating boundaries
boundaries = [(1, lambda x: np.isclose(x[0],0)),
              (2, lambda x: np.isclose(x[1],5)),
              (3, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],4),np.less(x[1],5))),
              (4, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],3),np.less(x[1],4))),
              (5, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],2),np.less(x[1],3))),
              (6, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],1),np.less(x[1],2))),
              (7, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],0),np.less(x[1],1))),
              (8, lambda x: np.isclose(x[1],0))]

facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
    facets = locate_entities(mesh, fdim, locator)
    facet_indices.append(facets)
    facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)

ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)

# Variational formulation
u, v = TrialFunction(V), TestFunction(V)
a = form(rho * C_p * u * v * dx + dt * inner(k * grad(u), grad(v)) * dx + 
     dt * h * u * v * (ds(1) + ds(2) + ds(3) + ds(5) + ds(7) + ds(8)))
L = form((rho * C_p * T_n + dt * Q_gen) * v * dx - dt * q_tabs * v * (ds(4) + ds(6)) +
     dt * h * T_inf * v * (ds(1) + ds(2) + ds(3) + ds(5) + ds(7) + ds(8)))

# Solving
A = assemble_matrix(a, bcs=[])
A.assemble()
b = create_vector(L)

solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)

for i in range(num_steps):
    t += dt
    
    # Update the right hand side reusing the initial vector
    with b.localForm() as loc_b:
        loc_b.set(0)
    assemble_vector(b, L)
    
    # Apply Dirichlet boundary condition to the vector
    apply_lifting(b, [a], [[]])
    b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
    set_bc(b, [])

    # Solve linear problem
    solver.solve(b, T_h.vector)
    T_h.x.scatter_forward()

    # Update solution at previous time step (u_n)
    T_n.x.array[:] = T_h.x.array

# Plotting the solution
pyvista.start_xvfb()
pyvista_cells, cell_types, geometry = vtk_mesh(V)
grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry)
grid.point_data["u"] = T_n.x.array
grid.set_active_scalars("u")

plotter = pyvista.Plotter()
plotter.add_text("uh", position="upper_edge", font_size=14, color="black")
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("robin_neumann_dirichlet.png")

The solution I got is:

[appana_3]

Hi, can someone help me find out why I am getting erroneous results?

[jkrokowski]

It’s not clear to anyone what you are trying to do and what you would expect to see. You need to provide more information on the problem you are trying to solve and what boundary conditions you are trying to implement. There is not enough information here for someone to help you. One place to start might be here

I would suggest simplifying your problem and clearly writing out (in math) what you are trying to do.

[appana_3]

Hello @jkrokowski ,

Thank you for pointing it out! Below is the information on the problem I am trying to solve.

The picture below shows the geometry and the boundary definitions of the problem.

The PDE I am trying to solve is the transient heat conduction equation

Boundary condition at edges 1, 2, 3, 4, 6, 8:
−𝑘 𝑑𝑇/𝑑𝑛= ℎ(𝑇−𝑇_𝑖𝑛𝑓)

Boundary condition at edges 5 and 7:
−𝑘 𝑑𝑇/𝑑𝑛=𝑞

I have assumed that the terms rho, c_p, k, Q_gen, Q_conv, h, T_inf, q are all constant for simplicity.

As I have shown in the code above, my variational formulation is:

a = form(rho * C_p * u * v * dx + dt * inner(k * grad(u), grad(v)) * dx + 
     dt * h * u * v * (ds(1) + ds(2) + ds(3) + ds(5) + ds(7) + ds(8)))
L = form((rho * C_p * T_n + dt * Q_gen) * v * dx - dt * q_tabs * v * (ds(4) + ds(6)) +
     dt * h * T_inf * v * (ds(1) + ds(2) + ds(3) + ds(5) + ds(7) + ds(8)))

where u is my trial function, v is my test function, T_n is the solution of the previous time step.

As you can see from the geometry and the equations, the solution should be symmetric about the central horizontal axis. But I am getting weird results.

I hope this makes my problem clear!

[appana_3]

Hello forum,
Can someone help me point out what the issue is?

[wanghao]

Hi appana_3!
I don’t know what kind of outcome you want.

Is it because your result boundary 1,2,3,4,6,8 has no temperature change？
if so, perhaps it’s because your initial temperature value is different from T_ inf is consistent, so there is no heat exchange, you just need to replace T_inf with other values can result in such a result:

[appana_3]

Hi @wanghao ,

Thank you for responding.

No, that is not my concern! I understand that the boundaries 1,2,3,4,6,8 will have temperature change if I change T_inf.

My concern is that: theoretically, T(x=a, y=b) should be equal to T(x=a, y=5-b). In simple terms, if I compare just 2 points, the temperature at the bottom right corner should be equal to the temperature at top right corner. Also, in the results you have shared, these temperatures are not equal to each other.

This is because of the symmetry in the problem. Let’s say we flip the geometry along y=2.5 line. Then everything in the new problem is just the same as earlier. But the new temperature of the earlier bottom right corner or the earlier top right corner are changing. I am not able to figure out why that is happening.

The solution I am expecting should have a symmetric temperature distribution about y=2.5 line. Hope it makes my problem clear.

[raboynics]

I don’t see symmetry between points 4 and 8 in your weak form. Therefore there’s no reason to expect a symmetry about y = 1.5. Am I missing something?

[appana_3]

Hello @raboynics , thank you for your response!

My apologies. I made a mistake in the label edges in the geometry I uploaded above.

Please find below the correctly labeled geometry. This is what I have implemented in my weak formulation. In this, you can see that there should be symmetry between edges 3 & 7 about y=2.5

[appana_3]

Hi everyone,

I was cross checking the implementation of boundary conditions using Paraview and I think I finally see the problem. It looks like the way I have designated ids for the boundaries is not working correctly.

# Creating boundaries
boundaries = [(1, lambda x: np.isclose(x[0],0)),
              (2, lambda x: np.isclose(x[1],5)),
              (3, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],4),np.less(x[1],5))),
              (4, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],3),np.less(x[1],4))),
              (5, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],2),np.less(x[1],3))),
              (6, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],1),np.less(x[1],2))),
              (7, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],0),np.less(x[1],1))),
              (8, lambda x: np.isclose(x[1],0))]

As I have shown in the geometry earlier, the boundary id ‘5’ should be for (x=5, 2<y<3). When I see in paraview (screenshot attached below), the boundary id ‘5’ extends from y=2 to y=5.

Can someone tell me how to define a boundary to be for (x=5, 2<y<3)?

I have tried (5, lambda x: np.logical_and(np.isclose(x[0],5),x[1]>=2,x[1]<3)). It didn’t work.

Thank you!

[spcvanschie]

My np.logical_and() function doesn’t like being given three arrays as input arguments. One way around this that I’ve quickly tried is to combine two of the input arrays into a separate np.logical_and() instance. For example:

(5, lambda x: np.logical_and(np.isclose(x[0],5), np.logical_and(np.greater_equal(x[1],2),np.less(x[1],3))))

This is obviously not an elegant way of writing things. If I were you I would write a function that wraps some of this code to make it more general (i.e. let it specify arbitrary boundary intervals). However, the above implementation should work (it did for me in a quick test at least).

Edit: A suggestion by @jkrokowski is to use np.logical_and.reduce(), as follows:

(5, lambda x: np.logical_and.reduce([np.isclose(x[0],5), np.greater_equal(x[1],2), np.less(x[1],3)]))

[appana_3]

It worked! My original problem also got resolved!

Thank you so much everyone

[dokken]

You can also use the bitwise operation:
https://numpy.org/doc/stable/reference/generated/numpy.bitwise_and.html
np.isclose(x[0],5) & np.greater_equal(x[1],2) & np.less(x[1],3)])

[appana_3]

Thank you @dokken!

I have one question. In the last command of this bitwise operation, should I use np.less or np.less_equal? Which is more accurate?

np.isclose(x[0],5) & np.greater_equal(x[1],2) & np.less(x[1],3))

or

np.isclose(x[0],5) & np.greater_equal(x[1],2) & np.less_equal(x[1],3))

This boundary extends from y=2 to y=3

[dokken]

@appana_3

I have one question. In the last command of this bitwise operation, should I use np.less or np.less_equal? Which is more accurate?

This is just a question about what you want to include in your boundary. If there is a node at 3, and you want to include it,you should use less_equal.

[French_fries]

你好，我也在求解导热问题，我想问一下如何将模型的某个边界设置为对称边界呢，fenics里面好像没有