You might try dividing through by S(h) and applying the chain rule: (1/S(h)) d/dx(T(h) dh/dx) = d/dx((T(h) / S(h)) dh/dx) + (T(h) / S(h)^2) dh/dx . dS(h)/dx
The second term on the RHS still doesn’t have very efficient representation in FiPy, either, but it’s at least not limited to first order backward Euler. Frankly, I would check that S(h) really sits outside the time derivative. The derivations I’ve seen that look like that usually (and often quietly) assume the coefficient (e.g., heat capacity) is uniform. On May 22, 2020, at 4:41 AM, Urzainqui Aramburu Iñaki (Luke) <inaki.urzain...@luke.fi<mailto:inaki.urzain...@luke.fi>> wrote: Dear developers, I am trying to solve the following diffusion PDE with fipy: S(h) * dh/dt = d/dx(T(h)*dh/dx) + P S and T are both nonlinear functions of h. (There is no closed analytical form to express them, but they could be parameterized somehow, e.g., S(h) = a*exp(b*h)) Reading the docs and following the examples, it seems to me that TransientTerm(coeff=S) would not do the job in FiPy, since this corresponds to d(S*h)/dt. One idea is to write the left hand side of the equation above as: S*dh/dt = d(S*h)/dt – h*dS/dt. Then, the equation would look like so: d(S*h)/dt = d/dx(T(h)*dh/dx) + P + h*dS/dt And this new version of the equation would be implemented as follows: TransientTerm(coeff=S) == DiffusionTerm(coeff=T) + P + ImplicitSourceTerm((S – S.old)/dt) However, this new term makes the equation very unstable. Is this the right way to approach this problem in FiPy? Are there any alternatives? Thanks for your help. Best regards, Iñaki _______________________________________________ fipy mailing list fipy@nist.gov<mailto:fipy@nist.gov> http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
_______________________________________________ fipy mailing list fipy@nist.gov http://www.ctcms.nist.gov/fipy [ NIST internal ONLY: https://email.nist.gov/mailman/listinfo/fipy ]
