CFD Online Discussion Forums

CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Programming & Development (http://www.cfd-online.com/Forums/openfoam-programming-development/)
-   -   inverse of a matrix or solution of implicit equation (http://www.cfd-online.com/Forums/openfoam-programming-development/100421-inverse-matrix-solution-implicit-equation.html)

romant April 25, 2012 08:51

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.

romant April 25, 2012 10:35

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


romant April 27, 2012 05:21

solution found
 
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

atoof August 6, 2013 08:02

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,

romant August 6, 2013 09:02

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,
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,

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()


All times are GMT -4. The time now is 08:56.