
[Sponsors] 
fvMatrix, fvOptions and SuSp: automatic implicit/explicit sourceterm treatment 

LinkBack  Thread Tools  Search this Thread  Display Modes 
December 31, 2016, 11:38 
fvMatrix, fvOptions and SuSp: automatic implicit/explicit sourceterm treatment

#1 
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 
fvMatrix, fvOptions and SuSp: automatic implicit/explicit sourceterm 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 ); … 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 } Code:
70 source.addSup(mtx, fieldI); So, initial equation Code:
fvm::ddt(T) + fvm::div(phi, T)  fvm::laplacian(DT, T) == fvOptions(T) Code:
fvMatrix_A == fvMatrix_B 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 } 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 } 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 } Code:
1065 lduMatrix::operator=(fvmv); Code:
//"lduMatrixOperations.C" ... 223 void Foam::lduMatrix::operator=(const lduMatrix& A) 224 { 225 if (A.diagPtr_) 226 { 227 diag() = A.diag(); 228 } ... Code:
1066 source_ = fvmv.source_; Code:
fvMatrix_A == fvMatrix_B fvMatrix_C = fvMatrix_A  fvMatrix_B where (in pseudocode): Code:
C.diag() = A.diag()  B.diag(); C.source() = A.source()  B.source(); 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 } 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] } } } } Code:
eqn += Su + fvm::SuSp(Sp, psi); Code:
eqn += fvm::SuSp(Sp, psi); 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_; ... Code:
eqn.diag() += (fvm::SuSp(Sp, psi)).diag(); eqn.source() += (fvm::SuSp(Sp, psi)).source() 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 } Code:
210 fvm.diag() += mesh.V()*max(susp.field(), scalar(0)); Code:
212 fvm.source() = mesh.V()*min(susp.field(), scalar(0)) 213 *vf.internalField(); 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 dictionaryfile 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); } Code:
tmC().diag() = tmA.diag()  tmOL.diag(); Code:
tmC().diag() = tmA.diag()  mesh.V()*Sp.field(); 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); Code:
fvMatrix_C.diag() = fvMatrix_A.diag()  mesh.V()*Sp.field; 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; Code:
fvMatrix_C.source() = fvMatrix_A.source() + (mesh.V()*(Sp.field())*T.internalField()); Code:
fvMatrix_C.source() = fvMatrix_A.source() + mesh.V()*Sp.field()*T.internalField(); Code:
fvMatrix_A + SuSp(Sp, T) == 0; 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); Thanks. Last edited by Zeppo; January 2, 2017 at 07:21. 

February 23, 2017, 09:05 

#2 
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); Code:
2) fvMatrix_A  SuSp(Sp, T); 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); Most matrix equations in OpenFOAM's solvers look like: Code:
fvMatrix_A == fvOptions(T); Code:
fvMatrix_A == SuSp(Sp, T); Code:
fvMatrix_A == SuSp(Sp, T); 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 } 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 } 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. 

April 9, 2019, 22:44 

#3 
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 13, 2019, 13:55 

#4  
Senior Member
Sergei
Join Date: Dec 2009
Posts: 261
Rep Power: 21 
Solution is given in my second post:
Quote:
Quote:
Quote:


April 23, 2020, 05:15 

#5 
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 hardcoded a semiimplicit 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 runtime 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 

April 23, 2020, 14:46 

#6 
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 } 

May 5, 2020, 13:30 

#7 
Senior Member
René Thibault
Join Date: Dec 2019
Location: Canada
Posts: 114
Rep Power: 6 
Hi everyone!
I tried the fvOptions with scalarCodedSource to link a source term on epsilon to the ke 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.3474e05; 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 #{ #}; } // ************************************************************************* // Last edited by Tibo99; May 5, 2020 at 16:07. 

December 15, 2021, 11:20 

#8 
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 

Tags 
fvmatrix, fvoptions, implicit/explicit, matrix diagonal dominance, susp 
Thread Tools  Search this Thread 
Display Modes  

