Information exchange in region coupled solvers
 

Hello Foamers
I have been delving into openFoam (v4.1 extend) codes for a couple of weeks to acquire a better understanding about coupledFvMatrix and the way it uses to transfer information between two regions IMPLICITLY. So, I choosed 'conjugateHeatFoam to study for a test case where chtRcTemperature is applied as boundary condition in interface of fluid and solid. I'd like to know:
1- What kind of equations does chtRcTemperature applies?
2- How does openFoam exchange the information in coupledFvmatrix?
I more interested to have some general information based on the class functionality, e.g. LduInterface, regionCouplingFvPatchField, and chtRcTemperatureFvPatchScalarField

Thanks In Advanced

coupling patch in conjugateHeatFoam
 

Hello Foamers
I asked a question about implementing coupling patch in conjugateHeatFoam two weeks ago, and still, I could see no answer there. I hope if the following question is vague, let me know to explain it more.
in chtRctemperaturescalarField.C there is source() function which is added to matrix[doaminI].source(). I can not understand what's the reason behind that, theoretically! I studied the codes for solution procedure in conjugateheatFoam; I found out matrixUpdateBuffer_, and coupleBoundaryCoeffs() associated to this patch adds another source term to the matrix[domainI].source(), later (after calling coupledldusolver). So, what's the reason for manipulateMatrix() in this boundary condition as below?
source()=kOwn*(Tw - TcOwn) - k*(TcNei - TcOwn);
matrix.source()+=source()*magSf

Thanks in advance

How to perform preconditioning of region coupled matrices?
 

1. The equations that chRcTemperature solves can be found in the boundary condition source files.

2. The openfoam exchange format for the coupledFvMatrix is pretty innovative IMHO. It creates two separate matrices using the rows of the coupledMatrix system. The coupling between the two happens through the interface boundary condition such as the derived boundary conditions from regionCouplingFvPatchField. In this interface boundary condition, with the two meshes coupled, the coupling happens at the time of matrix inversion. Using iterative solvers the information is transferred using the initMatrixInterfaces function and the updateMatrixInterfaces in coupledLduMatrix.C. This allows for implicit coupling during the AMul process.


I do have some concerns about this method as far as preconditioning is concerned. I am yet to find the appropriate section of code in the framework that considers the coupling coefficients for the matrix coupling in the preconditioner. I would appreciate if anybody on this forum can guide me with that?

egards,
Ripu

About matrix coupling, as far as I got from the source code, CouplingPatchField possess functions including valueInteranalCoeffs(), BoundaryInternalCoeffs(), gradientInternalCoeffs() and ngradientBoundaryCoeffs() which are added to source of the matrix in coupledFvMatrix.solve() function. internalCoeffs() are added to diag like other types of fvPatchFields but boundaryCoeffs() of couplingPatchField are multiplied by matrixbufferUpdate and added to source in solution loop

You are too kind to call it innovative, for me its rather an attempt at fixing a very serious design flaw in the software.

There should not be need for so many threads that crop up time to time trying to make multi-region implicit treatment for various physics models. Even this thread would not be needed if solution fixed the problem.

I guess, innovation happens when pushed against the wall :)

I do agree with you that had this "design flaw" not been there, there would be lesser posts on this topic. I hope someone is working on fixing this.

Do you have any opinion about the concern I raised below? - "I do have some concerns about this method as far as preconditioning is concerned. I am yet to find the appropriate section of code in the framework that considers the coupling coefficients for the matrix coupling in the preconditioner. I would appreciate if anybody on this forum can guide me with that?"

Thanks for reply
What do you mean by 'IMHO'? Can you introduce me a document to understand the procedure?

As for as I know, in the boundary condition "source()=kOwn*(Tw - TcOwn) - k*(TcNei - TcOwn);" equals to zero. so do nothing for contribution to the source term.


In file "coupledFvScalarMatrix.C", there is a line, i.e.

curMatrix.addBoundaryDiag(curMatrix.diag(), 0); to add the boundary contribution to diag(), as for the boundary contribution for off diag, I guess it is done through 'updateMatrixInterface'.


By the way, In file "coupledFvScalarMatrix.C", there is a code line, i.e.

curMatrix.addBoundarySource(source[rowI], 0); which really does nothing, because

its second parameter is false, so according to the definition of this function, it does nothing.


Please anyone make comments, to discuss if i am right. Thanks a lot

"source()=kOwn*(Tw - TcOwn) - k*(TcNei - TcOwn);" is not zero. I suggest to check it with either solution of a case or study code to get how KOwn is calculated.


For curMatrix.addBoundaryDiag(curMatrix.diag(), 0) you are right. But note the line:


if (!ptf.coupled())
{

addToInternalField(lduAddr().patchAddr(patchI), pbc, source);
}






means that if a patch is not coupled type, "source" will be added to "boundaryCoeffs_" of that patch.

For conjugate Heat this patch field is coupled, so the line you quoted would not be executed. So as I said it do nothing. I did invested the code, by now, for fluid and solid thermal couple, that source is equal to zero.

sorry I mean it do nothing for coupled patch. and source is zero. So this is the reason why I guess off diag contribution from coupled patch is from update Matrix

