How to include the values on the surfaces?

[CalebDmArcher]

Could you please teach me how to include the values on the surfaces (boundaries)?

The code below shows a cube with constant values on the left and the right surfaces. The values are 25 on the left and 0 on the right. No other boundary conditions on the rest of the surfaces.

The image below shows my result. As you can see, the max value is slightly smaller than 25 and the min value is slightly bigger than 0.
By checking the plot of the front surface value, I just realized that the values are located at the center of each mesh grid and therefore the values on the surfaces are not included. I assume that vol_slice = phi.faceValue[mesh.facesFront] does not exactly get the values on the surface, either. It just returns the values on the centers of the grid cells laid on the front surface. More specifically, for the front surface plot, it will be better if I can see the min value is 0 and the max value is 25 which are the values right on the left and the right surfaces.

Is there a way to get the values on the surfaces? I also need to extract the values on the surfaces for other cases.

# Solver Setup
import os
# os.environ["FIPY_SOLVERS"] = "scipy"
os.environ["FIPY_SOLVERS"] = "petsc"
print("FIPY_SOLVERS is set to:", os.environ["FIPY_SOLVERS"])
os.environ["OMP_NUM_THREADS"] = "1"
print("OMP_NUM_THREADS is set to:", os.environ["OMP_NUM_THREADS"])

# Solver definitions
import fipy as fp
# Import petsc4py and set PETSc options
from petsc4py import PETSc
PETSc.Options().setValue('-ksp_monitor_true_residual', '')
PETSc.Options().setValue('-ksp_converged_reason', '')
PETSc.Options().setValue('-pc_type', 'hypre')
PETSc.Options().setValue('-mg_levels_pc_type', 'sor')
PETSc.Options().setValue('-ksp_rtol', '1.e-12')
PETSc.Options().setValue('-ksp_max_it', '10000')
PETSc.Options().setValue('-ksp_type', 'fgmres')
print("The solver package is being used is", fp.solvers.solver)

from fipy import CellVariable, Grid3D, DiffusionTerm, FaceVariable, Viewer
import numpy as np
import matplotlib.pyplot as plt
plt.switch_backend('TkAgg')
from mpl_toolkits.mplot3d import Axes3D

# Define the mesh dimensions
L = 1  # Length in meters
w = 1  # Width in meters
t = 1  # Thickness in meters

nx = 50  # Number of grid points in x-direction (length)
ny = 50  # Number of grid points in y-direction (width)
nz = 50   # Number of grid points in z-direction (thickness)

dx = L / nx  # Size of each cell in x-direction
dy = w / ny  # Size of each cell in y-direction
dz = t / nz  # Size of each cell in z-direction

mesh = Grid3D(nx=nx, ny=ny, nz=nz, dx=dx, dy=dy, dz=dz)

# Define the solution variables
phi = CellVariable(name='phi', mesh=mesh, value=0.)

# Define the diffusion coefficient
gamma = 5.8e+7

# Set spatial variables for calculation and plot
xc, yc, zc = mesh.cellCenters
x, y, z = mesh.faceCenters

# Apply boundary condition
phi.constrain(25., where=mesh.facesLeft)
phi.constrain(0., where=mesh.facesRight)

# Define the equation for steady-state
eq = DiffusionTerm(coeff=gamma)

# Solve the equation
eq.solve(var=phi)

# Plot the results at the mid-plane in 2D (Voltage Distribution in mV)
if __name__ == '__main__':
    # Check the min and max values
    print("Min:", np.nanmin(phi.value))
    print("Max:", np.nanmax(phi.value))

    viewer = Viewer(vars=phi, datamin=0., datamax=25.)
    viewer.plot()

    # plot the values on the back face
    mask_slice = mesh.facesFront
    vol_slice = phi.faceValue[mask_slice]
    x_slice = x[mask_slice]
    y_slice = y[mask_slice]
    vol_slice = np.round(vol_slice, 2)
    # scatter plot
    plt.figure()
    plt.scatter(x_slice, y_slice, c=vol_slice, cmap='viridis', s=15)
    cbar = plt.colorbar(label='Solution Variable')
    plt.xlabel('x (m)')
    plt.ylabel('y (m)')
    plt.title('Values on Front Surface')

    plt.show()

[guyer]

Please stop putting questions in issues. That's what discussions is for.

[CalebDmArcher]

Sorry for the late reply. I had a procedure recently. I will post my questions in the discussion session from now on.

[guyer]

I assume that vol_slice = phi.faceValue[mesh.facesFront] does not exactly get the values on the surface, either. It just returns the values on the centers of the grid cells laid on the front surface. More specifically, for the front surface plot, it will be better if I can see the min value is 0 and the max value is 25 which are the values right on the left and the right surfaces.

Your interpretation is incorrect. phi.faceValue[mesh.facesFront] does give the values on the front faces. The issue is, that because you have specified no constraints on this boundary, it has a default condition of zero flux, which amounts to projecting the value of the adjacent cell center.

You may find this notebook instructive.

[CalebDmArcher]

Thank you for the response and the reference.

I tried print("FaceMin:", np.nanmin(phi.faceValue)) and print("FaceMax:", np.nanmax(phi.faceValue)) and then it returns 0 and 25.

