Dear all,

Using a locally compiled version of nektar release 5.6.0 (which runs pretty well) to solve fully 3d plane Poiseuille flow geometry,

I'm trying to add variable or temperature dependant viscosity, by defining a function in the nektar session file, to be evaluated and passed to VarCoeffMap D00, D11,D22 as suggested by M. R. Rouhanian, a while ago, indeed I've come to the conclusion that the Incompressible Navier-Stokes solvers use too much memory in a cumulative way without releasing it throughout the time advancement process.

key points:

- Correct behavior with constant viscosity. no problem, the IncNavierStokesSolver runs good

- Source file modified: nektar-v5.6.0/solvers/IncNavierStokesSolver/EquationSystems/VelocityCorrectionScheme.cpp (Line:890, VelocityCorrectionScheme::v_SolveViscous function)

- Monotonic memory usage increase when enabling variable viscosity through any of the possible 'nu' dependencies i.e.: nu(Temp),nu(x,y,z), nu(t), nu(x,y,z,t).

- The memory usage reaches easily 10 gigabytes after 20 to 30 timesteps, obviously evenly spread over the 4 x IncNavierStokesSolver spawned processes, compared to 4x140MB=560MB for the constant viscosity case

- The memory usage keeps increasing, never decreases, no release of the memory surplus from one timestep to the next, it just keeps adding up...

- No crash the calculations do carry on, with a noticeable slowdown due to viscosity array update, however not usable from a practicality view point as any serious case simulation would claim few terabytes of leaked memory

- Placing some "cout<<" printing arrays sizes show that the arrays involved do keep the same size no issues

It seems to be a garbage collection issue or a missing destructor call that I'm not able to pinpoint or fix.

Any help on this issue or suggestion to fix it are greatly appreciated,

Thank you all in advance, wishing you all the best,

Saad

Code details:

My VelocityCorrectionScheme::v_SolveViscous modified function (file attached herein) reads:

void VelocityCorrectionScheme::v_SolveViscous(
const Array<OneD, const Array<OneD, NekDouble>> &Forcing,
const Array<OneD, const Array<OneD, NekDouble>> &inarray,
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);
ComputeGJPNormalVelocity(inarray, varCoeffMap);

/*
switch (coordim)
{
case 3:
Exp->GetCoords(xc0, xc1, xc2);
break;
default:
ASSERTL0(false, "Coordim not valid");
break;
}*/

////////////////////////////////////////////////////////////////////////////////
/////////// Modifications to account for variable viscosity start here://///////

int nuflag=0; //initialiazed to 0, fall back to constant viscosity if nu(T) is not defined
time_t rawTime;
time(&rawTime);
NekDouble time = static_cast<NekDouble>(rawTime);//ensure time is of NekDouble type

// if the function nuOfT is defined in the session file, evaluate the function and set variable viscosity
if(m_session->DefinesFunction("nuOfT"))
{
Array<OneD, NekDouble> xxc0,xxc1,xxc2;// declare NekDouble arrays to hold coordinates

int nqv;// total points
int coordim=3;// default dimension

LibUtilities::EquationSharedPtr nuOfTfunc;
nqv=m_fields[0]->GetTotPoints();

// allocate 0.0 valued arrays of totalpoints size, to receive the space coordinates
xxc0 = Array<OneD, NekDouble>(nqv, 0.0);
xxc1 = Array<OneD, NekDouble>(nqv, 0.0);
xxc2 = Array<OneD, NekDouble>(nqv, 0.0);

m_fields[0]->GetCoords(xxc0,xxc1,xxc2);
coordim=m_fields[0]->GetCoordim(0);

Array<OneD,NekDouble> nuOfT(nqv,0.0);// 0.0 valued nekdouble array to hold viscosity value at each point of space
nuOfTfunc=m_session->GetFunction("nuOfT",0);

nuOfTfunc->Evaluate(xxc0,xxc1,xxc2,nuOfT);// evaluate the function nu(T)

varCoeffMap[StdRegions::eVarCoeffD00]=nuOfT;// assign the values of viscosity for the three diagonal components (same) evaluated
varCoeffMap[StdRegions::eVarCoeffD11]=nuOfT;// at every point of space using the decalred nu(T) function
varCoeffMap[StdRegions::eVarCoeffD22]=nuOfT;

nuflag=1;// set the variable viscosity flag

}

// Solve Helmholtz system and put in Physical space
for (int i = 0; i < m_nConvectiveFields; ++i)
{

// Add diffusion coefficient to GJP matrix operator (Implicit part)
if (m_useGJPStabilisation)
{
factors[StdRegions::eFactorGJP] = m_GJPJumpScale / m_diffCoeff[i];
}

// Setup coefficients for equation
factors[StdRegions::eFactorLambda] = 1.0 / aii_Dt / m_diffCoeff[i];

if (nuflag==1) factors[StdRegions::eFactorTime]=time;//if the nuOfT nu(T) function is defined set the tume factor,
//else continue with constant viscosity

m_fields[i]->HelmSolve(Forcing[i], m_fields[i]->UpdateCoeffs(), factors,
varCoeffMap, varFactorsMap);
m_fields[i]->BwdTrans(m_fields[i]->GetCoeffs(), outarray[i]);
}

}

Sample function node from the xml session file: