fvMatrix, fvOptions and SuSp: automatic implicit/explicit source-term treatment

Zeppo · December 31, 2016, 11:38

fvMatrix, fvOptions and SuSp: automatic implicit/explicit source-term treatment to facilitate the matrix diagonal dominance

The matrix equation to be solved in scalarTransportFoam reads:

Code:

//“scalarTransportFoam.C”
…
61             solve
62             (
63                 fvm::ddt(T)
64               + fvm::div(phi, T)
65               - fvm::laplacian(DT, T)
66              ==
67                 fvOptions(T)
68              );
…

Every term fvm::ddt(T), fvm::div(phi, T) and fvm::laplacian(DT, T) on the left hand side of the equation (before the sign “==”) is a fvMatrix object. Summation of the first and the second terms results in another fvMatrix from which the third term is subtracted resulting in the final fvMatrix object.

The object fvOptions on the right hand side is of IOoptionList type (which is derived from optionList). The class IOoptionList was designed to facilitate the addition of source terms to the matrix equation. It can handle multiple sources of different nature in a systematic way. fvOptions(T) expands into

Code:

//”fvOptionListTemplates.C”
38 template<class Type>
39 Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator()
40 (
41     GeometricField<Type, fvPatchField, volMesh>& field,
42     const word& fieldName
43 )
44 {
45     checkApplied();
46
47     const dimensionSet ds = field.dimensions()/dimTime*dimVolume;
48 
49     tmp<fvMatrix<Type> > tmtx(new fvMatrix<Type>(field, ds));
50     fvMatrix<Type>& mtx = tmtx();
51 
52     forAll(*this, i)
53     {
54         option& source = this->operator[](i);
…
70                 source.addSup(mtx, fieldI);
…
73     }
74 
75     return tmtx;
76 }

Basically it loops through the set of sources and calls addSup function-member for every source:

Code:

70                 source.addSup(mtx, fieldI);

A newly created fvMatrix is passed to the function addSup as an input-output parameter. Every source is responsible for making changes to the fvMatrix depending on what the source was supposed to serve for. After fvMatrix has made the rounds of all the sources, it returns to the place in the code where fvOptions(T) was invoked.

So, initial equation

Code:

   fvm::ddt(T)
 + fvm::div(phi, T)
 - fvm::laplacian(DT, T)
==
   fvOptions(T)

can be rewritten in a form where we have a single fvMatrix object on the left hand side and a single fvMatrix object on the right hand side of the equation (I mark them with letters “A” and “B”):

Code:

fvMatrix_A == fvMatrix_B

This expression is actually a call of the function operator==

Code:

//”fvMatrix.C”
…
1443 template<class Type>
1444 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator==
1445 (
1446     const fvMatrix<Type>& A,
1447     const fvMatrix<Type>& B
1448 )
1449 {
1450     checkMethod(A, B, "==");
1451     return (A - B);
1452 }

(A - B) is calculated from here:

Code:

//”fvMatrix.C”
…
1858 template<class Type>
1859 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator-
1860 (
1861     const fvMatrix<Type>& A,
1862     const fvMatrix<Type>& B
1863 )
1864 {
1865     checkMethod(A, B, "-");
1866     tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A));
1867     tC() -= B;
1868     return tC;
1869 }

tC() -= B is calculated from here:

Code:

//”fvMatrix.C”
..
1059 template<class Type>
1060 void Foam::fvMatrix<Type>::operator-=(const fvMatrix<Type>& fvmv)
1061 {
1062     checkMethod(*this, fvmv, "-=");
1063 
1064     dimensions_ -= fvmv.dimensions_;
1065     lduMatrix::operator-=(fvmv);
1066     source_ -= fvmv.source_;
1067     internalCoeffs_ -= fvmv.internalCoeffs_;
1068     boundaryCoeffs_ -= fvmv.boundaryCoeffs_;
…
1080 }

fvMatrix basically includes 5 ldu-addressing arrays (defined in lduMatrix which fvMatrix is derived from) and a source coefficient array. ldu-addressing arrays are "arrays representing the diagonal, upper, and lower coefficients and the lower and upper indices of the face owner and neighbor, respectively" (c). The source term fvOptions(T) in our equation can only affect the diagonal coefficients array and the source coefficients array. In void Foam::fvMatrix<Type>:

perator-=(const fvMatrix<Type>& fvmv) the diagonal coefficients array is modified here:

Code:

1065     lduMatrix::operator-=(fvmv);

And lduMatrix:

perator-=(const lduMatrix& A) is implemented in this way:

Code:

