# 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: Stockholm, Sweden Posts: 359 Rep Power: 11 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: Stockholm, Sweden Posts: 359 Rep Power: 11 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: Stockholm, Sweden Posts: 359 Rep Power: 11 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: 62 Rep Power: 7 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: Stockholm, Sweden
Posts: 359
Rep Power: 11
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

 Thread Tools Display Modes Linear Mode

 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 OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post cuteapathy CFX 14 March 20, 2012 07:45 colopolo CFX 13 October 4, 2011 22:03 mcaro Main CFD Forum 0 January 25, 2011 07:52 Phiper Main CFD Forum 4 December 10, 2010 02:58 bzz77 Main CFD Forum 0 March 25, 2010 17:31

All times are GMT -4. The time now is 04:37.