Is there a example for Robin boundary condition?

[CalebDmArcher]

Hi, I am trying to include a Robin boundary condition in my code. Are there any examples of it?

I am trying to simulate free air convection on a material boundary for a thermal model.
for the cells right on the surfaces, there is the boundary condition which can be illustrated by the equation 'k*(dT/dx) = h*(T-T_ambient)'.
T is the solution variable that the program needs to be calculated. Other variables in the equation (k, h, dx, T_ambient) are known. A brief explanation for this equation is that the left side of the equation is the heat conduction between the current cell and the adjacent cell, and the right side is the heat convection between the current cell and the ambient with a value T_ambient.

The above situation is neither Dirichlet nor Neumann boundary conditions, i.e. it is not a constant value on the boundary nor a constant derivative of the solution variable on the boundary. Instead, the derivative of the solution variable on the boundary depends on the value of the solution variable itself, which is dT/dx = (h/k)*(T-T_ambient).

Plus, since it is free air convection, the convection happens in every direction.

Could anyone help me?
I need to implement the thermal model for a 3D system. So a 2D or 3D example would be great!

[guyer]

Please see Applying Robin boundary conditions

[guyer]

You appear to have asked a question about how to apply this, but I no longer see it. Please repost if you still have questions.

[guyer]

I don't see results anything like this. What OS are you running on? What packages do you have installed for FiPy? If you run python setup.py test, it should report near the top (I suppose it might come out near the bottom) something like

Solver suite is petsc
{'python': '3.11.4 | packaged by conda-forge | (main, Jun 10 2023, 18:08:41) [Clang 15.0.7 ]', 'fipy': '3.4.4+2.gb3d357617', 'numpy': '1.24.3', 'pysparse': 'not installed', 'scipy': '1.10.1', 'matplotlib': '3.7.1', 'mpi4py': '3.1.4', 'petsc4py': '3.18.4', 'pyamgx': 'not installed', 'PyTrilinos': 'not installed', 'mayavi': '4.8.1', 'gmsh': '4.11.1', 'solver': 'petsc'}

What does it say for you?

[CalebDmArcher]

I am using Python 3.9 with PyCharm on Win11.
I do not know how to run python setup.py test which reports '[Errno 2] No such file or directory', but I can find some of them through PyCharm:

Package	Version
FiPy	3.4.4
gmsh	4.11.1
matplotlib	3.8.0
mayavi	4.8.1
numpy	1.26.0rc1
scipy	1.11.2

[CalebDmArcher]

I think I was using the default solver, LinearLUSolver, since I didn't specify another solver.

[guyer]

python setup.py test would be run from the command-line, within the installation directory for fipy. With an installed version of fipy, you can execute the python commands

import fipy
fipy.test()

[CalebDmArcher]

This is what I saw:

python      3.9.2 (tags/v3.9.2:1a79785, Feb 19 2021, 13:44:55) [MSC v.1928 64 bit (AMD64)]
fipy        3.4.4
numpy       1.26.0rc1
pysparse    not installed
scipy       1.11.2
matplotlib  3.8.0
mpi4py      not installed
petsc4py    not installed
pyamgx      not installed
PyTrilinos  not installed
mayavi      4.8.1
gmsh        4.11.1
solver      scipy

[guyer]

with some tiny modifications to the code:

D = FaceVariable(mesh=mesh, value=0.9)
D.setValue(400., where=z >= Lz/2)

You must have made more modifications than this.

• If you still have

  x, y, z = mesh.cellCenters

  then D.setValue() will fail.
• If you instead have

  x, y, z = mesh.faceCenters

  then phi_slice = phi.value[mask_slice] will later fail (as will subsequent steps).

Please don't make me guess what you're doing.

[CalebDmArcher]

Sorry for the inconvenience and confusion. I will try to post the full code from now on:

from fipy import CellVariable, Gmsh3D, DiffusionTerm, Viewer, FaceVariable, ImplicitSourceTerm
from fipy.tools import numerix
import numpy as np
import matplotlib.pyplot as plt
import time
from scipy.interpolate import griddata
from fipy import DefaultSolver


print(DefaultSolver)

Lx = Ly = Lz = 20
dx = dy = dz = 1.
nx = ny = nz = int(Lx / dx)


