Strong non-linear term treatment

Since Krw and Kro are based on projected quantities (the formulas contain Sw_, which is a projected quantity), they are Function objects, and lambda and fw become Function objects. There is no need to project a Function (unless you want to project it to another FunctionSpace). However, lambda is a keyword in Python so you the code will trigger a SyntaxError. The convention is to use a trailing underscore if you need keywords as variables: lambda_, class_, for_, from_, etc.

Thank you sir. Yes, lambda was called Dlambda, I just replaced it in the text for no good reason. Sorry about that.

Sw (and Sw_) ranges between [0,1] and a plot of Sw and a check of min/max shows this is correct. However, the concurrent fw lies between [-0.12, 1.2] and I cannot figure out a good reason why this would be the case. As the system progresses in time, Sw blows up because fw blows up.

Any ideas on how to force fw to be calculated more accurately? Am I doing something wrong here?

I would grab the array of dofs of Sw_, Krw, Kro and fw and write out the max and min value (.min() and .max() are the numpy class members to call) at every iteration and time step.

The formula looks correct, and if it is really true that Sw_ stays in [0,1], fw should also. It's quite common that Sw_ gets overshoots. It could help to clip the values, i.e., remove all values outside [0,1] in the numpy array and insert back in Sw_ before proceeding.

In the form, fw will be projected, i.e., the combined formula (fw expressed in terms of Sw_) will be projected. Some inaccuracies could appear here, but I would be very surprised if they led to blow up.

Thank you sir. I am trying to add something like:

        fw.array.vector()[np.asarray(np.where( fw.array.vector() < 0 ))[0]]=0
        fw.array.vector()[np.asarray(np.where( fw.array.vector() > 1 ))[0]]=1

in the form function of the PDESubSystem class but fw is of type ufl.algebra.Division (no array.vector() ). I cannot project it either. I tried adding the expression in the update function, I got the following:

*** Error: Unable to access vector of degrees of freedom.
*** Reason: Cannot access a non-const vector from a subfunction.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0

If it helps I made further modifications to the code. The PDESystem form function looks like this now:

class DarcyGlobal(PDESubSystem):
    def form(self, u, u_, u_1, v_u, p, p_, v_p, dt, Sw, Sw_, Sw_1, v_Sw, fw, v_fw, Krw, v_Krw, Kro, v_Kro, **kwargs):

        my_Krw = Sw_*Sw_*Sw_*Sw_
        my_Kro = (1.0-Sw_)*(1.0-Sw_)*(1.0-Sw_*Sw_)

        # some dependent variables, fractional flow
        Dlambda = my_Krw + my_Kro
        my_fw = my_Krw/(my_Krw + my_Kro)

        F = inner(u,v_u)*dx+inner(Dlambda*grad(p),v_u)*dx \
           + inner(u,grad(v_p))*dx \
           + (Sw-Sw_1)*v_Sw/dt*dx + div(fw*u)*v_Sw*dx \
           + fw*v_fw*dx-my_fw*v_fw*dx \
           + Krw*v_Krw*dx-my_Krw*v_Krw*dx \
           + Kro*v_Kro*dx-my_Kro*v_Kro*dx

        return F

and I plot Krw, Kro and fw as I added them as functions in the PDESystem.

The problem is (of course) that fw, Kwo, Kwr are no Function objects, but UFL expressions. Therefore, you cannot extract an array from fw and clip it.

BTW, your code is right, but I think it can be as short as fw.array().vector()[fw.array().vector() > 1] = 1 etc.

But before doing this, you need to either project fw or interpolate it and use the resulting Function in the forms. I'm not sure if interpolation of a UFL expression is as smart as just evaluating the expression at the nodes (will check), but that is the safest computation as your formula will be applied directly at all the nodes. With project you do not get the formula exactly, but some least-squares approximation.

Anyway, you should monitor the evolution of the range of fw, Krw, and Kro before starting to clip fw values. There should only be small overshoots (say max 30%) in two-phase porous media flow. However, this is true only if you use some kind of stabilization for hyperbolic problems. I cannot see any stabilization terms in your forms. A pure Galerkin method for a hyperbolic PDE will give oscillations, and these can for coarse meshes and large time steps completely destroy the solution. This artifact of the Galerkin might be the reason why fw blows up. Fine meshes and small time steps help. Monitor the range of fw in that case.

Thank you very much sir. Just to clarify, when you say stabilization terms you mean something like SUPG, correct?
Do you think I should go for a DG implementation sir or switch to SUPG?

I meant SUPG - and I think that is the easiest first step. Just some v + a*velocity*grad(v) type of weighting function.
But use fine mesh and small dt first to make sure that the implementation works. You don't need SUPG if the mesh is sufficiently fine.

Thank you very much sir.