[mm.abdollahzadeh]

Dear all

I was trying to implement the hybrid scheme but I have faced some problems.
I will appreciate any suggestions,

for the implementation I have added three more parameter to the limiter function to be able to define the Peclect number. then I pick the Quick scheme as the template for new scheme.

now the code is reading as :

    scalar limiter
    (
        const scalar cdWeight,
        const scalar faceFlux,
        const typename LimiterFunc::phiType& phiP,
        const typename LimiterFunc::phiType& phiN,
        const typename LimiterFunc::gradPhiType& gradcP,
        const typename LimiterFunc::gradPhiType& gradcN,
        const vector& d,
        const vector& Sf,
        const scalar& magSf ,
       const scalar Diffu
    ) const
    {

    scalar Pe = (faceFlux/magSf)*(d & Sf)/(Diffu*magSf);
    scalar QLimiter;
    QLimiter = pos(1-0.5*mag(Pe));

        // Limit the limiter between upwind and downwind
        return  QLimiter;
    }

the code compiles well. but then when I want to run a case it gives me floating point error. if I put the Qlimiter 0 or 1 it runs without problem. but as I use a function of Pe it wont run.

best
Mahdi

[wouter]

Dear Mahdi,
sorry if stating the obvious: prevent magSf and or Diffu from being zero or add SMALL to the divider

Hope anyone else has a better solution
Best regards
Wouter

[mm.abdollahzadeh]

Dear Wouter

your suggestion will make it work. But then I'm just worried what will happen to boundaries ( possible place that Diffu and Magsf could be zero?!)

best
Mahdi

[mm.abdollahzadeh]

now it runs

Im reading the Diffu in limitedscheme class in this way

const volScalarField& Diff =
phi.db().objectRegistry::lookupObject<volScalarFie ld>(Diffname_);
tmp<surfaceScalarField> tDiffu=fvc::interpolate(Diff);

const surfaceScalarField& Diffu =tDiffu();

best
mahdi