weak form of diffusion term for 1d homogeneous problem
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);
HI Zhicheng, Since Yan is very familiar with this bit of code I suggest that he liase with you on this issue. Cheers Spencer. On 22 May 2015, at 09:10, Zhicheng Wang <wangzhicheng09@gmail.com<mailto:wangzhicheng09@gmail.com>> wrote: 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); _______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052
H Zhicheng, I can’t get much useful information from your attached code, so I don’t know what’s your problem or mistakes happened in your implementation. Also, I don’t really understand why do you use explicit scheme for your stabilized term, As far as I know, the viscous term in Incom. solver is included at the last step by implicitly solving the Helmholtz equation. /Yan On 22 May 2015, at 16:00, Sherwin, Spencer J <s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk>> wrote: HI Zhicheng, Since Yan is very familiar with this bit of code I suggest that he liase with you on this issue. Cheers Spencer. On 22 May 2015, at 09:10, Zhicheng Wang <wangzhicheng09@gmail.com<mailto:wangzhicheng09@gmail.com>> wrote: 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); _______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052
hi, Yan, thank you very much for reply. the reason why i'd like to use explicit scheme for the additional stabilized term is that the artificial viscosity is variable with quadratue points and time. if this term is discretesized implicitly, the Helmohltz matrix will be assembled at each time step, which is a time consuming procedure. the code i attached is computing the diffusion term with variable viscosity for the case of 3DHelmholtz1D case. On 2015年5月24日周日 00:58 Bao, Yan <y.bao@imperial.ac.uk> wrote:
H Zhicheng,
I can’t get much useful information from your attached code, so I don’t know what’s your problem or mistakes happened in your implementation.
Also, I don’t really understand why do you use explicit scheme for your stabilized term, As far as I know, the viscous term in Incom. solver is included at the last step by implicitly solving the Helmholtz equation.
/Yan
On 22 May 2015, at 16:00, Sherwin, Spencer J <s.sherwin@imperial.ac.uk <mailto:s.sherwin@imperial.ac.uk>> wrote:
HI Zhicheng,
Since Yan is very familiar with this bit of code I suggest that he liase with you on this issue.
Cheers Spencer.
On 22 May 2015, at 09:10, Zhicheng Wang <wangzhicheng09@gmail.com<mailto: wangzhicheng09@gmail.com>> wrote:
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); _______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk<mailto:Nektar-users@imperial.ac.uk> https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
Spencer Sherwin McLaren Racing/Royal Academy of Engineering Research Chair, Professor of Computational Fluid Mechanics, Department of Aeronautics, Imperial College London South Kensington Campus London SW7 2AZ
s.sherwin@imperial.ac.uk<mailto:s.sherwin@imperial.ac.uk> +44 (0) 20 759 45052
_______________________________________________ Nektar-users mailing list Nektar-users@imperial.ac.uk https://mailman.ic.ac.uk/mailman/listinfo/nektar-users
participants (3)
-
Bao, Yan
-
Sherwin, Spencer J
-
Zhicheng Wang