# inverse of a matrix or solution of implicit equation

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

 April 25, 2012, 08:51 inverse of a matrix or solution of implicit equation #1 Senior Member     Roman Thiele Join Date: Aug 2009 Location: Eindhoven, NL Posts: 374 Rep Power: 20 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. __________________ ~roman

 April 25, 2012, 10:35 #2 Senior Member     Roman Thiele Join Date: Aug 2009 Location: Eindhoven, NL Posts: 374 Rep Power: 20 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);``` __________________ ~roman

 April 27, 2012, 05:21 solution found #3 Senior Member     Roman Thiele Join Date: Aug 2009 Location: Eindhoven, NL Posts: 374 Rep Power: 20 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 Tushar@cfd and huangxianbei like this. __________________ ~roman Last edited by romant; April 27, 2012 at 05:32. Reason: addition

 August 6, 2013, 08:02 #4 Member   Hossein Join Date: Apr 2010 Posts: 65 Rep Power: 16 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,

August 6, 2013, 09:02
#5
Senior Member

Roman Thiele
Join Date: Aug 2009
Location: Eindhoven, NL
Posts: 374
Rep Power: 20
Quote:
 Originally Posted by atoof 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()
__________________
~roman