//"lduMatrixOperations.C"
...
223 void Foam::lduMatrix::operator-=(const lduMatrix& A)
224 {
225     if (A.diagPtr_)
226     {
227         diag() -= A.diag();
228     }
...

The source coefficient array is modified here

Code:

1066     source_ -= fvmv.source_;

All in all, the expression

Code:

fvMatrix_A == fvMatrix_B

is equivalent to creating a new fvMatrix object (I gave it name “C”) defined as

fvMatrix_C = fvMatrix_A - fvMatrix_B

where (in pseudocode):

Code:

C.diag() = A.diag() - B.diag();
C.source() = A.source() - B.source();

Let’s put it all off for a moment and switch over to the specific sources we can use in a matrix equation. There are many types of sources defined in OpenFOAM for a user to choose from. Consider, for example, SemiImplicitSource<Type>. Its function-member addSup looks like this:

Code:

//” SemiImplicitSource.C”
…
133 template<class Type>
134 void Foam::fv::SemiImplicitSource<Type>::addSup
135 (
136     fvMatrix<Type>& eqn,
137     const label fieldI
138 )
139 {
... 
146     const GeometricField<Type, fvPatchField, volMesh>& psi = eqn.psi();
147 
148     DimensionedField<Type, volMesh> Su
149     (
...
159         dimensioned<Type>
160         (
161             "zero",
162             eqn.dimensions()/dimVolume,
163             pTraits<Type>::zero
164         ),
165         false
166     );
167 
168     UIndirectList<Type>(Su, cells_) = injectionRate_[fieldI].first()/VDash_;
169
170     DimensionedField<scalar, volMesh> Sp
171     (
...
181         dimensioned<scalar>
182         (
183             "zero",
184             Su.dimensions()/psi.dimensions(),
185             0.0
186         ),
187         false
188     );
189 
190     UIndirectList<scalar>(Sp, cells_) = injectionRate_[fieldI].second()/VDash_;
191
192     eqn += Su + fvm::SuSp(S	p, psi);
193 }

Su and Sp are actually constants which are read from a user-defined dictionary "system/fvOptions":

Code:

//"system/fvOptions"	
options
{
    energySource_1
    {
        type            scalarSemiImplicitSource;
        active          yes;

        scalarSemiImplicitSourceCoeffs
        {
            timeStart       0;
            duration        100000000;
            selectionMode   cellZone; // all; points; cellSet; cellZone;
            cellZone        cellZone_2;
            volumeMode      specific; // specific; absolute
            injectionRateSuSp
            {
                T           (0.0 10.0); //- absolute [W]; specific [w/m^3]
            }
        }
	}
}

In this example Su = 0.0 and Sp = 10.0. We can omit Su from

Code:

eqn += Su + fvm::SuSp(Sp, psi);

and have only

Code:

eqn += fvm::SuSp(Sp, psi);

eqn is a fvMatrix initially (at the point where it was constructed) initialized with zero. Function fvm::SuSp(Sp, psi) returns a fvMatrix too. The second is added to the first with operator+= :

Code:

// ”fvMatrix.C”
..
1025 template<class Type>
1026 void Foam::fvMatrix<Type>::operator+=(const fvMatrix<Type>& fvmv)
1027 {
1028     checkMethod(*this, fvmv, "+=");
1029 
1030     dimensions_ += fvmv.dimensions_;
1031     lduMatrix::operator+=(fvmv);
1032     source_ += fvmv.source_;
1033     internalCoeffs_ += fvmv.internalCoeffs_;
1034     boundaryCoeffs_ += fvmv.boundaryCoeffs_;
...

which leads to (in pseudocode):

Code:

eqn.diag() += (fvm::SuSp(Sp, psi)).diag();
eqn.source() += (fvm::SuSp(Sp, psi)).source()

The definition of function fvm::SuSp(Sp, psi) can be found in ”fvmSup.C”:

Code:

// ”fvmSup.C”
190 template<class Type>
191 Foam::tmp<Foam::fvMatrix<Type> >
192 Foam::fvm::SuSp
193 (
194     const DimensionedField<scalar, volMesh>& susp,
195     const GeometricField<Type, fvPatchField, volMesh>& vf
196 )
197 {
198     const fvMesh& mesh = vf.mesh();
199 
200     tmp<fvMatrix<Type> > tfvm
201     (
202         new fvMatrix<Type>
203         (
204             vf,
205             dimVol*susp.dimensions()*vf.dimensions()
206         )
207     );
208     fvMatrix<Type>& fvm = tfvm();
209 
210     fvm.diag() += mesh.V()*max(susp.field(), scalar(0));
211 
212     fvm.source() -= mesh.V()*min(susp.field(), scalar(0))
213         *vf.internalField();
214 
215     return tfvm;
216 }

Its role is to perform the "automatic implicit/explicit source-term treatment to facilitate the matrix diagonal dominance" (c). When Sp is positive, it is treated as implicit, modifying the diagonal coefficients in fvMatrix:

Code:

210     fvm.diag() += mesh.V()*max(susp.field(), scalar(0));

For the case when Sp is negative, it is treated as explicit, modifying the source coefficients in fvMatrix:

Code:

212     fvm.source() -= mesh.V()*min(susp.field(), scalar(0))
213         *vf.internalField();

All in all, omitting irrelevant details, the whole code can be summed up in a rather compact piece of pseudocode:

Code:

{
    //- transient+convection+diffusion matrix
    fvMatrix<scalar> tmA = 
      fvm::ddt(T)
    + fvm::div(phi, T)
    - fvm::laplacian(DT, T);
    
    //- uniform implicit source.
    //- Formally defined here, but actually read from a dictionary-file
    DimensionedField Sp(..., dimensionedScalar(..., 10.0))
    
    //- matrix defined in fvm::SuSp(…)
    fvMatrix<scalar> tmSp (...);
    tmSp.diag() += mesh.V()*Sp.field();
    
    //- matrix defined in optionList::operator()
    fvMatrix<scalar> tmOL (...);
    tmOL.diag() += tmSp.diag();
    
    tmp<fvMatrix<scalar> > tmC (...);
    tmC().diag() = tmA.diag() - tmOL.diag();
    
    solve(tmC);
}

The line

Code:

tmC().diag() = tmA.diag() - tmOL.diag();

expands into

Code:

tmC().diag() = tmA.diag() - mesh.V()*Sp.field();

As you can see the positive term mesh.V()*Sp.field() is subtracted from the diagonal coefficients of the matrix. This means that the diagonal dominance deteriorates!

Is it a bug? If it is not, can anyone tell me where I went wrong in my reasoning?

(Kind of) solution

The problem comes from the fact that the initial equation is formulated in the form:

Code:

fvMatrix_A == fvOptions(T);

Code:

fvMatrix_A == SuSp(Sp, T);

which for the positive Sp results in a new fvMatrix_C with the diagonal coefficients being defined:

Code:

fvMatrix_C.diag() = fvMatrix_A.diag() - mesh.V()*Sp.field;

If we reformulate it, carrying the term fvOptions(T) from the rhs to the lhs and bringing the sign “-“ into the brackets SuSp(-Sp, psi) we have

Code:

fvMatrix_A == fvOptions(T);

Code:

fvMatrix_A == SuSp(Sp, T);

Code:

fvMatrix_A - SuSp(Sp, T) == 0;

Code:

fvMatrix_A + SuSp(-Sp, T) == 0;

which for the positive Sp results in a fvMatrix_C with the source coefficients being defined:

Code:

fvMatrix_C.source() = 
fvMatrix_A.source() + -(mesh.V()*(-Sp.field())*T.internalField());

or

Code:

fvMatrix_C.source() = 
fvMatrix_A.source() + mesh.V()*Sp.field()*T.internalField();

The diagonal coefficients are left intact (to not deteriorate the diagonal dominance) and only the source coefficients are altered. So, the reformulated equation

Code:

fvMatrix_A + SuSp(-Sp, T) == 0;

is discretized in a perfectly correct way.

The problem is that solvers in OpenFOAM are designed so that fvOptions is on the rhs which looks quite reasonable from the math base:

Code:

fvMatrix_A == fvOptions(T);

But, as you just saw, it leads to the “incorrect” result. This “incorrectness” isn’t a severe error. It just makes the matix less robust and sometimes may lead to divergence in matrix solition.

Thanks.

Zeppo · February 23, 2017, 09:05

The final part of the previous post (called "(Kind of) solution") was incorrect and should be ignored. The info presented below is what should have actually went in place of that wrong part.

(Kind of) solution
The problem comes from the fact that the choice between implicit/explicit formulation of a source term SuSp(Sp, T) depends not only on the Sp being positive/negative but also on the way this source term interacts with other terms in the matrix equation. For example,

Code:

1) fvMatrix_A + SuSp(Sp, T);

and

Code:

2) fvMatrix_A - SuSp(Sp, T);

