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()
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.