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