Hi all,

   I'm trying to implement a stabilized method for  incompressible  N-S equation. Basically, the form of the stabilized term  looks like nu*L_u, where nu is the the artificial viscosity which is varying across quadratue points, L_u is the Laplacian operation of a velocity component.  My plan is using the explicit scheme to discretize this stabilized term to avoid assembling the system matrix at every time step.  I'd like to add the stabilized term to the right hand side in the viscous correction stage of the velocity correction scheme. 

   However I'm not very clear how to implement it using Nektar++  for the case of 3DHomogeneous1D . Here I attached a piece of my current code,  hope someone could help me to correct the mistakes.

 

// velocity_array[0] stored the first component of velocity vector in Fourier space

// z direction is homogeneous

// nu stored artificial viscosity

// np is number of quadrature points in a 2D plane

// nc is number of SEM coefficients in a 2D plane

 

Vmath::Vadd(np*num_planes, nu ,1, velocity_array[0],1,velocity_array[0],1);

 

m_fields[1]->PhysDeriv(0,velocity_array[0],velocity_derv[0]); // gradient in x direction

m_fields[1]->PhysDeriv(1,velocity_array[0],velocity_derv[1]); // gradient in y direction

            

cnt  = 0;

cnt1 = 0;

for(unsigned int i=0; i<num_planes; ++i)

  {

double beta  = 2*M_PI*m_trans->GetK(i)/m_homoLen;

       beta *= beta;

 

    m_fields[0]->GetPlane(i)->IProductWRTDerivBase(0, &velocity_derv[0]+cnt, &wkCoef1[0] +cnt1);                                                         m_fields[1]->GetPlane(i)->IProductWRTDerivBase(1, &velocity_derv[1]+cnt, &wkCoef1[0]+cnt1);

 

    m_fields[1]->GetPlane(i)->IProductWRTBase(tmp1 = &velocity_array[0]+i*np, &wkCoef2[0]+i*nc,MultiRegions::eGlobal);

    Vmath::Smul(nc, beta, &wkCoef2[0] + i*nc, 1, &wkCoef2[0]+ i*nc, 1);

 

cnt  += np;

cnt1 += nc;

   }

  Vmath::Vadd(NumTotCoeffs,wkCoef0,1,wkCoef1,1,visc,1);

  Vmath::Vadd(NumTotCoeffs,wkCoef0,1,wkCoef2,1,visc,1);