CFD Online Discussion Forums - inverse of a matrix or solution of implicit equation

CFD Online Discussion Forums (https://www.cfd-online.com/Forums/)

- OpenFOAM Programming & Development (https://www.cfd-online.com/Forums/openfoam-programming-development/)

- - inverse of a matrix or solution of implicit equation (https://www.cfd-online.com/Forums/openfoam-programming-development/100421-inverse-matrix-solution-implicit-equation.html)

inverse of a matrix or solution of implicit equation

Hej,

I am looking for a way to solve the following equation

Code:

a_i = c1 * a_j * du_i / dx_j + c2

for i=1 this looks like

Code:

a_1 = c1 * a_1 * du_1/dx_1 + c1 * a_2 * du_1 / dx_2 + c1 * a_3 * du_1 / dx_3 + c2

where a is a vector, u is a vector and i and j run from 1-3. This can be solved by either matrix inversion where A_ij = (I- (c * du_i / dx_j)) is the coefficient matrix in every cell (I is the identity matrix).

The idea is to rewrite this into an explicit form, by taking a_1 to the other side.

I was also looking into a form to solve this with the build in fvMatrix of OF, unfortunately you can basically only insert div, grad (or other derivative forms) and source. However, source (fvm::Sp(scalar, solvableVariable))cannot take a tensor as it would be required from du_i/dx_j as a multiplication factor for vector a

The other idea is to take the coefficient matrix (A_ij) in every cell and divide the right hand side by this.

Code:

A_ij*a_j = S_i => a_j = A_ij^-1 * S_i

Does anybody know how to do this in OF?

OpenFOAM version used is 2.1.

Further checking this, I found a function inv(const tensorField&), which I believe takes the inverse of a tensor, but I am unsure for what it takes the inverse, either the tensor in every cell or the overall tensor (the field itself).

I now assumed that it takes the inverse of every tensor in the cells and the following code came out for the following equation

(I-A_ij) * a_j = S_i
-> a_j = (I-A_ij)^1 * S_i

Code:

        // velocity gradient matrix, M_ij

      const volTensorField coeffMatrix(-sThetauCoeff*Cthetau2*fvc::grad(U));

        // I_ij - M_ij

      const volTensorField negCoeffMatrix(I-coeffMatrix);

        // (I_ij - M_ij)^1

      const volTensorField invNegCoeffMatrix(inv(negCoeffMatrix));



      thetau = invNegCoeffMatrix & (gSource+gradTSource);

The solution was a little bit easier than first assumed. OpenFOAM has a great operator definition which lets you divide by a tensor :)

So in case you have an equation

A_ij x_j = s_i

and you want to solve for x_j, where A is a second order tensor (for example the gradient of velocity, grad(U) ) and x and s are vectors, the one simply divides by A_ij

as in

Code:

  volVectorField s(g);  // g is the gravity field

  volTensorField A(fvc::grad(U));

  volVectorField x = s / A;

I guess you could do the same with

Code:

 volVectorField s(g); // g is the gravity field

 volTensorField A(fvc::grad(U));

 volTensorField invA(inv(A));  // taking the inverse of A

 volVectorField x = (1/det(A)*invA) & s; // det(A) takes the determinant

the previous solution is of course much simpler and more elegant

Dear Roman,

I need to calculate a vector using anothor vector and a matrix by formula:

Code:

a=[X^T X]^-1 X^T y,

where y is a vector and X is a square matrix (T means transpose and -1 means inverse of a matrix), but the matrix and vectors are not volScalarField.

Is there any possible to do that by openFOAM?

Regards,

Quote:

Originally Posted by atoof (Post 444099)

Dear Roman,

I need to calculate a vector using anothor vector and a matrix by formula:

Code:

a=[X^T X]^-1 X^T y,

X should be a volTensorField I think and then y is a volVectorField, there is an operator for transpose, which is X.T() for a tensor X, so an inverse (^1) you get with inv(X) and a transpose with X.T()