CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM

Extracting the matrix data from an fvVectorMatrix

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

Like Tree2Likes
  • 1 Post By Jesper_Roland
  • 1 Post By Jesper_Roland

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   February 22, 2016, 06:30
Default Extracting the matrix data from an fvVectorMatrix
  #1
New Member
 
Steve Moore
Join Date: Sep 2012
Posts: 3
Rep Power: 10
stevemoore1981 is on a distinguished road
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 labelUList& neighbor = M.lduAddr().lowerAddr();
const labelUList& owner = M.lduAddr().upperAddr();
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
File Type: pdf Comparison.pdf (184.7 KB, 97 views)
stevemoore1981 is offline   Reply With Quote

Old   April 16, 2016, 22:49
Default
  #2
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 452
Rep Power: 20
santiagomarquezd will become famous soon enough
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
santiagomarquezd is offline   Reply With Quote

Old   November 12, 2017, 15:09
Default
  #3
New Member
 
Join Date: Mar 2009
Location: Sao Jose dos Campos, Brazil
Posts: 29
Rep Power: 13
piccinini is on a distinguished road
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();
piccinini is offline   Reply With Quote

Old   August 25, 2018, 04:51
Default
  #4
New Member
 
wei leng
Join Date: Mar 2014
Posts: 3
Rep Power: 8
lengweee is on a distinguished road
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<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();











    //
    //
    // Output Boundary condition
    //
    //
    for (int k = 0; k < 3; k++) {
	char name[1000];
	char surfix[3][100] = {
	    "x", "y", "z"
	};
	sprintf(name, "Bdry%d_%s.m", timeIndex, surfix[k]);
	FILE *fp = fopen(name, "w");
	fprintf(fp, "bc%s = [\n", surfix[k]);
	    
	forAll(M.internalCoeffs(), patchI)
	    {
		const label*  fPtr = M.lduAddr().patchAddr(patchI).begin();
		Field<double> intF(M.internalCoeffs()[patchI].component(k));
		Field<double> 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)));
lengweee is offline   Reply With Quote

Old   March 22, 2019, 21:02
Default fvMesh
  #5
New Member
 
-- Country other than USA or Canada --
Join Date: Mar 2019
Posts: 3
Rep Power: 3
kkamau is on a distinguished road
I am pretty sure you are interchanging owners with neighbours in the above code. May be that is the source of the error
kkamau is offline   Reply With Quote

Old   October 14, 2020, 12:15
Default
  #6
New Member
 
Jesper R. K. Qwist
Join Date: Dec 2017
Posts: 18
Rep Power: 4
Jesper_Roland is on a distinguished road
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<scalar> wdiag
(
IOobject
(
"diag",
mesh.time().timeName(),
mesh,
IOobject::NO_READ
),
diag
);

const scalarField& upper = p_rghEqn.upper();
IOField<scalar> wupper
(
IOobject
(
"upper",
mesh.time().timeName(),
mesh,
IOobject::NO_READ
),
upper
);

const scalarField& source = p_rghEqn.source();
IOField<scalar> 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 06:16. Reason: Added alternative solution
Jesper_Roland is offline   Reply With Quote

Old   October 19, 2020, 14:12
Default
  #7
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 110
Rep Power: 8
dlahaye is on a distinguished road
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.
dlahaye is offline   Reply With Quote

Old   October 20, 2020, 08:50
Default
  #8
New Member
 
Jesper R. K. Qwist
Join Date: Dec 2017
Posts: 18
Rep Power: 4
Jesper_Roland is on a distinguished road
I do not know how it works for parallel computations.

Best,

Jesper
dlahaye likes this.
Jesper_Roland is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
[General] Extracting ParaView Data into Python Arrays Jeffzda ParaView 19 June 11, 2020 16:54
compiling fails: UDF for extracting mixture fraction data Nekrokeks Fluent UDF and Scheme Programming 5 June 22, 2015 08:37
Extracting Raw Data from STARCCM+ rks171 STAR-CCM+ 3 December 30, 2011 18:10
Problems with extracting data Alicelin OpenFOAM Post-Processing 0 January 25, 2010 09:54
Extracting Stream Function from nodal u,v data Abhijit Tilak Main CFD Forum 5 March 8, 2004 17:13


All times are GMT -4. The time now is 06:26.