should be treated differently. If Sp>0 then SuSp(Sp, T) should be discretized implicitly for the first expression and explicitly for the second one. If Sp<0 then things should be the other way around: explicit discretization for the first expression and implicit discretization for the second one. Function SuSp() just doesn’t know how fvMatrix it produces is going to be used later on. The client code (code which calls SuSp()) does know!

The second expression should be rewritten by bringing the sign “-“ into the brackets (i.e. effectively, multiplying Sp by -1: SuSp(-Sp, T)) :

Code:

2) fvMatrix_A + SuSp(-Sp, T);

Now it is correct: for positive Sp the term SuSp(-Sp, T) is treated implicitly and for negative Sp - explicitly.

Most matrix equations in OpenFOAM's solvers look like:

Code:

fvMatrix_A == fvOptions(T);

which for a source of type scalarSemiImplicitSource expands into

Code:

fvMatrix_A == SuSp(Sp, T);

which (as you just saw) leads to the incorrect behaviour of the code and should be rewritten in this way:

Code:

fvMatrix_A == -SuSp(-Sp, T);

According to the “study” presented above the line of code in SemiImplicitSource.C should be rewritten from

Code:

// ”SemiImplicitSource.C”
…
133 template<class Type>
134 void Foam::fv::SemiImplicitSource<Type>::addSup
135 (
136     fvMatrix<Type>& eqn,
137     const label fieldI
138 )
139 {
... 
192     eqn += Su + fvm::SuSp(Sp, psi);
193 }

to

Code:

// ”SemiImplicitSource.C”
…
133 template<class Type>
134 void Foam::fv::SemiImplicitSource<Type>::addSup
135 (
136     fvMatrix<Type>& eqn,
137     const label fieldI
138 )
139 {
... 
192     eqn += Su - fvm::SuSp(-Sp, psi);
193 }

Now it looks correct to me.

P.S. I believe something like lazy evaluation technique can also solve this sort of problems. In which case the actual calculation of a final fvMatrix can be deffered until the very last operation in the series of matrix operations. Unfortunately, within a frame of current OpenFOAM implementation this is not yet posible.

gtarang · April 9, 2019, 22:44

