Hi all,

Please consider me as a beginner in editing the Nektar++ solvers.

I'd like to implement Sub-Grid-Scale Smagorinsky model into IncNavierStokesSolver. This would require:
1. local element size,

2. calculation of the velocity gradient at each point in physical space,

3. adding turbulent viscosity nu_turb=func( size_element, grad(v) ) to the viscosity in physical space.

I had a look into the code focusing on my requirements and have a few questions:

1. Re documentation, is there any developer guide for the IncNavierStokesSolver as a starting point?

2. Re the implementation, I'm using the velocity correction scheme.

Please suggest if the strategy below makes sense (performance-wise):

    2.a) For the solution points separation in space,  repeat the same as in SVVVarDiffCoeff() function to calculate h/p:

                {    h = max(h, exp3D->GetGeom3D()->GetEdge(i)->GetVertex(0)->dist(
                             *(exp3D->GetGeom3D()->GetEdge(i)->GetVertex(1))));
                }
                for(int i = 0; i < 3; ++i)
                {
                    p = max(p,exp3D->GetBasisNumModes(i)-1);
                }

    2.b) Obtain grad(v) using  PhysDeriv(space, var, gradvar). However, I need some help here. Is it possible to access different components of the gradients (du_dx, du_dy, du_dz etc.) of the advected fields without going into the loops as done in ProcessModules/ProcessGrad.cpp?

    2.c) Add the above changes in VelocityCorrectionScheme.cpp as commented:

    /**
     * Solve velocity system
     */
    void   VelocityCorrectionScheme::v_SolveViscous(
        const Array<OneD, const Array<OneD, NekDouble> > &Forcing,
        Array<OneD, Array<OneD, NekDouble> > &outarray,
        const NekDouble aii_Dt)
    {
        StdRegions::ConstFactorMap factors;
        StdRegions::VarCoeffMap varCoeffMap = StdRegions::NullVarCoeffMap;
        MultiRegions::VarFactorsMap varFactorsMap =
            MultiRegions::NullVarFactorsMap;

        AppendSVVFactors(factors,varFactorsMap);

        // // // Smagorinsky model Start // // // 

        // Find ele_size i.e. h/p as in SVVVarDiffCoeff()
        // Find du_dxi, dv_dxi, dw_dxi, dAnotherField_dxi, where i=1,2,3
        // nu_t = sq(const*ele_size)*sqrt(2.0*S_ijS_ij)
        // add nu_t field to m_diffCoeff

        // // // Smagorinsky model End // // //

        // Solve Helmholtz system and put in Physical space
        for(int i = 0; i < m_nConvectiveFields; ++i)
        {
            // Setup coefficients for equation
            factors[StdRegions::eFactorLambda] = 1.0/aii_Dt/m_diffCoeff[i];
            m_fields[i]->HelmSolve(Forcing[i], m_fields[i]->UpdateCoeffs(),
                                   NullFlagList,  factors, varCoeffMap,
                                   varFactorsMap);
            m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(),outarray[i]);
        }
    }

3. Since the velocity field is not C1 continuous at element boundaries, does anyone see a problem due to that?

General comments are welcome on the methodology.

Sorry for a long email.

BR,

Vishal

---
Vishal SAINI

PhD Student

Loughborough University