# Extracting the matrix data from an fvVectorMatrix

 Register Blogs Members List Search Today's Posts Mark Forums Read

February 22, 2016, 05:30
Extracting the matrix data from an fvVectorMatrix
#1
New Member

Steve Moore
Join Date: Sep 2012
Posts: 3
Rep Power: 11
Hi

I'm currently trying to modify the pimpleFoam solver to incorporate some optimal control theory into it. In doing so, I'm trying to obtain a matrix that describes the converged solution at each time step... meaning something like A U = b, where U is the converged discrete velocity field at each time step. I was assuming that if the solution has converged at a given time step, then one could create an fvVectorMatrix:

fvVectorMatrix UEqnTest(fvm::ddt(U) + fvm::div(phi,U) == fvm::laplacian(nu,U) - fvc::grad(p));

and then extract the diag, upper, and lower data, and the source term. I was passing this object into a function (along with the converged velocity field):

write("System", runTime.timeIndex(), UEqnTest, U);

at each time step and writing out the matrix data (and source term) into an .m file as a list of {row, col, val} triplets that I could use to create a sparsematrix object in Matlab, solve the linear system there (i.e. U=A\b), and compare to the velocity field written out by OpenFOAM. The write function looks like:

void write(word matrixName, label timeIndex, fvVectorMatrix M, volVectorField& U)
{
std::stringstream fileName;
char timeStep[8];

sprintf(timeStep, "%05d", timeIndex);

fileName << matrixName << timeStep << ".m";

std::fstream file(fileName.str().c_str(), std::ios:ut);

const scalarField& diag = M.diag();
const vectorField& source = M.source();
const scalarField& lower = M.lower();
const scalarField& upper = M.upper();

file << "Data = [" << std::endl;
for(label i=0; i<diag.size(); i++)
{
file << i+1 << "\t" << i+1 << "\t" << diag[i] << std::endl;
}
for(label f=0; f<upper.size(); f++)
{
file << owner[f]+1 << "\t" << neighbor[f]+1 << "\t" << lower[f] << std::endl;
file << neighbor[f]+1 << "\t" << owner[f]+1 << "\t" << upper[f] << std::endl;
}
file << "];\n " << matrixName << "Matrix = sparse(Data(:,1), Data(:,2), Data(:,3));\n";
file << "Source = [" << std::endl;
for(label i=0; i<source.size(); i++)
{
file << i+1 << "\t" << source[i][0] << "\t" << source[i][1] << "\t" << source[i][2] << std::endl;
}
file << "];\n";
file << "Solution = [" << std::endl;
for(label i=0; i<U.size(); i++)
{
file << i+1 << "\t" << U[i][0] << "\t" << U[i][1] << "\t" << U[i][2] << std::endl;
}
file << "];\n";
file.close();

return;

}

What I've found is that the OpenFOAM solution and the Matlab solution (computing U=A\b for the A and b that I wrote out from this fvVectorMatrix) are similar, but there is a significant difference. I've attached a .pdf plotting the two solutions for the Ux component of the solution for a simple 2D channel flow (generated in blockMesh as a simple rectangular domain, 1 block with 100 x 15 x 1 cells).