Thanks a lot!

Today I also came across this issue and was really hoping for some solution. You have nicely describe the problem. Any solution or explanation you found? If not I will also put a problem from solidificationMeltingSource.C where it calculated Sp (sink) term for velocity, and I have found that eqn.Sp() is negative at the time when it should provide sink to the velocity equation.

Regards
gtarang

Zeppo · April 13, 2019, 13:55

Solution is given in my second post:

Quote:

Originally Posted by Zeppo

the line of code in SemiImplicitSource.C should be rewritten from

Code:

// ”SemiImplicitSource.C”
…
133 template<class Type>
134 void Foam::fv::SemiImplicitSource<Type>::addSup
135 (
136     fvMatrix<Type>& eqn,
137     const label fieldI
138 )
139 {
... 
192     eqn += Su + fvm::SuSp(Sp, psi);
193 }

to

Code:

// ”SemiImplicitSource.C”
…
133 template<class Type>
134 void Foam::fv::SemiImplicitSource<Type>::addSup
135 (
136     fvMatrix<Type>& eqn,
137     const label fieldI
138 )
139 {
... 
192     eqn += Su - fvm::SuSp(-Sp, psi);
193 }

Rereading the text I have spotted a mistake in the line

Quote:

Originally Posted by Zeppo

Now it is correct: for positive Sp the term SuSp(-Sp, T) is treated implicitly and for negative Sp - explicitly.

It should read

Quote:

Now it is correct: for positive Sp the term SuSp(-Sp, T) is treated explicitly and for negative Sp - implicitly.

Diro7 · April 23, 2020, 05:15

Hi Sergei,

Thank you for your brilliant post!
I'm relatively new to working with OF source code and I didn't expect to find such kind of issues in my own code.

Basically I hard-coded a semi-implicit temperature source in my incompressible solver.
Everything worked fine until I tried to achieve the same result using fvOptions and scalarSemiImplicitSource.

I struggled for several days, with my solver was repeatedly exploding.
Then I recalled having found this thread and I tried to change the sign to both the fvOptions term in TEqn and the Su/Sp coefficients in the fvOptions dictionary file.
Magically, now everything works just fine as before. It is quite annoying indeed.

As I understand the solution is either to change that single line in semiImplicitSource.C (which is the most correct approach I think, but quite overkill IMHO since I should compile my own version of the template class) or, as I did, change signs in the specific equation and in the corresponding dictionary file (but I admit it is ugly, inconsistent with other equations and would behave badly if other kinds of fvOptions are added to the same equation, I believe).
Another solution could be maybe to replace fvm::SuSp() with fvm::Sp() which should enforce implicit treatment in all cases, if I got it right. I didn't try yet, though.

Anyway, I acknowledge that after some reasoning one can recognize and solve this issue, but I believe that one must have at least some experience with the source code to find out, understand and correct it.
If fvOptions are intended to provide all users (more basic ones included) with some run-time tool to add simple sources to equations (with semiImplicitSource being the simplest and most common one), I believe that such behaviour can be regarded as a bug.
Basically, if you try to use it as explained in the official documentation things won't simply work in many cases with no apparent reason and no possibility of solution without touching the source code.
Did you try to file some kind of official issue report about this?

Andrea

Zeppo · April 23, 2020, 14:46

Hi Andrea.

No, I haven't submitted any bug report on this issue. I believe that the most correct, consistent and robust way to overcome it is to build your own version of SemiImplicitSource class with its function addSup being modified:

Code:

// ”SemiImplicitSource.C”
…
133 template<class Type>
134 void Foam::fv::SemiImplicitSource<Type>::addSup
135 (
136     fvMatrix<Type>& eqn,
137     const label fieldI
138 )
139 {
... 
192     //eqn += Su + fvm::SuSp(Sp, psi); //- comment out this line
          eqn += Su - fvm::SuSp(-Sp, psi); //- add this line
193 }

Tibo99 · May 5, 2020, 13:30

Hi everyone!

I tried the fvOptions with scalarCodedSource to link a source term on epsilon to the k-e turbulence model for 1D (x=1mm, y=150mm, z=1mm) ABL simulation and it seem that I got similar issue.

If you got time and knowledge, It would be appreciated to know your insight/suggestion/tips about the code I did?

By the results, the calculation by the formula are good since I printed out the results in a text file and verified the validity but the source term doesn't seem to be applied to the domain(see picture). The boundary wall seem to work (because in the epsilon patch I choose the BC "epsilonWallFunction") and then after, the values are null. The simulation stop after 11 step of iteration...

Thank again for your comment and your help.

Best Regards,

Code:

/*--------------------------------*- C++ -*----------------------------------*\
| =========                 |                                                 |
| \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox           |
|  \\    /   O peration     | Version:  v1906                                 |
|   \\  /    A nd           | Web:      www.OpenFOAM.com                      |
|    \\/     M anipulation  |                                                 |
\*---------------------------------------------------------------------------*/
FoamFile
{
    version     v1906;
    format      ascii;
    class       dictionary;
    location    "system";
    object      fvOptions;
}
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

