Hi Daniel, I updated your code a bit and it seems to do better. See https://pastebin.com/51xfqWH3 for the full version.
I changed the "eqX" equation to param = k_c * (P - Peq) maxX = 0.9 largeValue = 1e9 constraint = largeValue * (X > maxX) big_value = 1e9 eqX = TransientTerm(var = X) == param - ImplicitSourceTerm(param, var=X) + maxX * constraint - ImplicitSourceTerm(constraint, var=X) to have a constraint on the maximum value of X. I'm not sure whether this is physical, but it does constrain the max value of X. I changed the "eqY" equation to be, eqT = TransientTerm(rho_b_MH * Cp_b, var=T) == DiffusionTerm(coeff = L_eff, var=T) + rho_b_MH*(dH/M)*rc to have the transient coefficient inside the term. Again, not sure if that matters, but did it anyway. I changed the constraint to be constraint_value = FaceVariable(mesh=mesh) T.faceGrad.constrain(constraint_value,mesh.facesLeft) and then updated the constraint in the loop with for sweep in range(sweeps): constraint_value[:] = h * ((60. + 273.15) - T.faceValue) If FiPy worked correctly that wouldn't be necessary, but evidently it doesn't. Anyway that all seems to work at least as far as I can tell. Cheers, Daniel -- Daniel Wheeler _______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]