Boundary constrain leads to internal loop, what happens here and why?

Question

My Script ceases to work for my easy conduction problem. Could somebody explain to me why the following line of code T.faceValue.constrain(alp/lam*(Tu-T.faceValue),where=mesh.exteriorFaces) # Boundary Condition for Solver results in fipy giving up on me?

Full Code:

cv=900.
lam=5.
alp=300.

T0 = 25.
Tu = 400.


cellSize = 0.05
radius = 1.


mesh = Gmsh2D('''
              cellSize = %(cellSize)g;
              radius = %(radius)g;
              Point(1) = {0, 0, 0, cellSize};
              Point(2) = {-radius, 0, 0, cellSize};
              Point(3) = {0, radius, 0, cellSize};
              Point(4) = {radius, 0, 0, cellSize};
              Point(5) = {0, -radius, 0, cellSize};
              Circle(6) = {2, 1, 3};
              Circle(7) = {3, 1, 4};
              Circle(8) = {4, 1, 5};
              Circle(9) = {5, 1, 2};
              Line Loop(10) = {6, 7, 8, 9};
              Plane Surface(11) = {10};
              ''' % locals()) # doctest: +GMSH

T = CellVariable(name = "HeatingUp",mesh = mesh,value = T0)


viewer = None
if __name__ == '__main__':
     try:
         viewer = Viewer(vars=T, datamin=T0, datamax=Tu)
         viewer.plotMesh()
         input("Irregular circular mesh. Press <return> to proceed") # doctest: +GMSH
     except:
         print("Unable to create a viewer for an irregular mesh (try Matplotlib2DViewer or MayaviViewer)")
# =============================================================================



eq = TransientTerm(coeff=rho*cv)==DiffusionTerm(coeff=lam)

T.faceValue.constrain(alp/lam*(Tu-T.faceValue),where=mesh.exteriorFaces) # Boundary Condition for Solver



timeStepDuration = 0.1
steps = 10
for step in range(steps):
    eq.solve(var=T, dt=timeStepDuration) # doctest: +GMSH
    if viewer is not None:
        viewer.plot() # doctest: +GMSH
``

score 1 · Accepted Answer · answered Jan 12 '20 at 13:43

You've written that T.faceValue depends on T.faceValue, which depends on T.faceValue, which depends on T.faceValue, ... FiPy has dutifully provided you with the infinite loop that you requested.

Just write T.faceValue.constrain(Tu * (alp/lam) / (1 + alp/lam), where=mesh.exteriorFaces).

In the more likely event that you wanted to relate the gradient to the value at the boundary, please see the discussion on Robin conditions.

Boundary constrain leads to internal loop, what happens here and why?

1 Answers1