epsilonSource
{
    type            		scalarCodedSource;
    selectionMode   		all;
    active			true;

    fields          		(epsilon);
    name			epsilonSource;
	
    
    codeInclude
    	#{
		#include "OFstream.H"
	#};
    
    codeAddSup
    	#{
	const Time& time = mesh_.time();
				
        const scalar us_ = 0.76;
        const scalar D1k_ = 1.44;
        const scalar D2k_ = 1.92;
        const scalar kappa_ = 0.41;
        const scalar sigmaEps_ = 1.17;
        const scalar Cmu_ = 0.09;
        const scalar z0_ = 3.3474e-05;
        
        const scalarField& y = mesh_.C().component(1);
        
	scalarField& epsilonSource = eqn.source();
				
        // Apply the source
        forAll(y, i)
        {
          
            {
               epsilonSource[i] = (pow(us_,4)/(pow((y[i]+z0_),2)))
               *(((D2k_-D1k_)*sqrt(Cmu_)/pow(kappa_,2))-(1/sigmaEps_));
            }   
                               
			   			
			std::ofstream file;
			file.open ("Epsilonterm.txt", std::ofstream::out | std::ofstream::app);
			file << time.value()<< "	" << epsilonSource[i] << "	" << y[i] << endl << "\n";
			file.close();
						
        
        };
        
        
    #};
    
    codeCorrect
    #{
    #};
    
    codeConstrain
    #{
    #};
}
// ************************************************************************* //

nipinl · December 15, 2021, 11:20

I made a small test case to see the effects of the source term to fvMatrix. In my case (heat conduction, basically laplacian solver) since diag() is negative, the following equation(given in the last line) increases the diagonal dominance, not deteriorating it.

I feel mostly diag() will be having -ve entries and the original implementation can be retained, isn't it?

tmC().diag() = tmA.diag() - mesh.V()*Sp.field();

