[mrishi]

I am working to change the default divDevRhoReff() defined in turbulenceModels/linearViscousStress.C. The change requires replacing the laplacian(U) operation with div(grad(U)). I change the content of divDevRhoReff(U) as follows:

 
template<class BasicTurbulenceModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff
(
    volVectorField& U
) const
{
/*
    //first change
    Info<<"broken into two: nu + nut" <<endl;
    return
    (
      - fvc::div((this->alpha_*this->rho_*this->nu())*dev2(T(fvc::grad(U))))    //split
      - fvc::div((this->alpha_*this->rho_*this->nut())*dev2(T(fvc::grad(U))))   //split
      - fvm::laplacian(this->alpha_*this->rho_*this->nu(), U)   //not changed
      - fvm::laplacian(this->alpha_*this->rho_*this->nut(), U)
     );
//no issues as expected; some overhead due to repeated grad calculations.

 */
    //second change
    Info<<"changed fvm::laplacian to fvm::div(fvc::grad) for nut term" <<endl;
    const volTensorField A = fvc::grad(U); //to create type required by fvm::div
    const surfaceScalarField N = fvc::interpolate(this->nut()); //surfacescalarField required by fvm::div, created out of volScalarField nut
    return
    (
//      - fvc::div((this->alpha_*this->rho_*this->nu())*dev2(T(fvc::grad(U))))  //split
//      - fvc::div((this->alpha_*this->rho_*this->nut())*dev2(T(fvc::grad(U)))) //split
//        - fvm::laplacian(this->alpha_*this->rho_*this->nu(), U)       //not changed

        - fvm::div(N, fvc::grad(U)())                           // changed, removed alpha,rho to remove clutter.
//      - fvm::div(this->alpha_*this->rho_*(fvc::interpolate(this->nut())),A) //changed to div(grad(U)) instead of laplacian(U) , same as above
    );
}

In the code marked second change, I have commented out the non-problematic components so as to isolate the issue. It generates following error:

../turbulenceModels/lnInclude/linearViscousStress.C:120:5: error: could not convert ‘Foam::operator-(const Foam::tmp<Foam::fvMatrix<Type> >&) [with Type = Foam::Tensor<double>]()’ from ‘Foam::tmp<Foam::fvMatrix<Foam::Tensor<double> > >’ to ‘Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > >’

The issue is type incompatibility of the fvMatrix thus created, which is of type Tensor instead of Vector. While the div operator is rank-reducing, I suppose doing this in an implicit form using fvm may create such an issue. How can I circumvent this situation while still creating a proper fvVectorMatrix which adds correctly to my UEqn.H?
PS: The reason I need to do this is because I need to carry out some operations on the grad(U) before calculating its final divergence. In other words, my nu cannot be taken out of the second derivative directly, so I cannot use laplacian(nut,U).

[simrego]

Hi!


First of all I think you should read about fvc and fvm. fvc is an explicit operator, while fmv is an implicit and based on your code I guess that you are not familiar with these.
Instead of fvm::div(fvc::grad(<your stuff>)) you should use: fvc::div(fvc::grad(<your stuff>)) and treat this term explicitly.

[Santiago]

Originally Posted by mrishi

I am working to change the default divDevRhoReff() defined in turbulenceModels/linearViscousStress.C. The change requires replacing the laplacian(U) operation with div(grad(U)). I change the content of divDevRhoReff(U) as follows:

 
template<class BasicTurbulenceModel>
Foam::tmp<Foam::fvVectorMatrix>
Foam::linearViscousStress<BasicTurbulenceModel>::divDevRhoReff
(
    volVectorField& U
) const
{
/*
    //first change
    Info<<"broken into two: nu + nut" <<endl;
    return
    (
      - fvc::div((this->alpha_*this->rho_*this->nu())*dev2(T(fvc::grad(U))))    //split
      - fvc::div((this->alpha_*this->rho_*this->nut())*dev2(T(fvc::grad(U))))   //split
      - fvm::laplacian(this->alpha_*this->rho_*this->nu(), U)   //not changed
      - fvm::laplacian(this->alpha_*this->rho_*this->nut(), U)
     );
//no issues as expected; some overhead due to repeated grad calculations.

 */
    //second change
    Info<<"changed fvm::laplacian to fvm::div(fvc::grad) for nut term" <<endl;
    const volTensorField A = fvc::grad(U); //to create type required by fvm::div
    const surfaceScalarField N = fvc::interpolate(this->nut()); //surfacescalarField required by fvm::div, created out of volScalarField nut
    return
    (
//      - fvc::div((this->alpha_*this->rho_*this->nu())*dev2(T(fvc::grad(U))))  //split
//      - fvc::div((this->alpha_*this->rho_*this->nut())*dev2(T(fvc::grad(U)))) //split
//        - fvm::laplacian(this->alpha_*this->rho_*this->nu(), U)       //not changed

        - fvm::div(N, fvc::grad(U)())                           // changed, removed alpha,rho to remove clutter.
//      - fvm::div(this->alpha_*this->rho_*(fvc::interpolate(this->nut())),A) //changed to div(grad(U)) instead of laplacian(U) , same as above
    );
}

In the code marked second change, I have commented out the non-problematic components so as to isolate the issue. It generates following error:

../turbulenceModels/lnInclude/linearViscousStress.C:120:5: error: could not convert ‘Foam::operator-(const Foam::tmp<Foam::fvMatrix<Type> >&) [with Type = Foam::Tensor<double>]()’ from ‘Foam::tmp<Foam::fvMatrix<Foam::Tensor<double> > >’ to ‘Foam::tmp<Foam::fvMatrix<Foam::Vector<double> > >’

The issue is type incompatibility of the fvMatrix thus created, which is of type Tensor instead of Vector. While the div operator is rank-reducing, I suppose doing this in an implicit form using fvm may create such an issue. How can I circumvent this situation while still creating a proper fvVectorMatrix which adds correctly to my UEqn.H?
PS: The reason I need to do this is because I need to carry out some operations on the grad(U) before calculating its final divergence. In other words, my nu cannot be taken out of the second derivative directly, so I cannot use laplacian(nut,U).

I dont understand the reason why you would replace the laplacian with div-grad operator, assuming it was possible. It will be a slower and less precise operation.

[simrego]

Originally Posted by Santiago

I dont understand the reason why you would replace the laplacian with div-grad operator, assuming it was possible. It will be a slower and less precise operation.

I guess this is not his "final" code since he wrote:

"The reason I need to do this is because I need to carry out some operations on the grad(U) before calculating its final divergence."

[Santiago]

Originally Posted by simrego

I guess this is not his "final" code since he wrote:

"The reason I need to do this is because I need to carry out some operations on the grad(U) before calculating its final divergence."

then first what you need to manipulate is the snGrad class, and then use surfaceIntegrate to obtain the divergenge of said snGrad