I got confused by print("Min:", np.nanmin(phi.value)) and print("Max:", np.nanmax(phi.value)) returning 0.24999999999948305 and 24.750000000000778 previously. They are the values on the cell centers, not the surfaces.

[CalebDmArcher]

I assume that vol_slice = phi.faceValue[mesh.facesFront] does not exactly get the values on the surface, either. It just returns the values on the centers of the grid cells laid on the front surface. More specifically, for the front surface plot, it will be better if I can see the min value is 0 and the max value is 25 which are the values right on the left and the right surfaces.

Your interpretation is incorrect. phi.faceValue[mesh.facesFront] does give the values on the front faces. The issue is, that because you have specified no constraints on this boundary, it has a default condition of zero flux, which amounts to projecting the value of the adjacent cell center.

You may find this notebook instructive.

From there, I worry about my other simulation which I need to extract the value on the surface where only the Robin boundary condition is applied (no constrain() functions are used on this surface). I always see the simulation results I got were slightly smaller than the theoretical values. Could you have a look at it? The code is attached below.

A brief description of this simulation is that there is a structure that has a constant temperature on the bottom and a constant heat source on the top surface (only on the surface which means that no thickness is associated with this power). So only temp.constrain(T_ambient, bottom_faces) is applied on the bottom surface. The free air convection is applied on the top surface as well as the heat power is applied on the top surface so the Robin boundary condition is used. The Robin boundary condition is based on Applying Robin boundary conditions. No specific boundary conditions are applied on the other surfaces.

A good thing about this setup is that I can calculate the results by hand (because it can be considered as a 1D problem) and compare them with the results from the simulation. The manually calculated results and the simulation results are shown below. The plot is the temperature on the top surface vs. one of the parameters (convection coefficient, which is h in the code).

As you may see, the top surface temperature value is always slightly smaller than the hand-calculated values for every h value with different mesh sizes. As I make the mesh smaller, the differences go smaller. (I also checked the value on the bottom boundary which is set as 25 by using the constrain() method and I can see that the values are indeed always 25 which is great.) I initially thought the difference in the values on the top surface between the theoretical value and the simulation value was because of the mesh property. But your last post showed me that the solver might just use the value of the adjacent cells for the surface if no specific constraints are set for that surface. Do you think this is the case (as constrain() function is not used)?

example08_20_24.zip

[CalebDmArcher]

I assume that vol_slice = phi.faceValue[mesh.facesFront] does not exactly get the values on the surface, either. It just returns the values on the centers of the grid cells laid on the front surface. More specifically, for the front surface plot, it will be better if I can see the min value is 0 and the max value is 25 which are the values right on the left and the right surfaces.

Your interpretation is incorrect. phi.faceValue[mesh.facesFront] does give the values on the front faces. The issue is, that because you have specified no constraints on this boundary, it has a default condition of zero flux, which amounts to projecting the value of the adjacent cell center.
You may find this notebook instructive.

From there, I worry about my other simulation which I need to extract the value on the surface where only the Robin boundary condition is applied (no constrain() functions are used on this surface). I always see the simulation results I got were slightly smaller than the theoretical values. Could you have a look at it? The code is attached below.

A brief description of this simulation is that there is a structure that has a constant temperature on the bottom and a constant heat source on the top surface (only on the surface which means that no thickness is associated with this power). So only temp.constrain(T_ambient, bottom_faces) is applied on the bottom surface. The free air convection is applied on the top surface as well as the heat power is applied on the top surface so the Robin boundary condition is used. The Robin boundary condition is based on Applying Robin boundary conditions. No specific boundary conditions are applied on the other surfaces.

A good thing about this setup is that I can calculate the results by hand (because it can be considered as a 1D problem) and compare them with the results from the simulation. The manually calculated results and the simulation results are shown below. The plot is the temperature on the top surface vs. one of the parameters (convection coefficient, which is h in the code).

As you may see, the top surface temperature value is always slightly smaller than the hand-calculated values for every h value with different mesh sizes. As I make the mesh smaller, the differences go smaller. (I also checked the value on the bottom boundary which is set as 25 by using the constrain() method and I can see that the values are indeed always 25 which is great.) I initially thought the difference in the values on the top surface between the theoretical value and the simulation value was because of the mesh property. But your last post showed me that the solver might just use the value of the adjacent cells for the surface if no specific constraints are set for that surface. Do you think this is the case (as constrain() function is not used)?

example08_20_24.zip

To further address whether the difference between the theoretical values and the simulated results is because no specific boundary conditions are applied, I also tried to capture the temperature on the bottom side of the top copper layer. Based on the notebook, only the values on the surface approximated to the values of the adjacent cells and values inside the structure should be the 'real' values. However, right after calculating the theoretical values on the overlapped face (between the top copper layer and the second layer), I realized that these between-layer theoretical values are even greater than the surface values from the simulation.

This might imply that the mesh quality is a bigger issue leading to inaccuracy.

May I know what you think about it?

[guyer]

FiPy numerically solves PDEs on cell centers, not on face centers. It is therefore not particularly meaningful to talk about differences between values at face centers and some analytical solution. It is meaningful to compare the cell center values to your analytical solution. If refining the mesh gets you closer to the answer you want, then refine the mesh.