December 31, 2016, 11:38	fvMatrix, fvOptions and SuSp: automatic implicit/explicit source-term treatment	#1
Zeppo Senior Member Sergei Join Date: Dec 2009 Posts: 261 Rep Power: 21	fvMatrix, fvOptions and SuSp: automatic implicit/explicit source-term treatment to facilitate the matrix diagonal dominance The matrix equation to be solved in scalarTransportFoam reads: Code: //“scalarTransportFoam.C” … 61 solve 62 ( 63 fvm::ddt(T) 64 + fvm::div(phi, T) 65 - fvm::laplacian(DT, T) 66 == 67 fvOptions(T) 68 ); … Every term fvm::ddt(T), fvm::div(phi, T) and fvm::laplacian(DT, T) on the left hand side of the equation (before the sign “==”) is a fvMatrix object. Summation of the first and the second terms results in another fvMatrix from which the third term is subtracted resulting in the final fvMatrix object. The object fvOptions on the right hand side is of IOoptionList type (which is derived from optionList). The class IOoptionList was designed to facilitate the addition of source terms to the matrix equation. It can handle multiple sources of different nature in a systematic way. fvOptions(T) expands into Code: //”fvOptionListTemplates.C” 38 template<class Type> 39 Foam::tmp<Foam::fvMatrix<Type> > Foam::fv::optionList::operator() 40 ( 41 GeometricField<Type, fvPatchField, volMesh>& field, 42 const word& fieldName 43 ) 44 { 45 checkApplied(); 46 47 const dimensionSet ds = field.dimensions()/dimTimedimVolume; 48 49 tmp<fvMatrix<Type> > tmtx(new fvMatrix<Type>(field, ds)); 50 fvMatrix<Type>& mtx = tmtx(); 51 52 forAll(this, i) 53 { 54 option& source = this->operator[](i); … 70 source.addSup(mtx, fieldI); … 73 } 74 75 return tmtx; 76 } Basically it loops through the set of sources and calls addSup function-member for every source: Code: 70 source.addSup(mtx, fieldI); A newly created fvMatrix is passed to the function addSup as an input-output parameter. Every source is responsible for making changes to the fvMatrix depending on what the source was supposed to serve for. After fvMatrix has made the rounds of all the sources, it returns to the place in the code where fvOptions(T) was invoked. So, initial equation Code: fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(DT, T) == fvOptions(T) can be rewritten in a form where we have a single fvMatrix object on the left hand side and a single fvMatrix object on the right hand side of the equation (I mark them with letters “A” and “B”): Code: fvMatrix_A == fvMatrix_B This expression is actually a call of the function operator== Code: //”fvMatrix.C” … 1443 template<class Type> 1444 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator== 1445 ( 1446 const fvMatrix<Type>& A, 1447 const fvMatrix<Type>& B 1448 ) 1449 { 1450 checkMethod(A, B, "=="); 1451 return (A - B); 1452 } (A - B) is calculated from here: Code: //”fvMatrix.C” … 1858 template<class Type> 1859 Foam::tmp<Foam::fvMatrix<Type> > Foam::operator- 1860 ( 1861 const fvMatrix<Type>& A, 1862 const fvMatrix<Type>& B 1863 ) 1864 { 1865 checkMethod(A, B, "-"); 1866 tmp<fvMatrix<Type> > tC(new fvMatrix<Type>(A)); 1867 tC() -= B; 1868 return tC; 1869 } tC() -= B is calculated from here: Code: //”fvMatrix.C” .. 1059 template<class Type> 1060 void Foam::fvMatrix<Type>::operator-=(const fvMatrix<Type>& fvmv) 1061 { 1062 checkMethod(this, fvmv, "-="); 1063 1064 dimensions_ -= fvmv.dimensions_; 1065 lduMatrix::operator-=(fvmv); 1066 source_ -= fvmv.source_; 1067 internalCoeffs_ -= fvmv.internalCoeffs_; 1068 boundaryCoeffs_ -= fvmv.boundaryCoeffs_; … 1080 } fvMatrix basically includes 5 ldu-addressing arrays (defined in lduMatrix which fvMatrix is derived from) and a source coefficient array. ldu-addressing arrays are "arrays representing the diagonal, upper, and lower coefficients and the lower and upper indices of the face owner and neighbor, respectively" (c). The source term fvOptions(T) in our equation can only affect the diagonal coefficients array and the source coefficients array. In void Foam::fvMatrix<Type>:perator-=(const fvMatrix<Type>& fvmv) the diagonal coefficients array is modified here: Code: 1065 lduMatrix::operator-=(fvmv); And lduMatrix:perator-=(const lduMatrix& A) is implemented in this way: Code: //"lduMatrixOperations.C" ... 223 void Foam::lduMatrix::operator-=(const lduMatrix& A) 224 { 225 if (A.diagPtr_) 226 { 227 diag() -= A.diag(); 228 } ... The source coefficient array is modified here Code: 1066 source_ -= fvmv.source_; All in all, the expression Code: fvMatrix_A == fvMatrix_B is equivalent to creating a new fvMatrix object (I gave it name “C”) defined as fvMatrix_C = fvMatrix_A - fvMatrix_B where (in pseudocode): Code: C.diag() = A.diag() - B.diag(); C.source() = A.source() - B.source(); Let’s put it all off for a moment and switch over to the specific sources we can use in a matrix equation. There are many types of sources defined in OpenFOAM for a user to choose from. Consider, for example, SemiImplicitSource<Type>. Its function-member addSup looks like this: Code: //” SemiImplicitSource.C” … 133 template<class Type> 134 void Foam::fv::SemiImplicitSource<Type>::addSup 135 ( 136 fvMatrix<Type>& eqn, 137 const label fieldI 138 ) 139 { ... 146 const GeometricField<Type, fvPatchField, volMesh>& psi = eqn.psi(); 147 148 DimensionedField<Type, volMesh> Su 149 ( ... 159 dimensioned<Type> 160 ( 161 "zero", 162 eqn.dimensions()/dimVolume, 163 pTraits<Type>::zero 164 ), 165 false 166 ); 167 168 UIndirectList<Type>(Su, cells_) = injectionRate_[fieldI].first()/VDash_; 169 170 DimensionedField<scalar, volMesh> Sp 171 ( ... 181 dimensioned<scalar> 182 ( 183 "zero", 184 Su.dimensions()/psi.dimensions(), 185 0.0 186 ), 187 false 188 ); 189 190 UIndirectList<scalar>(Sp, cells_) = injectionRate_[fieldI].second()/VDash_; 191 192 eqn += Su + fvm::SuSp(S p, psi); 193 } Su and Sp are actually constants which are read from a user-defined dictionary "system/fvOptions": Code: //"system/fvOptions" options { energySource_1 { type scalarSemiImplicitSource; active yes; scalarSemiImplicitSourceCoeffs { timeStart 0; duration 100000000; selectionMode cellZone; // all; points; cellSet; cellZone; cellZone cellZone_2; volumeMode specific; // specific; absolute injectionRateSuSp { T (0.0 10.0); //- absolute [W]; specific [w/m^3] } } } } In this example Su = 0.0 and Sp = 10.0. We can omit Su from Code: eqn += Su + fvm::SuSp(Sp, psi); and have only Code: eqn += fvm::SuSp(Sp, psi); eqn is a fvMatrix initially (at the point where it was constructed) initialized with zero. Function fvm::SuSp(Sp, psi) returns a fvMatrix too. The second is added to the first with operator+= : Code: // ”fvMatrix.C” .. 1025 template<class Type> 1026 void Foam::fvMatrix<Type>::operator+=(const fvMatrix<Type>& fvmv) 1027 { 1028 checkMethod(this, fvmv, "+="); 1029 1030 dimensions_ += fvmv.dimensions_; 1031 lduMatrix::operator+=(fvmv); 1032 source_ += fvmv.source_; 1033 internalCoeffs_ += fvmv.internalCoeffs_; 1034 boundaryCoeffs_ += fvmv.boundaryCoeffs_; ... which leads to (in pseudocode): Code: eqn.diag() += (fvm::SuSp(Sp, psi)).diag(); eqn.source() += (fvm::SuSp(Sp, psi)).source() The definition of function fvm::SuSp(Sp, psi) can be found in ”fvmSup.C”: Code: // ”fvmSup.C” 190 template<class Type> 191 Foam::tmp<Foam::fvMatrix<Type> > 192 Foam::fvm::SuSp 193 ( 194 const DimensionedField<scalar, volMesh>& susp, 195 const GeometricField<Type, fvPatchField, volMesh>& vf 196 ) 197 { 198 const fvMesh& mesh = vf.mesh(); 199 200 tmp<fvMatrix<Type> > tfvm 201 ( 202 new fvMatrix<Type> 203 ( 204 vf, 205 dimVolsusp.dimensions()vf.dimensions() 206 ) 207 ); 208 fvMatrix<Type>& fvm = tfvm(); 209 210 fvm.diag() += mesh.V()max(susp.field(), scalar(0)); 211 212 fvm.source() -= mesh.V()min(susp.field(), scalar(0)) 213 vf.internalField(); 214 215 return tfvm; 216 } Its role is to perform the "automatic implicit/explicit source-term treatment to facilitate the matrix diagonal dominance" (c). When Sp is positive, it is treated as implicit, modifying the diagonal coefficients in fvMatrix: Code: 210 fvm.diag() += mesh.V()max(susp.field(), scalar(0)); For the case when Sp is negative, it is treated as explicit, modifying the source coefficients in fvMatrix: Code: 212 fvm.source() -= mesh.V()min(susp.field(), scalar(0)) 213 vf.internalField(); All in all, omitting irrelevant details, the whole code can be summed up in a rather compact piece of pseudocode: Code: { //- transient+convection+diffusion matrix fvMatrix<scalar> tmA = fvm::ddt(T) + fvm::div(phi, T) - fvm::laplacian(DT, T); //- uniform implicit source. //- Formally defined here, but actually read from a dictionary-file DimensionedField Sp(..., dimensionedScalar(..., 10.0)) //- matrix defined in fvm::SuSp(…) fvMatrix<scalar> tmSp (...); tmSp.diag() += mesh.V()Sp.field(); //- matrix defined in optionList::operator() fvMatrix<scalar> tmOL (...); tmOL.diag() += tmSp.diag(); tmp<fvMatrix<scalar> > tmC (...); tmC().diag() = tmA.diag() - tmOL.diag(); solve(tmC); } The line Code: tmC().diag() = tmA.diag() - tmOL.diag(); expands into Code: tmC().diag() = tmA.diag() - mesh.V()Sp.field(); As you can see the positive term mesh.V()Sp.field() is subtracted from the diagonal coefficients of the matrix. This means that the diagonal dominance deteriorates! Is it a bug? If it is not, can anyone tell me where I went wrong in my reasoning? (Kind of) solution The problem comes from the fact that the initial equation is formulated in the form: Code: fvMatrix_A == fvOptions(T); Code: fvMatrix_A == SuSp(Sp, T); which for the positive Sp results in a new fvMatrix_C with the diagonal coefficients being defined: Code: fvMatrix_C.diag() = fvMatrix_A.diag() - mesh.V()Sp.field; If we reformulate it, carrying the term fvOptions(T) from the rhs to the lhs and bringing the sign “-“ into the brackets SuSp(-Sp, psi) we have Code: fvMatrix_A == fvOptions(T); Code: fvMatrix_A == SuSp(Sp, T); Code: fvMatrix_A - SuSp(Sp, T) == 0; Code: fvMatrix_A + SuSp(-Sp, T) == 0; which for the positive Sp results in a fvMatrix_C with the source coefficients being defined: Code: fvMatrix_C.source() = fvMatrix_A.source() + -(mesh.V()(-Sp.field())T.internalField()); or Code: fvMatrix_C.source() = fvMatrix_A.source() + mesh.V()Sp.field()T.internalField(); The diagonal coefficients are left intact (to not deteriorate the diagonal dominance) and only the source coefficients are altered. So, the reformulated equation Code: fvMatrix_A + SuSp(-Sp, T) == 0; is discretized in a perfectly correct way. The problem is that solvers in OpenFOAM are designed so that fvOptions is on the rhs which looks quite reasonable from the math base: Code: fvMatrix_A == fvOptions(T); But, as you just saw, it leads to the “incorrect” result. This “incorrectness” isn’t a severe error. It just makes the matix less robust and sometimes may lead to divergence in matrix solition. Thanks. kiski, Zaphod'sSecondHead, Daniel_Khazaei and 19 others like this. Last edited by Zeppo; January 2, 2017 at 07:21.

