CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Programming & Development (
-   -   implicit div(k*grad(a*u)) (

gunnar2 August 8, 2011 14:45

implicit div(k*grad(a*u))
I have asked in "Running / Solving / CFD (" before, but now this seems to be a programming question.

How can I treat \nabla\cdot k \nabla \alpha u implicitly?

k - scalar field (known diffusivity)
\alpha - scalar field (known mass fraction)
u - scalar field (unknown)

fvm::laplacian(k,alpha*u) generates a compiler error.

I think a possible solution is to right-multiply the fvMatrix fvm::laplacian(k,u) with a diagonal matrix D_\alpha made up of the values of \alpha on it's diagonal.

But how do I write this in C++?

Thank you

gunnar2 August 9, 2011 15:22

Is it really not possible to build a diagonal matrix?

Another possibility can be first to construct

fvMatrix M=fvm::laplacian(k,u)

In the next step I would have to multiply each row i of the matrix M with alpha[i]. However I don't understand the storage format of fvMatrix. May I ask if somebody would share a few lines of code, how to multiply one row of a fvMatrix with a scalar value?

Thank you

gunnar2 August 11, 2011 17:40

anybody any idea?

or at least a reason, why there cannot be an implicit scheme for \nabla\cdot k\nabla\alpha u?

Thank you

marupio August 12, 2011 09:52

Can you use the chain rule to seperate alpha and u?

jakobh August 12, 2011 12:41

Hi Gunnar!
You can use this function:


namespace Foam{ namespace fvm{
  template<class Type, class GType>
  tmp<fvMatrix<Type> >
    const GeometricField<GType, fvPatchField, volMesh>& gamma,
    GeometricField<Type, fvPatchField, volMesh>& vf,
    GeometricField<Type, fvPatchField, volMesh>& alpha
    tmp<fvMatrix<Type> > Laplacian=fvm::laplacian(gamma,vf);
    return Laplacian;
} }

In your case: fvm::laplacian(k,u,alpha)


ziemowitzima August 21, 2011 15:57

Hi Jakob,
I read your post and it seems that you can know the answer for my question with laplacian operator.
If you dont mind please read my post at:

All times are GMT -4. The time now is 01:41.