My questions essentially boils down to whether I am correct in assembling a system of equations this way (i.e. the fvVectorMatrix object that I construct inside the time marching loop, but after the pimple.loop() while loop) and assuming that solving that system should give me the converged velocity field for that time step?
Attached Files
 Comparison.pdf (184.7 KB, 125 views)

 April 16, 2016, 21:49 #2 Senior Member     Santiago Marquez Damian Join Date: Aug 2009 Location: Santa Fe, Santa Fe, Argentina Posts: 452 Rep Power: 21 You are not considering the BC's contributions, take a look of this: https://openfoamwiki.net/index.php/Contrib_gdbOF Regards __________________ Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC) - CONICET/UNL Tel: 54-342-4511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe - Argentina. http://www.cimec.org.ar

 November 12, 2017, 14:09 #3 New Member   Join Date: Mar 2009 Location: Sao Jose dos Campos, Brazil Posts: 29 Rep Power: 15 Hello, Will that work (with bc contributions) If I export the matrix after solve() is called? That is, after pEqn.solve() is called in the code bellow? Thanks is advance. Code: fvScalarMatrix pEqn ( fvm::laplacian(-Mf,p) + fvc::div(phiG) ); pEqn.solve();

 August 25, 2018, 03:51 #4 New Member   wei leng Join Date: Mar 2014 Posts: 3 Rep Power: 10 Thanks a lot for the piece of code to export foam internal matrix to matlab. I add the boundary matrix output, based on your code, please find in the code below. Also attached is a simple matlab script to vertify that A * x = b. The function is based on vector, scalar case is similar. Code: void write(word matrixName, label timeIndex, fvVectorMatrix M, volVectorField& U) { std::stringstream fileName; char timeStep[8]; sprintf(timeStep, "%05d", timeIndex); fileName << matrixName << timeStep << ".m"; std::fstream file(fileName.str().c_str(), std::ios::out); file.precision(15); const scalarField& diag = M.diag(); const vectorField& source = M.source(); const unallocLabelList& neighbor = M.lduAddr().lowerAddr(); const unallocLabelList& owner = M.lduAddr().upperAddr(); const scalarField& lower = M.lower(); const scalarField& upper = M.upper(); file << "Data = [" << std::endl; for(label i=0; i intF(M.internalCoeffs()[patchI].component(k)); Field bouF(M.boundaryCoeffs()[patchI].component(k)); forAll(M.lduAddr().patchAddr(patchI), faceI) { label fCell = fPtr[faceI]; //diag(fCell) -= intF[faceI]; fprintf(fp, "%5d %30.15e %30.15e\n", fCell+1, intF[faceI], bouF[faceI] ); } } fprintf(fp, "];\n"); fclose(fp); } return; } Code: %% Read solution UEqn00001; Bdry1_x; Bdry1_y; Bdry1_z; n=length(Source); I = Source(:,1); bx = ones(n,1); by = bx; bz = bx; bx(I) = Source(:,2); by(I) = Source(:,3); bz(I) = Source(:,4); I = Solution(:,1); ux = ones(n,1); uy = ux; uz = ux; ux(I) = Solution(:,2); uy(I) = Solution(:,3); uz(I) = Solution(:,4); %% Modify with BC A0 = UEqnMatrix; % x A = A0; fx = bx; I=bcx(:,1); for i = 1:length(I) A(I(i),I(i)) = A(I(i),I(i)) + diag(bcx(i,2)); fx(I(i)) = fx(I(i)) + bcx(i,3); end fprintf('resx: %e\n', max(abs(A * ux - fx)));

 March 22, 2019, 20:02 fvMesh #5 New Member   -- Country other than USA or Canada -- Join Date: Mar 2019 Posts: 3 Rep Power: 5 I am pretty sure you are interchanging owners with neighbours in the above code. May be that is the source of the error

 October 14, 2020, 11:15 #6 New Member   Jesper R. K. Qwist Join Date: Dec 2017 Posts: 20 Rep Power: 6 owners and neighbours have been interchanged in code by lengweee. Reference: fvMesh.H - Line 375 - 387 //- Internal face owner. Note bypassing virtual mechanism so // e.g. relaxation always gets done using original addressing const labelUList& owner() const { return fvMesh::lduAddr().lowerAddr(); } //- Internal face neighbour const labelUList& neighbour() const { return fvMesh::lduAddr().upperAddr(); } Alternative solution to extract matrix. Put this code in solver code to extract pressure equation matrix system. You can modify it to work for a vector matrix as well. const scalarField& diag = p_rghEqn.diag(); IOField wdiag ( IOobject ( "diag", mesh.time().timeName(), mesh, IOobject::NO_READ ), diag ); const scalarField& upper = p_rghEqn.upper(); IOField wupper ( IOobject ( "upper", mesh.time().timeName(), mesh, IOobject::NO_READ ), upper ); const scalarField& source = p_rghEqn.source(); IOField wsource ( IOobject ( "source", mesh.time().timeName(), mesh, IOobject::NO_READ ), source ); // Assigning contribution from BC forAll(p_rgh.boundaryField(), patchI) { const fvPatch &vfp = p_rgh.boundaryField()[patchI].patch(); forAll(vfp, faceI) { label cellI = vfp.faceCells()[faceI]; wdiag[cellI] += p_rghEqn.internalCoeffs()[patchI][faceI]; wsource[cellI] += p_rghEqn.boundaryCoeffs()[patchI][faceI]; } } // write to file wdiag.write(); wsource.write(); wupper.write(); Now you can load the data from the text files into you favourite data processing tool and assemble the matrix and solve it. dlahaye likes this. Last edited by Jesper_Roland; October 19, 2020 at 05:16. Reason: Added alternative solution

 October 19, 2020, 13:12 #7 Senior Member   Domenico Lahaye Join Date: Dec 2013 Posts: 431 Blog Entries: 1 Rep Power: 13 Wonderful, thanks! Suppose now targeting the parallel case, i.e., suppose the linear system A u = b to be decomposed into 2 non-overlapping domains such that the vectors can be decomposed as u = [u_0; u_1], b = [b_0; b_1] and the matrix to be decomposed as A = [A_{00} A_{01} ; A_{10} A_{11}]. Can u_i, b_i and A_{ij} be printed seperately? Thanks again. Domenico Lahaye.

 October 20, 2020, 07:50 #8 New Member   Jesper R. K. Qwist Join Date: Dec 2017 Posts: 20 Rep Power: 6 I do not know how it works for parallel computations. Best, Jesper dlahaye likes this.