# modify the Gmesh file
# cube_template.geo uses hexahedron mesh
# cube_template2.geo uses tetrahedron mesh
def generate_geo_from_template_tetra(Lx, Ly, Lz, dx, template_path="cube_template.geo", output_path="cube.geo"):
    with open(template_path, 'r') as file:
        content = file.read()

    content = content.format(Lx=Lx, Ly=Ly, Lz=Lz, dx=dx, dy=dy, dz=dz, nx=nx, ny=ny, nz=nz)  # replace placeholders

    with open(output_path, 'w') as file:
        file.write(content)


generate_geo_from_template_tetra(Lx, Ly, Lz, dx)

mesh = Gmsh3D("cube.geo")

# Parameters
h = 10.0  # Convective heat transfer coefficient (W/m^2·K).
T_ambient = 25.  # Ambient temperature (°C).
maskSurfaces = mesh.exteriorFaces

phi = CellVariable(name="solution variable", mesh=mesh, value=T_ambient)

# Set the boundary conditions
# phi.constrain(0., mesh.exteriorFaces)

# get cell centers and face centers
xc, yc, zc = mesh.cellCenters
x, y, z = mesh.faceCenters

# Define the spatially varying diffusion coefficient
# D was a CellVariable previously when not using Robin boundary condition
# D = CellVariable(mesh=mesh, value=0.9)
D = FaceVariable(mesh=mesh, value=0.9)
D.setValue(400., where=z >= Lz/2)

# if the diffusion coefficient is a constant
# D = 1.

# Set the source term for the 5x5x5 cube in the center of the top part
source = CellVariable(mesh=mesh, value=0.0)
mask = ((xc > 7.5) & (xc < 12.5) &
        (yc > 7.5) & (yc < 12.5) &
        (zc > 12.5) & (zc < 17.5))
source.setValue(10., where=mask)

# Robin boundary condition part
# Gamma = FaceVariable(mesh=mesh, value=D)
Gamma = D
Gamma.setValue(0., where=maskSurfaces)
dPf = FaceVariable(mesh=mesh,
                   value=mesh._faceToCellDistanceRatio * mesh.cellDistanceVectors)

n = mesh.faceNormals

a = FaceVariable(mesh=mesh, value=h*n, rank=1)
# b = FaceVariable(mesh=mesh, value=D, rank=0)
b = D
g = FaceVariable(mesh=mesh, value=h*T_ambient, rank=0)

RobinCoeff = maskSurfaces * D * n / (-dPf.dot(a) + b)

# equation with source term
eq = DiffusionTerm(coeff=Gamma) + source + (RobinCoeff * g).divergence - ImplicitSourceTerm(coeff=(RobinCoeff * a.dot(n)).divergence)

# equation without source term
# eq = DiffusionTerm(coeff=Gamma) + (RobinCoeff * g).divergence - ImplicitSourceTerm(coeff=(RobinCoeff * a.dot(n)).divergence)

# start the timer
start_time = time.process_time()
start_time_e = time.perf_counter()

# Solve the equation by sweep()
# timeStepDuration = 0.9 * dx**2 / (2 * D) *10
# res = 1e+10
# while res > 1e-6:
#     res = eq.sweep(var=phi, dt=timeStepDuration)
#     print(res)

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

# end the timer
end_time_e = time.perf_counter()
end_time = time.process_time()
elapsed_time = end_time_e - start_time_e
cpu_time = end_time - start_time
print(f"CPU time: {cpu_time:.4f}s; elapsed time: {elapsed_time:.4f}s.")

if __name__ == '__main__':
    # plot
    viewer = Viewer(vars=phi, datamin=0., datamax=40.)
    viewer.plot()

    # 2D slice at x = Lx/2
    mask_slice = abs(xc - Lx / 2) < dx / 2
    phi_slice = phi.value[mask_slice]
    y_slice = yc[mask_slice]
    z_slice = zc[mask_slice]
    print(np.nanmin(phi_slice), np.nanmax(phi_slice))
    # scatter plot
    plt.figure()
    plt.scatter(y_slice, z_slice, c=phi_slice, cmap='viridis', s=15)
    plt.colorbar(label='Phi Value')
    plt.xlabel('Y')
    plt.ylabel('Z')
    plt.title(f'Slice at x = {Lx / 2}')

    plt.show()