For this special case you are right and this lines can be ignored. But maybe coupledFvMatric.C is written for more general coupled equations.
In coupledFvMatrix.C, there are two lines:


#curMatrix.completeAssembly();

# source.set(rowI, new scalarField(curMatrix.source()));


For the solver we are talking about, chtRcTemperatureFvPatchScalarField::Manipulatematr ix() is called form the first line which evaluates source term for chtRcTemperatureFvPatchScalarField (a coupled type) patch. However, In general, we might have other non-coupled patches that have source terms like what Manipulatematrix() does.



Did you have any radiation at interface or other source except conduction heat transfer? Look at this line in chtRcThermalDiffusivityFvPatchScalarField.C:

#k[facei] = max(min(k[facei], 100*kHarm[facei]), 0.01*kHarm[facei]);


Meaning that KOwn is not exactly kHarm all the time.

Hi, I only investigate fluid and solid thermal couple without radiation, which in this case all terms related radiation should be zero. for this specific case without radiation, the only thing confused me is this line

K.originalPatchField()/(1 - p.weights())/mld.magDelta(p.index()) in file chtRcTemperatureFvPatchScalarField; or similar relative line in file chtRcThermalDiffusivityFvPatchScalarField.C.

if we assume .originalPatchField() return field at the cell center next to the coupled face, every thing is fine, which leads the source term equal to zero ( not considering radiation), then the diag contribution from the coupled Boundary is obtained as I said in previous post, and I guess the off diag contribution from the coupled Boundary is from updateMatrix.
However, when assuming .originalPatchField() returns just the field at the coupled face, then in this case, every thing from source() function and from calculate thermal conductivity in very strange, as I investigated. And I have tried the case, it proves that calculated value for thermal conductivity at the coupled face is not affected by setting different values at the coupled face in the case folder, which proves ,I think, that .originalPatchField() returns the field at the cell center next to the coupled face.


Please make any comments . Thanks a lot

By the way, by saying
#k[facei] = max(min(k[facei], 100*kHarm[facei]), 0.01*kHarm[facei]);

I think because kHarm has been multiplied by 100 or 0.01, so kHarm is used to, like, bound the value of k[facei], so in most cases, the value should be k[facei] ?


what do you think about this line?

const scalarField kHarm = weights*fOwn + (1.0 - weights)*fNei;
I am confused about fOwn and fNei, as I as in the previous post. It looks like the field at the coupled face, however, if it is the field at the cell next to coupled face, every thing in the code not confused me.

I checked it before and concluded that:
# K.originalPatchField()/(1 - p.weights())/mld.magDelta(p.index())
means original K for a considering domain (not interpolated K) divided by "distance from cell to face for patch". Regarding ".originalPatchField()", from my experience, it returns cell value of K. However you can find it in regionCouplingFvPatchField.C for a better understanding.


This is my idea about effect of regionCoupling patches in solution process:



internalCoeffs_: added to diag() of neighbor cell in the matrix

boundaryCoeffs_: It will be then multiplied by matrixUpdateBuffer_ of the shadow cell and the result is added to matrix.source()

manipulateMatrix(): It will add a new source to matrix.source()

Regarding fOwn and fNei, these are K for Own and neighbor domains. These are used to evaluate kHar via harmonic weights:


# const scalarField& fOwn = owner.originalPatchField();
# const scalarField& lfNei = neighbour.originalPatchField();


lfNei is then interpolated from neighbor to own patch to give fNei.

Thanks for reply. I agree with you that '.originalPatchField()' returns cell value (based on my investigation)


And I agree with your ideas. To make further, for region coupled face, manipulateMatrix() will add source() multipled by surfaceAera to the Matrix source term.
And source() returns kOwn*(Tw - TcOwn) - k*(TcNei - TcOwn); this equation is

( thermalConductivityAtFluidCell / distanceFromCellNodeToWall ) * (Twall - TFluidCell) -
( thermalConductivityAtRegionCoupleFace / distanceFromFluidCellToSolidCell ) * (TSolidCell -TFluidCell). If we don't consider radiation, based on continuous heat flux, this equation will be zero.


This is my understanding

I agree with you. Except, I am not quite sure for zero value of source term, even without radiation. I have to check it

Yes,

I think in most case, kHar is not the value returned by this function. This member function will return the value calculated by :

k = kOwn*(TwOwn*(kNei*(TcNei - TcOwn) + Qr + fourQro) - TcOwn*fourQro);
k /= stabilise((fourQro + TwOwn*(kOwn + kNei))*(TcNei - TcOwn), SMALL);
k /= p.deltaCoeffs();
this code can be simplified as :

k= distanceFromFluidCellToSolidCell /

( distanceFromFluidCellToWall / ThermalConductivityAtFluidCell + distanceFromSolidCellToWall / ThermalConductivityAtSolidCell ),


And this equation is actually based on continuous Temperature and continuous heat flux at the wall. That is :
distanceFromFluidCellToSolidCell / thermalConductivityAtTheWall =

distanceFromFluidCellToWall / thermalConductivityAtFluidCell +
distanceFromSolidCellToWall / thermalConductivityAtSolidCell
This is actually what returns, we may call it as effective thermal conductivity at the wall.
(note: all I don;t consider radiation. )

And also, if don't consider radiation, kHarm will equal to k, please check.