137 lines
9.0 KiB
C
137 lines
9.0 KiB
C
#include "macrodef.h"
|
|
#include "share_func.h"
|
|
|
|
/*
|
|
* kodis_sh — Kreiss-Oliger dissipation on shell patches.
|
|
* Same stencil coefficients as Cartesian kodis. Uses symmetry_stbd.
|
|
*/
|
|
extern "C" void kodis_sh_(const int ex[3],
|
|
const double *X, const double *Y, const double *Z,
|
|
const double *f, double *f_rhs,
|
|
const double SoAi[2],
|
|
int Symmetry, double eps, int sst)
|
|
{
|
|
(void)sst;
|
|
const double ZEO=0.0;
|
|
const int ex1=ex[0], ex2=ex[1], ex3=ex[2];
|
|
const double dX=X[1]-X[0], dY=Y[1]-Y[0], dZ=Z[1]-Z[0];
|
|
const int imaxF=ex1, jmaxF=ex2, kmaxF=ex3;
|
|
const double SoA[2]={SoAi[0],SoAi[1]};
|
|
|
|
#if (ghost_width == 2)
|
|
{
|
|
const int ord=2, r=2;
|
|
const double cof=16.0, F4=4.0, F6=6.0;
|
|
const int NO_SYMM=0, OCTANT=2;
|
|
int iminF=1,jminF=1,kminF=1;
|
|
if(Symmetry==OCTANT&&fabs(X[0])<dX)iminF=-1;
|
|
if(Symmetry==OCTANT&&fabs(Y[0])<dY)jminF=-1;
|
|
if((sst==2||sst==4)&&fabs(Y[0])<dY)jminF=-1;
|
|
|
|
const size_t nx=(size_t)ex1+2*ord,ny=(size_t)ex2+2*ord,nz=(size_t)ex3,fh_size=nx*ny*nz;
|
|
double *fh=(double*)malloc(fh_size*sizeof(double));if(!fh)return;
|
|
symmetry_stbd(ord,ex,f,fh,SoA);
|
|
|
|
const int i0_lo=(iminF+1>0)?iminF+1:0,j0_lo=(jminF+1>0)?jminF+1:0,k0_lo=2;
|
|
const int i0_hi=imaxF-3,j0_hi=jmaxF-3,k0_hi=kmaxF-3;
|
|
if(!(i0_lo>i0_hi||j0_lo>j0_hi||k0_lo>k0_hi)){
|
|
for(int k0=k0_lo;k0<=k0_hi;++k0){const int kF=k0+1;
|
|
for(int j0=j0_lo;j0<=j0_hi;++j0){const int jF=j0+1;
|
|
for(int i0=i0_lo;i0<=i0_hi;++i0){const int iF=i0+1;const size_t p=idx_ex(i0,j0,k0,ex);
|
|
const double Dx=((fh[idx_fh_stbd(iF-2,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+2,jF,kF,ord,ex)])-F4*(fh[idx_fh_stbd(iF-1,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+1,jF,kF,ord,ex)])+F6*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dX;
|
|
const double Dy=((fh[idx_fh_stbd(iF,jF-2,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+2,kF,ord,ex)])-F4*(fh[idx_fh_stbd(iF,jF-1,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+1,kF,ord,ex)])+F6*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dY;
|
|
const double Dz=((fh[idx_fh_stbd(iF,jF,kF-2,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+2,ord,ex)])-F4*(fh[idx_fh_stbd(iF,jF,kF-1,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+1,ord,ex)])+F6*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dZ;
|
|
f_rhs[p]-=(eps/cof)*(Dx+Dy+Dz);
|
|
}}}
|
|
}
|
|
free(fh);return;
|
|
}
|
|
#elif (ghost_width == 3)
|
|
{
|
|
const int ord=3, r=3;
|
|
const double cof=64.0,SIX=6.0,FIT=15.0,TWT=20.0;
|
|
const int NO_SYMM=0,OCTANT=2;
|
|
int iminF=1,jminF=1,kminF=1;
|
|
if(Symmetry==OCTANT&&fabs(X[0])<dX)iminF=-2;
|
|
if(Symmetry==OCTANT&&fabs(Y[0])<dY)jminF=-2;
|
|
if((sst==2||sst==4)&&fabs(Y[0])<dY)jminF=-2;
|
|
|
|
const size_t nx=(size_t)ex1+2*ord,ny=(size_t)ex2+2*ord,nz=(size_t)ex3,fh_size=nx*ny*nz;
|
|
double *fh=(double*)malloc(fh_size*sizeof(double));if(!fh)return;
|
|
symmetry_stbd(ord,ex,f,fh,SoA);
|
|
|
|
const int i0_lo=(iminF+2>0)?iminF+2:0,j0_lo=(jminF+2>0)?jminF+2:0,k0_lo=3;
|
|
const int i0_hi=imaxF-4,j0_hi=jmaxF-4,k0_hi=kmaxF-4;
|
|
if(!(i0_lo>i0_hi||j0_lo>j0_hi||k0_lo>k0_hi)){
|
|
for(int k0=k0_lo;k0<=k0_hi;++k0){const int kF=k0+1;
|
|
for(int j0=j0_lo;j0<=j0_hi;++j0){const int jF=j0+1;
|
|
for(int i0=i0_lo;i0<=i0_hi;++i0){const int iF=i0+1;const size_t p=idx_ex(i0,j0,k0,ex);
|
|
const double Dx=((fh[idx_fh_stbd(iF-3,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+3,jF,kF,ord,ex)])-SIX*(fh[idx_fh_stbd(iF-2,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+2,jF,kF,ord,ex)])+FIT*(fh[idx_fh_stbd(iF-1,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+1,jF,kF,ord,ex)])-TWT*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dX;
|
|
const double Dy=((fh[idx_fh_stbd(iF,jF-3,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+3,kF,ord,ex)])-SIX*(fh[idx_fh_stbd(iF,jF-2,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+2,kF,ord,ex)])+FIT*(fh[idx_fh_stbd(iF,jF-1,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+1,kF,ord,ex)])-TWT*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dY;
|
|
const double Dz=((fh[idx_fh_stbd(iF,jF,kF-3,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+3,ord,ex)])-SIX*(fh[idx_fh_stbd(iF,jF,kF-2,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+2,ord,ex)])+FIT*(fh[idx_fh_stbd(iF,jF,kF-1,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+1,ord,ex)])-TWT*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dZ;
|
|
f_rhs[p]+=(eps/cof)*(Dx+Dy+Dz);
|
|
}}}
|
|
}
|
|
free(fh);return;
|
|
}
|
|
#elif (ghost_width == 4)
|
|
{
|
|
const int ord=4, r=4;
|
|
const double cof=256.0,F8=8.0,F28=28.0,F56=56.0,F70=70.0;
|
|
const int NO_SYMM=0,OCTANT=2;
|
|
int iminF=1,jminF=1,kminF=1;
|
|
if(Symmetry==OCTANT&&fabs(X[0])<dX)iminF=-3;
|
|
if(Symmetry==OCTANT&&fabs(Y[0])<dY)jminF=-3;
|
|
if((sst==2||sst==4)&&fabs(Y[0])<dY)jminF=-3;
|
|
|
|
const size_t nx=(size_t)ex1+2*ord,ny=(size_t)ex2+2*ord,nz=(size_t)ex3,fh_size=nx*ny*nz;
|
|
double *fh=(double*)malloc(fh_size*sizeof(double));if(!fh)return;
|
|
symmetry_stbd(ord,ex,f,fh,SoA);
|
|
|
|
const int i0_lo=(iminF+3>0)?iminF+3:0,j0_lo=(jminF+3>0)?jminF+3:0,k0_lo=4;
|
|
const int i0_hi=imaxF-5,j0_hi=jmaxF-5,k0_hi=kmaxF-5;
|
|
if(!(i0_lo>i0_hi||j0_lo>j0_hi||k0_lo>k0_hi)){
|
|
for(int k0=k0_lo;k0<=k0_hi;++k0){const int kF=k0+1;
|
|
for(int j0=j0_lo;j0<=j0_hi;++j0){const int jF=j0+1;
|
|
for(int i0=i0_lo;i0<=i0_hi;++i0){const int iF=i0+1;const size_t p=idx_ex(i0,j0,k0,ex);
|
|
const double Dx=((fh[idx_fh_stbd(iF-4,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+4,jF,kF,ord,ex)])-F8*(fh[idx_fh_stbd(iF-3,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+3,jF,kF,ord,ex)])+F28*(fh[idx_fh_stbd(iF-2,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+2,jF,kF,ord,ex)])-F56*(fh[idx_fh_stbd(iF-1,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+1,jF,kF,ord,ex)])+F70*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dX;
|
|
const double Dy=((fh[idx_fh_stbd(iF,jF-4,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+4,kF,ord,ex)])-F8*(fh[idx_fh_stbd(iF,jF-3,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+3,kF,ord,ex)])+F28*(fh[idx_fh_stbd(iF,jF-2,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+2,kF,ord,ex)])-F56*(fh[idx_fh_stbd(iF,jF-1,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+1,kF,ord,ex)])+F70*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dY;
|
|
const double Dz=((fh[idx_fh_stbd(iF,jF,kF-4,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+4,ord,ex)])-F8*(fh[idx_fh_stbd(iF,jF,kF-3,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+3,ord,ex)])+F28*(fh[idx_fh_stbd(iF,jF,kF-2,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+2,ord,ex)])-F56*(fh[idx_fh_stbd(iF,jF,kF-1,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+1,ord,ex)])+F70*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dZ;
|
|
f_rhs[p]-=(eps/cof)*(Dx+Dy+Dz);
|
|
}}}
|
|
}
|
|
free(fh);return;
|
|
}
|
|
#elif (ghost_width == 5)
|
|
{
|
|
const int ord=5, r=5;
|
|
const double cof=1024.0,F10=10.0,F45k=45.0,F120=120.0,F210=210.0,F252=252.0;
|
|
const int NO_SYMM=0,OCTANT=2;
|
|
int iminF=1,jminF=1,kminF=1;
|
|
if(Symmetry==OCTANT&&fabs(X[0])<dX)iminF=-4;
|
|
if(Symmetry==OCTANT&&fabs(Y[0])<dY)jminF=-4;
|
|
if((sst==2||sst==4)&&fabs(Y[0])<dY)jminF=-4;
|
|
|
|
const size_t nx=(size_t)ex1+2*ord,ny=(size_t)ex2+2*ord,nz=(size_t)ex3,fh_size=nx*ny*nz;
|
|
double *fh=(double*)malloc(fh_size*sizeof(double));if(!fh)return;
|
|
symmetry_stbd(ord,ex,f,fh,SoA);
|
|
|
|
const int i0_lo=(iminF+4>0)?iminF+4:0,j0_lo=(jminF+4>0)?jminF+4:0,k0_lo=5;
|
|
const int i0_hi=imaxF-6,j0_hi=jmaxF-6,k0_hi=kmaxF-6;
|
|
if(!(i0_lo>i0_hi||j0_lo>j0_hi||k0_lo>k0_hi)){
|
|
for(int k0=k0_lo;k0<=k0_hi;++k0){const int kF=k0+1;
|
|
for(int j0=j0_lo;j0<=j0_hi;++j0){const int jF=j0+1;
|
|
for(int i0=i0_lo;i0<=i0_hi;++i0){const int iF=i0+1;const size_t p=idx_ex(i0,j0,k0,ex);
|
|
const double Dx=((fh[idx_fh_stbd(iF-5,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+5,jF,kF,ord,ex)])-F10*(fh[idx_fh_stbd(iF-4,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+4,jF,kF,ord,ex)])+F45k*(fh[idx_fh_stbd(iF-3,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+3,jF,kF,ord,ex)])-F120*(fh[idx_fh_stbd(iF-2,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+2,jF,kF,ord,ex)])+F210*(fh[idx_fh_stbd(iF-1,jF,kF,ord,ex)]+fh[idx_fh_stbd(iF+1,jF,kF,ord,ex)])-F252*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dX;
|
|
const double Dy=((fh[idx_fh_stbd(iF,jF-5,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+5,kF,ord,ex)])-F10*(fh[idx_fh_stbd(iF,jF-4,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+4,kF,ord,ex)])+F45k*(fh[idx_fh_stbd(iF,jF-3,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+3,kF,ord,ex)])-F120*(fh[idx_fh_stbd(iF,jF-2,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+2,kF,ord,ex)])+F210*(fh[idx_fh_stbd(iF,jF-1,kF,ord,ex)]+fh[idx_fh_stbd(iF,jF+1,kF,ord,ex)])-F252*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dY;
|
|
const double Dz=((fh[idx_fh_stbd(iF,jF,kF-5,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+5,ord,ex)])-F10*(fh[idx_fh_stbd(iF,jF,kF-4,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+4,ord,ex)])+F45k*(fh[idx_fh_stbd(iF,jF,kF-3,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+3,ord,ex)])-F120*(fh[idx_fh_stbd(iF,jF,kF-2,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+2,ord,ex)])+F210*(fh[idx_fh_stbd(iF,jF,kF-1,ord,ex)]+fh[idx_fh_stbd(iF,jF,kF+1,ord,ex)])-F252*fh[idx_fh_stbd(iF,jF,kF,ord,ex)])/dZ;
|
|
f_rhs[p]+=(eps/cof)*(Dx+Dy+Dz);
|
|
}}}
|
|
}
|
|
free(fh);return;
|
|
}
|
|
#else
|
|
#error "kodiss_sh_c.C: unsupported ghost_width"
|
|
#endif
|
|
}
|