Robin boundary conditions have issues

[guyer]

As discussed on the mailing list, there seem to be some problems with Robin boundary conditions.

• contraints don't seem to automatically update, e.g.,

  phi.faceGrad.constrain(-hh * (phi.faceValue - T0)/k, mesh.facesRight)
  while res > restol:
   res = eq.sweep(var=phi)

  behaves differently from

  while res > restol:
   phi.faceGrad.constrain(-hh * (phi.faceValue - T0)/k, mesh.facesRight)
   res = eq.sweep(var=phi)

• underRelaxation may be necessary and should be illustrated, if so
• Robin example robin.py seems to have issues:
  • C.faceGrad.constrain([-P + P * C.faceValue], mesh.facesLeft) will cast away Variableness because a list containing a FiPy Variable is not a Variable. I think this should be C.faceGrad.constrain((-P + P * C.faceValue) * mesh.faceNormals, mesh.facesLeft) to properly assign the desired scalar value to the normal component of the gradient vector. Unfortunately, this change causes the numerical and analytical solutions to no longer agree (see below)
  • The analytical solution does not appear to be a solution to the stated equation and boundary conditions.

[raybsmith]

Just to add, I also agree that the analytical solution in the Robin example is incorrect. I'll try to submit a pull request including the correct analytical solution in the example soon.

[wd15]

Just to note that the analytical solution and the numerical solution seem to be in quite good agreement on the "develop" branch.

Regarding the constraints, I beginning to think that using source terms is clearer and easier. The problem with the constraints is that they really are hard to debug and understand. I want to discuss this further before making any changes.

[wd15]

This is the way to do the Robin boundary condition implicitly.

http://nbviewer.ipython.org/gist/wd15/f04b22c6ec80b1eda8e3

[nlooije]

@guyer, @wd15: What is the status of this bug?

I am running into this problem in a 1D diffusion problem with a reacting boundary condition. Updating the boundary condition at every sweep drives the system to the correct solution but it is taking much longer than i feel is necessary.

[wd15]

@nlooije, thanks for the heads up that this is an issue for you. I'm not sure when we can address it right now.

[nlooije]

@wd15 - I was able to solve the problem using the method you proposed. I have scaled this up to a 2D convection-diffusion problem with reacting boundary conditions and it is working nicely. No longer updating the boundary condition at every sweep has sped up the simulation considerably. Thx!

I have attached the code in case anyone is interested, besides the reactive robin BC it also deals with spatially dependent convection coefficients because of a velocity profile and changing the viewer settings like linestyle and use of markers, etc.

[wd15]

@nlooije, thanks for posting your code. If we ever get around to addressing this problem then I'm sure it will be useful.

[scbarton]

I've run @nlooije's code and I think it has an issue. A spot check of the results indicate that, at the Robin boundary, the face value and the cell center value are equal, which should not be the case in the presence of a gradient. To illustrate, attached is a plot of the solution at x=4.05.

[romarro]

@scbarton this post explains how to make a 1st order general boundary condition: https://stackoverflow.com/questions/42769749/general-boundary-conditions

[aaaaaacd]

Can you show me the code about Robin.py by using the method @nlooije ,the code I make seems to have some issues. `import fipy as fp nx=100. Lx=1. mesh=fp.Grid1D(dx=0.01,nx=100.) C=fp.CellVariable(mesh=mesh,name='C') A=1.0 D=2.0 P=3.0

C.faceGrad.constrain(0,mesh.facesRight) DiffFace=fp.FaceVariable(mesh=mesh,value=A) DiffFace.setValue(0.,where=mesh.facesLeft)

x = mesh.faceCenters convValue = (P,0) Convface=fp.FaceVariable(mesh=mesh,value=convValue) Convface.setValue(0.,where=mesh.facesLeft) mask=(mesh.facesLeftmesh.faceNormals).divergence Af=mesh._cellAreas[mesh.facesLeft][0] dPR=mesh._cellDistances[mesh.facesLeft][0] mask=maskAAf/(1+dPR) eq=0==fp.DiffusionTerm(coeff=DiffFace,var=C)-\ fp.ExponentialConvectionTerm(Convface,var=C)-\ fp.ImplicitSourceTerm(D,var=C)+fp.ImplicitSourceTerm(coeff=Pmask,var=C)-P*mask res=1 restol=1e-4 viewer=fp.Viewer(vars=C) while res>restol: res=eq.sweep(var=C) print('residual: %s' % res) viewer.plot()`

[aaaaaacd]

@wd15 Can you fix the problem above?I don't know where goes wrong.And by the way, Are there more exemplary codes or book about fipy(best in heat transfer)? thanks a lot.

[aaaaaacd]

@jeguyer Can you fix the problem above?I don't know where goes wrong.And by the way, Are there more exemplary codes or book about fipy(best in heat transfer)? thanks a lot.