February 23, 2017, 09:05		#2
Zeppo Senior Member Sergei Join Date: Dec 2009 Posts: 261 Rep Power: 21	The final part of the previous post (called "(Kind of) solution") was incorrect and should be ignored. The info presented below is what should have actually went in place of that wrong part. (Kind of) solution The problem comes from the fact that the choice between implicit/explicit formulation of a source term SuSp(Sp, T) depends not only on the Sp being positive/negative but also on the way this source term interacts with other terms in the matrix equation. For example, Code: 1) fvMatrix_A + SuSp(Sp, T); and Code: 2) fvMatrix_A - SuSp(Sp, T); should be treated differently. If Sp>0 then SuSp(Sp, T) should be discretized implicitly for the first expression and explicitly for the second one. If Sp<0 then things should be the other way around: explicit discretization for the first expression and implicit discretization for the second one. Function SuSp() just doesn’t know how fvMatrix it produces is going to be used later on. The client code (code which calls SuSp()) does know! The second expression should be rewritten by bringing the sign “-“ into the brackets (i.e. effectively, multiplying Sp by -1: SuSp(-Sp, T)) : Code: 2) fvMatrix_A + SuSp(-Sp, T); Now it is correct: for positive Sp the term SuSp(-Sp, T) is treated implicitly and for negative Sp - explicitly. Most matrix equations in OpenFOAM's solvers look like: Code: fvMatrix_A == fvOptions(T); which for a source of type scalarSemiImplicitSource expands into Code: fvMatrix_A == SuSp(Sp, T); which (as you just saw) leads to the incorrect behaviour of the code and should be rewritten in this way: Code: fvMatrix_A == -SuSp(-Sp, T); According to the “study” presented above the line of code in SemiImplicitSource.C should be rewritten from Code: // ”SemiImplicitSource.C” … 133 template<class Type> 134 void Foam::fv::SemiImplicitSource<Type>::addSup 135 ( 136 fvMatrix<Type>& eqn, 137 const label fieldI 138 ) 139 { ... 192 eqn += Su + fvm::SuSp(Sp, psi); 193 } to Code: // ”SemiImplicitSource.C” … 133 template<class Type> 134 void Foam::fv::SemiImplicitSource<Type>::addSup 135 ( 136 fvMatrix<Type>& eqn, 137 const label fieldI 138 ) 139 { ... 192 eqn += Su - fvm::SuSp(-Sp, psi); 193 } Now it looks correct to me. P.S. I believe something like lazy evaluation technique can also solve this sort of problems. In which case the actual calculation of a final fvMatrix can be deffered until the very last operation in the series of matrix operations. Unfortunately, within a frame of current OpenFOAM implementation this is not yet posible. randolph, Kummi, DoQuocVu and 6 others like this.

