Solving two equations with a coupledSolver

[FoamScience]

This issue is not directly related to one of your solvers, but it's related to the interface to PETsc

Hi, First of all, thanks for this marvelous piece of software you're maintaining!

I'm attempting to solve a coupled system (for two scalars T1, T2, to which I already have a working segregated solver):

From what I gather from rheoFoam.C, I should:

autoPtr<coupledSolver> cs (new coupledSolver(...));
cs->insertMesh(mesh);
cs->insertField(T1);
cs->insertField(T2);

// Start time loop

// S1 and S2 are matrices for T2
// Define equations
auto T1Eqn( fvm::ddt(C1, T1) + fvc::(phi1G) );              // gravity treated explicitely
auto T2InT1Eqn ( - fvm::laplacian(-M1, T2)  + S1 );
cs->insertEquation(T1.name(), T1.name(), T1Eqn);
cs->insertEquation(T1.name(), T2.name(), T2InT1Eqn);

// Define T2Eqn and T1InT2Eqn similarly and insert them into the matrix
auto T2Eqn( fvm::laplacian(-M2, T2) + fvc::(phi2G) + S2);
auto T1InT2Eqn ( - fvm::ddt(C2, T1)  );
cs->insertEquation(T2.name(), T2.name(), T1Eqn);
cs->insertEquation(T2.name(), T1.name(), T1InT2Eqn);

// Then
cs->solve();

// End time loop

But this doesn't really work; So I decided to check here if the interface is general enough to support such circumstances.

1. I'm not sure about coupling ddt terms in this way (Should it "just work"?)
2. My (working) segregated solver actually solves Eq1 then (Eq1+Eq2) but I don't think it's relevant

I'm new to PETsc. Appologies if this seems trivial

[fppimenta]

The time derivative of the 2nd eq seems a little strange. Should it be in T2? Assuming that the equations are correctly typed, there are wrong signs in the eqs inserted. Are there additional coupling terms in sources S1_2?

[FoamScience]

Thanks for the super-quick reply.

1. The time derivative of the 2nd eq is indeed in "T1"
2. No additional coupling terms in S1 and S2. They are matrices for T2

Should S1 in T2InT1Eqn have a negative sign? auto T2InT1Eqn ( - fvm::laplacian(-M1, T2) - S1 );? or should it be auto T2InT1Eqn ( fvm::laplacian(-M1, T2) + S1 );

I was under the impression I need to insert the opposite of the coupled terms.

[fppimenta]

When the equations are equated to 0 on the same side, they should be inserted as they are written, as in your case.

I was asking about the source terms, because you may achieve additional implicit coupling through them.

[FoamScience]

Thanks for clarifying things; I must have messed with something else then. I'll have to debug the coupled system :)

Yes, I know about coupling the source terms (some experience from foam-extend), in my case though, they both depend on T2 only but I would love to find examples on how to achieve implicit coupling through them (would be useful for future equations). The rheo* solvers don't seem to present examples.

[fppimenta]

The Nernst-Planck equations are the closer I can see in rheoTool from your case. Take a look at this (you will find theory in the user-guide; you can also check the signs of these equations):

https://github.com/fppimenta/rheoTool/tree/master/of70/src/libs/EDFModels/models/NernstPlanck/NernstPlanckCoupled

There are additional cases of coupling, but they are more complex and related with new operators. The theory is here: https://www.sciencedirect.com/science/article/abs/pii/S0045793019302427 and you can find the code in the repo (look inside classes, not only inside solvers).

Re-open issue if needed.