[stevemoore1981]

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?

[santiagomarquezd]

You are not considering the BC's contributions, take a look of this:

https://openfoamwiki.net/index.php/Contrib_gdbOF

Regards

[piccinini]

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.

        fvScalarMatrix pEqn
                (
                    fvm::laplacian(-Mf,p) + fvc::div(phiG)
                );

        pEqn.solve();

[lengweee]

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.

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;

}

%% 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)));

[kkamau]

I am pretty sure you are interchanging owners with neighbours in the above code. May be that is the source of the error

[Jesper_Roland]

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]

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.

[Jesper_Roland]

I do not know how it works for parallel computations.

Best,

Jesper