(A side question: do you think it is easier for you if I attach the code file instead of pasting the code here?)

[guyer]

Troubleshooting this with relaxation [*], I was able to identify a sign error in our documentation of the Robin condition (#967).

diff --git a/cube_geos/source3D.py b/cube_geos/source3D.py
index 539c61f..ede517e 100644
--- a/cube_geos/source3D.py
+++ b/cube_geos/source3D.py
@@ -75,7 +74,7 @@ a = FaceVariable(mesh=mesh, value=h*n, rank=1)
 b = D
 g = FaceVariable(mesh=mesh, value=h*T_ambient, rank=0)
 
-RobinCoeff = maskSurfaces * D * n / (-dPf.dot(a) + b)
+RobinCoeff = maskSurfaces * D * n / (dPf.dot(a) + b)
 
 # equation with source term

I've submitted a pull request to fix the documentation (#968).

Additionally, there's an issue with how you've defined Gamma. By writing

Gamma = D

and

b = D

you end up with three Python variables, D, Gamma, and b, that all point to the same thing, so when you write

Gamma.setValue(0., where=maskSurfaces)

all three of them end up zero on the boundary and the Robin condition doesn't work.

Instead, you can do this to zero out the diffusion coefficient on the boundary:

diff --git a/cube_geos/source3D.py b/cube_geos/source3D.py
index 539c61f..ede517e 100644
--- a/cube_geos/source3D.py
+++ b/cube_geos/source3D.py
@@ -63,8 +63,7 @@ source.setValue(10., where=mask)
 
 # Robin boundary condition part
 # Gamma = FaceVariable(mesh=mesh, value=D)
-Gamma = D
-Gamma.setValue(0., where=maskSurfaces)
+Gamma = D * ~maskSurfaces
 dPf = FaceVariable(mesh=mesh,
                    value=mesh._faceToCellDistanceRatio * mesh.cellDistanceVectors)

[*] In general, I always recommend introducing time dependence when steady state solutions don't make sense. Sometimes they're never going to converge without help, but even when they can (like here), when there are errors, the steady state solutions are just nonsense.

[guyer]

These are the solutions I get

[CalebDmArcher]

Thank you so much. It works well!

I have two more follow up question:

1. Previously, you mentioned the equation (eq) does not have non-linear term. Is it because no coefficient depends on the value of the solution variable?
2. Seems like you used 'petsc' as the solver and I am using 'scipy'. Could you show me where I could find the info about the solvers and how to switch the solver manually? I am curious about their accuracies and speeds.

[guyer]

1. Previously, you mentioned the equation (eq) does not have non-linear term. Is it because no coefficient depends on the value of the solution variable?

Correct

2. Seems like you used 'petsc' as the solver and I am using 'scipy'. Could you show me where I could find the info about the solvers and how to switch the solver manually? I am curious about their accuracies and speeds.

I tried both, ultimately. When I used scipy, before I found the sign error, I got the same failure modes that you did, which was helpful to know.

There is some information about the different solvers that FiPy can use here: https://pages.nist.gov/fipy/en/latest/SOLVERS.html

I'm am currently working on a significant refactoring of the solvers to try to harmonize how things like "tolerance" are interpreted.

Short answer, though, is that:

• in serial, PySparse tends to be quite a bit faster than the other solvers. But it requires Python 2.7
• PETSc and Trilinos are comparable in performance, and one or the other is required for parallel solves. PETSc is less of a headache to install.
• SciPy is by far the slowest, but it's pretty easy to get.
• I think it's unlikely that any accuracy differences between them are something you should worry about.

I have some insights into why PySparse is faster when we use it and I hope to be able to bring PETSc and Trilinos on par. Making SciPy faster is beyond my ken.

[CalebDmArcher]

Hi, thank you for your previous suggestions of PETSc.
After some tries, I finally installed PETSc on my Windows PC by using WSL with Ubuntu, and seems Fipy reorganized PETSc as a solver and the speed is indeed faster than before.

But I still have a question: on the page of PETSC, it mentions 'FiPy does not implement any precoditioner objects for PETSc. Simply pass one of the PCType strings in the precon= argument when declaring the solver.' I am not very sure what it says. Is this step mandatory?

Also, would really appreciate it if you could let me know how to switch solvers. I think now Fipy just automatically uses PETSc.

[guyer]

Setting a preconditioner is not mandatory. Using a preconditioner might help with convergence.

Solver selection is discussed at the beginning of https://pages.nist.gov/fipy/en/latest/SOLVERS.html.

[CalebDmArcher]

After another painful PETSc installation process, I finally successfully ran the code with PETSc on my PC.
And You are right, it is much faster.

PETSc 3.18.4 on Ubuntu on Win11 costs:
CPU time: 17.7270s; elapsed time: 2.9414s for GMRES bjacobi precon
CPU time: 17.9074s; elapsed time: 4.0591s for GMRES jacobi precon
CPU time: 17.4173s; elapsed time: 2.7906s for GMRES no precon

I am so happy and cannot imagine how fast it could be if I ran it on a Macbook with M3 chip.

[guyer]

Glad to hear it.

Why painful? conda install --channel conda-forge "petsc==3.18.4"

[CalebDmArcher]

This reply is NOT important.
///////////////////////////////////////////////////////

I used Conda on the server but it took me lots of effort to make it work (as you may know) plus I have less control over the server so I tried PETSc 3.18.4 on my PC locally (Ubuntu on WSL on Win11).
With Ubuntu on WSL on Win11, I did not use Conda. When I installed PETSc 3.20 last time, I followed the command found on their official website official website by using sudo apt install petsc-dev and python -m pip install mpi4py petsc petsc4py (This method directly worked without any wheel cannot be built error). So I thought I would follow a similar path.
ChatGPT said since I wanted to install an archive version, things might get more complex than last time. I followed ChatGPT's suggestion and did the following command which not only does not work at the final pip step but broke fipy (always reported segfaults on multiple memory locations when running my code even with scipy):

sudo apt remove petsc-dev
python -m pip uninstall mpi4py petsc petsc4py
wget https://gitlab.com/petsc/petsc/-/archive/v3.18.4/petsc-v3.18.4.tar.gz
tar -xzf petsc-v3.18.4.tar.gz
cd petsc-v3.18.4
./configure --with-cc=gcc --with-fc=gfortran --with-cxx=g++ --download-fblaslapack --download-mpich --download-metis --download-parmetis --download-superlu
make all
make install
export PETSC_DIR=/path/to/petsc-v3.18.4
export PETSC_ARCH=arch-linux-c-opt
python -m pip install mpi4py petsc petsc4py

After spending certain amount of time, I realized that I might not be able to fix it. So I take ChatGPT's second suggestion which is to use a virtual environment (on Ubuntu on WSL on Windows 11. Environment in an environment in an environment LOL) I honestly have no idea why I did not try Conda at this time. Maybe I was too frustrated and my brain stopped working.

python -m venv fipy-env
source fipy-env/bin/activate
python -m pip install fipy petsc petsc4py mpi4py

which failed to work since fipy cannot detect PETSc. I have tried changing system variables and some other stuff that I cannot recall. But none of them work.

The story ended up with me deleting the virtual environment, uninstalling fipy and PETSc, then re-installing them only using the sudo pip command. And they just suddenly worked nice and fast.

Based on a thread post on PETSc's forum somewhere, if directly using pip command to download PETSc, it will report no wheel can be built error first then solve the problem by itself, and then successfully install. I feel kind of sad when everything could be done in one line of command but I had tons of drama.

[guyer]

I think it's probably safest to always start from the assumption that ChatGPT is an idiot.

[guyer]

As of 7befdbc, run times with PETSc 3.19 and 3.20 should be good again.

[CalebDmArcher]

As of 7befdbc, run times with PETSc 3.19 and 3.20 should be good again.

Is this a fix from PETSc side?
If so it would be fine if using Fipy 3.4.4 with PETSc 3.19/3.20 without manually doing anything, right?

[guyer]

No, the fixes are in FiPy to work around changes in PETSc. At the moment, there is no FiPy release with these changes; you'd need to run source on the master branch.