April 23, 2020, 14:46		#6
Zeppo Senior Member Sergei Join Date: Dec 2009 Posts: 261 Rep Power: 21	Hi Andrea. No, I haven't submitted any bug report on this issue. I believe that the most correct, consistent and robust way to overcome it is to build your own version of SemiImplicitSource class with its function addSup being modified: Code: // ”SemiImplicitSource.C” … 133 template<class Type> 134 void Foam::fv::SemiImplicitSource<Type>::addSup 135 ( 136 fvMatrix<Type>& eqn, 137 const label fieldI 138 ) 139 { ... 192 //eqn += Su + fvm::SuSp(Sp, psi); //- comment out this line eqn += Su - fvm::SuSp(-Sp, psi); //- add this line 193 }

December 15, 2021, 11:20		#8
nipinl New Member Nipin L Join Date: Nov 2012 Location: Canada Posts: 23 Rep Power: 13	I made a small test case to see the effects of the source term to fvMatrix. In my case (heat conduction, basically laplacian solver) since diag() is negative, the following equation(given in the last line) increases the diagonal dominance, not deteriorating it. I feel mostly diag() will be having -ve entries and the original implementation can be retained, isn't it? tmC().diag() = tmA.diag() - mesh.V()Sp.field(); Last edited by nipinl; December 15, 2021 at 11:24. Reason: adding link to the test cases*

April 9, 2019, 22:44		#3
gtarang Member Tarang Join Date: Feb 2011 Location: Delhi, India Posts: 47 Rep Power: 15	Thanks a lot! Today I also came across this issue and was really hoping for some solution. You have nicely describe the problem. Any solution or explanation you found? If not I will also put a problem from solidificationMeltingSource.C where it calculated Sp (sink) term for velocity, and I have found that eqn.Sp() is negative at the time when it should provide sink to the velocity equation. Regards gtarang

April 23, 2020, 05:15		#5
Diro7 Member Andrea Di Ronco Join Date: Nov 2016 Location: Milano, Italy Posts: 57 Rep Power: 9	Hi Sergei, Thank you for your brilliant post! I'm relatively new to working with OF source code and I didn't expect to find such kind of issues in my own code. Basically I hard-coded a semi-implicit temperature source in my incompressible solver. Everything worked fine until I tried to achieve the same result using fvOptions and scalarSemiImplicitSource. I struggled for several days, with my solver was repeatedly exploding. Then I recalled having found this thread and I tried to change the sign to both the fvOptions term in TEqn and the Su/Sp coefficients in the fvOptions dictionary file. Magically, now everything works just fine as before. It is quite annoying indeed. As I understand the solution is either to change that single line in semiImplicitSource.C (which is the most correct approach I think, but quite overkill IMHO since I should compile my own version of the template class) or, as I did, change signs in the specific equation and in the corresponding dictionary file (but I admit it is ugly, inconsistent with other equations and would behave badly if other kinds of fvOptions are added to the same equation, I believe). Another solution could be maybe to replace fvm::SuSp() with fvm::Sp() which should enforce implicit treatment in all cases, if I got it right. I didn't try yet, though. Anyway, I acknowledge that after some reasoning one can recognize and solve this issue, but I believe that one must have at least some experience with the source code to find out, understand and correct it. If fvOptions are intended to provide all users (more basic ones included) with some run-time tool to add simple sources to equations (with semiImplicitSource being the simplest and most common one), I believe that such behaviour can be regarded as a bug. Basically, if you try to use it as explained in the official documentation things won't simply work in many cases with no apparent reason and no possibility of solution without touching the source code. Did you try to file some kind of official issue report about this? Andrea