Problem about viscous flux calculation at SPs

[YaojieYu]

Hi Jacob,

I am a bit confused about the calculation of viscous flux at SPs. If the motion flag is active, the divergence would be calculated by Liang-Miyaji's method and the viscous fluxes don't need to transform back to reference domain, just as the invicid flux function do. The following is the original code.

solver::calcViscousFlux_spts(void)
{
 double tempF[3][5];
#pragma omp parallel for collapse(2)
 for (uint spt = 0; spt < nSpts; spt++) {
   for (uint e = 0; e < nEles; e++) {
     for (uint dim = 0; dim < nDims; dim++)
       for (uint k = 0; k < nFields; k++)
         tempDU(dim,k) = dU_spts(dim,spt,e,k);

     if (params->equation == NAVIER_STOKES)
       viscousFlux(&U_spts(spt,e,0), tempDU, tempF, params);
     else
       viscousFluxAD(tempDU, tempF, params);
     /* Add physical inviscid flux at spts */
     for (uint dim = 0; dim < nDims; dim++)
       for (uint k = 0; k < nFields; k++)
         tempF[dim][k] += F_spts(dim,spt,e,k);

     /* --- Transform back to reference domain --- */
     for (uint dim1 = 0; dim1 < nDims; dim1++)
       for (uint k = 0; k < nFields; k++)
         F_spts(dim1,spt,e,k) = 0.;

     for (uint dim1 = 0; dim1 < nDims; dim1++) {
       for (uint k = 0; k < nFields; k++) {
         for (uint dim2 = 0; dim2 < nDims; dim2++) {
           F_spts(dim1,spt,e,k) += JGinv_spts(dim1,spt,e,dim2)*tempF[dim2][k];
         }
       }
     }
   }
 }
}

I added a motion branch and designed a special test (moveAx= moveAy=moveFx=moveFy=0) to check it. In this test, the residual was the same as static one. After fix:

void solver::calcViscousFlux_spts(void)
{
	double tempF[3][5];
	#pragma omp parallel for collapse(2)
	for (uint spt = 0; spt < nSpts; spt++) {
		for (uint e = 0; e < nEles; e++) {
			for (uint dim = 0; dim < nDims; dim++)
				for (uint k = 0; k < nFields; k++)
					tempDU(dim,k) = dU_spts(dim,spt,e,k);

			if (params->equation == NAVIER_STOKES)
				viscousFlux(&U_spts(spt,e,0), tempDU, tempF, params);
			else
				viscousFluxAD(tempDU, tempF, params);

			/* Add physical inviscid flux at spts */
			for (uint dim = 0; dim < nDims; dim++)
				for (uint k = 0; k < nFields; k++)
					tempF[dim][k] += F_spts(dim,spt,e,k);

			/* --- Transform back to reference domain --- */
			for (uint dim1 = 0; dim1 < nDims; dim1++)
				for (uint k = 0; k < nFields; k++)
					F_spts(dim1,spt,e,k) = 0;
										   
			if (params->motion) {
				for (uint dim1 = 0; dim1 < nDims; dim1++) 
					for (uint k = 0; k < nFields; k++) 
							F_spts(dim1,spt,e,k) =tempF[dim1][k];					
					}
				
			else {
				for (uint dim1 = 0; dim1 < nDims; dim1++) 
					for (uint k = 0; k < nFields; k++) 
						for (uint dim2 = 0; dim2 < nDims; dim2++) 
							F_spts(dim1,spt,e,k) += JGinv_spts(dim1,spt,e,dim2)*tempF[dim2][k];					 					
				}                    
			}
		}
}

So, I guess maybe something wrong with the function. Please check it.