On 07/11/2018 19:52, Justin Chang wrote:
Okay so I had a slight error in my code. I forgot to transpose the u0 expression by some value such that it won't be negative anywhere (I thought I did it earlier, but I guess I forgot when I did the FEniCS to Firedrake conversion for the MWE). It fixed my issue with CG1, but the IP formulation still does not work.
Having a quick look at your code, I think the penalty term in your IP formulation: + alpha/h_avg*dot(jump(v, n), jump(u_, n))*dS \ needs scaling by D_u. Ideally there should be a factor corresponding to the local min/max ratio of D_u as well. Have a look at http://dx.doi.org/10.1016/j.cam.2006.08.029 for reference. In general, I don't think IP works very well for strongly varied diffusion coefficients, the penalty is either too small and things go unstable or too large and the system becomes very stiff Cheers Stephan