1 #ifndef __NAVIERSTOKES_RHS_H__ 2 #define __NAVIERSTOKES_RHS_H__ 32 #define CONSTANT_CFL 1 41 template<
typename Gr
idT,
typename StateT,
typename OperatorT>
54 using rhs_base_t::globalPtr;
58 rhs() : rhs_base_t(), useAlternateRHS(false), exchangeFlux(false),
59 exchangeDV(false), artificialDissipation(false),useWENO(false)
69 int Initialize(GridType &inGrid, StateType &inState, StateType &inParam,
70 OperatorType &inOperator)
76 std::cout <<
"NavierStokesRHS:Fatal error: Init failed." << std::endl;
81 std::cout <<
"NavierStokesRHS:Fatal error: Set Grid failed." << std::endl;
85 if(SetParamState(inParam)){
86 std::cout <<
"NavierStokesRHS:Fatal error: Set params failed." << std::endl;
93 globalPtr->FunctionEntry(
"InitializeState");
96 if(InitializeState(inState)){
97 std::cout <<
"NavierStokesRHS:Fatal error: Initialize state failed." << std::endl;
103 globalPtr->FunctionExit(
"InitializeState");
104 globalPtr->FunctionEntry(
"InitMinSpacing");
113 globalPtr->FunctionExit(
"InitMinSpacing");
114 globalPtr->FunctionEntry(
"HaloSetup");
118 if(HaloSetup(inGrid.Halo())){
119 std::cout <<
"NavierStokesRHS:Fatal error: Halo setup failed." << std::endl;
126 globalPtr->FunctionExit(
"HaloSetup");
127 globalPtr->FunctionEntry(
"SetupOperators");
131 if(SetupOperators(inOperator)){
132 std::cout <<
"NavierStokesRHS:Fatal error: Operator setup failed." << std::endl;
138 globalPtr->FunctionExit(
"SetupOperators");
150 StateType &rhsState,
int myThreadId = 0)
158 errorCode = SetHalo(inGrid.Halo());
168 errorCode = SetState(inState);
177 errorCode = SetRHSState(rhsState);
195 &fortranBufferInterval(fortranBufferIntervals[threadId]);
202 globalPtr->FunctionEntry(
"ZeroRHS");
207 for(
int iEquation = 0;iEquation < numEquations;iEquation++){
210 &fortranBufferInterval[0],&zero,myRHSBuffers[iEquation]);
217 globalPtr->FunctionExit(
"ZeroRHS");
227 PrepareBuffers(threadId);
231 if(ComputeDV(threadId))
238 if(ComputeTV(threadId))
241 if(refRe > 0 || artificialDissipation){
243 if(ComputeVelGrad(threadId))
246 if(artificialDissipation){
247 if(ComputeVelDiv(threadId))
251 if(ComputeStressTensor(threadId))
254 if(ComputeTemperatureGrad(threadId))
258 if(ComputeHeatFlux(threadId))
263 if(ComputeViscidEnergyFlux(threadId))
276 for(
int iDim = 0;iDim < numDim;iDim++){
279 if(InviscidFlux(iDim,threadId))
285 ViscidFlux(iDim,threadId);
290 if(ExchangeState(threadId))
294 if(ExchangeBuffer(fluxMessageId,numEquations,inviscidFluxBuffer,threadId))
297 int wenoStatus = ApplyWENO(iDim,threadId);
304 std::ostringstream statusMessage;
305 statusMessage <<
"ApplyWENO returned error code (" 306 << wenoStatus <<
"), exiting." 308 globalPtr->StdOut(statusMessage.str());
317 AccumulateToRHS(m1,dFluxBuffer,threadId);
321 if(ExchangeBuffer(fluxMessageId,numEquations,viscidFluxBuffer,threadId))
324 if(ApplyOperator(iDim,numEquations,viscidFluxBuffer,dFluxBuffer,threadId))
328 AccumulateToRHS(m2,dFluxBuffer,threadId);
332 AccumulateViscousFluxes(-1.0,viscidFluxBuffers,inviscidFluxBuffers,threadId,!exchangeFlux);
335 if(ExchangeBuffer(fluxMessageId,numEquations,inviscidFluxBuffer,threadId))
354 if(ApplyOperator(iDim,numEquations,inviscidFluxBuffer,dFluxBuffer,threadId))
375 AccumulateToRHS(m1,dFluxBuffer,threadId);
397 GridScaleRHS(threadId);
425 globalPtr->FunctionEntry(
"HandleBCs");
428 int returnCode = HandleBoundaryConditions(threadId);
433 globalPtr->FunctionExit(
"HandleBCs");
437 return(returnCode+100);
456 returnCode = HandleSources(threadId);
460 return(returnCode+200);
478 globalPtr->FunctionEntry(
"ApplyWENO");
482 int winWidth = coeffsWENO.hiAll - coeffsWENO.loAll + 1;
484 double fluxWindow[winWidth*numEquations];
485 double projFluxWindow[winWidth*numEquations];
486 double projFluxReconst[numEquations];
488 double stateL[numEquations];
489 double stateR[numEquations];
490 double deltaState[numEquations];
491 double projDeltaState[numEquations];
492 double fluxL[numEquations];
493 double fluxR[numEquations];
495 double tmat[numEquations*numEquations];
496 double tinv[numEquations*numEquations];
497 double lambda[numEquations*numEquations];
498 double lambdaL[numEquations*numEquations];
499 double lambdaR[numEquations*numEquations];
500 double normVec[numDim];
501 for (
int i=0; i<numDim; i++) {
502 if (i == iDim) normVec[i] = 1;
506 const std::vector<pcpp::IndexIntervalType>
507 &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
514 int jDim = (iDim + 1)%numDim;
515 int kDim = (iDim + 2)%numDim;
517 size_t js = regionInterval[jDim].first;
518 size_t je = regionInterval[jDim].second;
524 ks = regionInterval[kDim].first;
525 ke = regionInterval[kDim].second;
528 for (
size_t k=ks; k<=ke; k++) {
534 opInterval[kDim].first = k;
535 opInterval[kDim].second = k;
538 for (
size_t j=js; j<=je; j++) {
539 opInterval[jDim].first = j;
540 opInterval[jDim].second = j;
542 size_t is = regionInterval[iDim].first;
543 size_t ie = regionInterval[iDim].second;
545 opInterval[iDim].first = is + coeffsWENO.loAll - 1;
546 opInterval[iDim].second = ie + coeffsWENO.hiAll;
548 std::vector<size_t> indices;
552 for (
int loc=0; loc<winWidth; loc++) {
553 for (
int eqn=0; eqn<numEquations; eqn++) {
555 fluxWindow[eqn*winWidth + loc] = inviscidFluxBuffer[eqn*
numPointsBuffer + indices[loc]];
560 for (
size_t ii=-coeffsWENO.loAll; ii<indices.size()-coeffsWENO.hiAll; ii++) {
562 for (
int eqn=0; eqn<numEquations; eqn++) {
563 stateL[eqn] = myStateBuffers[eqn][indices[ii]];
564 stateR[eqn] = myStateBuffers[eqn][indices[ii+1]];
565 deltaState[eqn] = stateR[eqn] - stateL[eqn];
571 double gamma = eosPtr->GetGamma();
573 (&numDim, &stateL[0], &stateL[0], &gamma, &normVec[0], &tmat[0], &tinv[0], &lambdaL[0]);
575 (&numDim, &stateR[0], &stateR[0], &gamma, &normVec[0], &tmat[0], &tinv[0], &lambdaR[0]);
580 (&numDim, &stateL[0], &stateR[0], &gamma, &normVec[0], &tmat[0], &tinv[0], &lambda[0]);
582 WENO::Project(numEquations, winWidth, &tinv[0], &fluxWindow[0], &projFluxWindow[0]);
585 WENO::Project(numEquations, 1, &tinv[0], &deltaState[0], &projDeltaState[0]);
588 for (
int eqn=0; eqn<numEquations; eqn++) {
589 int diag = eqn*(numEquations+1);
590 int group = (1 + (int) copysign(1.0, lambda[diag]))/2;
592 int startLoc = coeffsWENO.LoGroup(group) - coeffsWENO.loAll;
593 WENO::ReconstPointVal(&projFluxWindow[eqn*winWidth + startLoc], coeffsWENO, group, projFluxReconst[eqn]);
596 double lambdaAbs = std::abs(lambda[diag]);
597 double entropyFix = (eta > lambdaAbs ? 0.5*(eta - lambdaAbs) : 0);
598 entropyFix *= projDeltaState[eqn];
607 projFluxReconst[eqn] -= entropyFix;
610 WENO::Project(numEquations, 1, &tmat[0], &projFluxReconst[0], &fluxR[0]);
612 for (
int eqn=0; eqn<numEquations && ii > -coeffsWENO.loAll; eqn++) {
613 dFluxBuffer[eqn*
numPointsBuffer + indices[ii]] = fluxR[eqn] - fluxL[eqn];
617 for (
int loc=1; loc<winWidth; loc++) {
618 for (
int eqn=0; eqn<numEquations; eqn++) {
619 int currIndex = eqn*winWidth + loc;
620 int newIndex = eqn*winWidth + (loc-1);
621 fluxWindow[newIndex] = fluxWindow[currIndex];
626 for (
int eqn=0; ((eqn<numEquations) && (ii < (indices.size()-coeffsWENO.hiAll-1))); eqn++) {
628 fluxWindow[eqn*winWidth + (winWidth-1)] = inviscidFluxBuffer[eqn*
numPointsBuffer + indices[ii+coeffsWENO.hiAll+1]];
631 for (
int eqn=0; eqn<numEquations; eqn++) {
632 fluxL[eqn] = fluxR[eqn];
643 globalPtr->FunctionExit(
"ApplyWENO");
654 size_t *dOpInterval = &fortranPartitionBufferIntervals[threadId][0];
655 std::vector<double> &
gridMetric(gridPtr->Metric());
662 globalPtr->FunctionEntry(
"ArtificialDissipation");
666 double dilRamp = 10.0;
667 double sigmaDil = inputSigmaDil;
669 double dilCutoff = -1.5;
676 for(
int iDim = 0;iDim < numDim;iDim++){
678 int dissDir = iDim + 1;
687 for(
int iEquation= 0;iEquation < numEquations;iEquation++){
690 double *stateBuf = myStateBuffers[iEquation];
698 ExchangeBuffer(fluxMessageId,numEquations,inviscidFluxBuffer,threadId);
703 ApplyDissipationOperator(iDim,numEquations,inviscidFluxBuffer,dFluxBuffer,
false,threadId);
705 for(
int iEquation= 0;iEquation < numEquations;iEquation++){
714 (&numDim,&
numPointsBuffer,&bufferSizes[0],dOpInterval,alphaDissipation,
719 ExchangeBuffer(fluxMessageId,numEquations,inviscidFluxBuffer,threadId);
724 ApplyDissipationOperator(iDim,numEquations,inviscidFluxBuffer,dFluxBuffer,
true,threadId);
727 for(
int iEquation= 0;iEquation < numEquations;iEquation++){
730 double *rhsPtr = myRHSBuffers[iEquation];
743 globalPtr->FunctionExit(
"ArtificialDissipation");
757 if(artificialDissipation){
758 returnCode = ArtificialDissipation(myThreadId);
776 else if (!subRegionsPtr || !domainBoundariesPtr)
779 std::vector<BoundaryType> &domainBoundaries(*domainBoundariesPtr);
780 std::vector<BCType> &domainBCs(*domainBCsPtr);
781 std::vector<GridRegionType> &gridRegions(*subRegionsPtr);
782 std::vector<double> &
gridMetric(gridPtr->Metric());
785 int numBCs = domainBCs.size();
786 int numGridRegions = gridRegions.size();
787 int numDomainBoundaries = domainBoundaries.size();
791 if(numDomainBoundaries == 0)
794 typename std::vector<BoundaryType>::iterator dbIt = domainBoundaries.begin();
798 double gamma = eosPtr->GetGamma();
799 gasParams[0] = refRe;
800 gasParams[1] = gamma;
801 gasParams[2] = specGasConst;
802 gasParams[3] = temperatureScale;
803 gasParams[4] = refSndSpd;
804 gasParams[5] = std::max(gamma/refPr,5.0/3.0);
806 while(dbIt != domainBoundaries.end()){
808 BoundaryType &domainBoundary(*dbIt++);
809 GridRegionType &gridRegion(gridRegions[domainBoundary.
GridRegionID()]);
811 &threadPartBufferIndices(gridRegion.threadPartitionBufferIndices[threadId]);
812 size_t numBoundaryPoints = threadPartBufferIndices.size();
814 if(numBoundaryPoints > 0){
816 BCType &boundaryCondition(domainBoundary.
BC());
818 int bcType = boundaryCondition.BCType();
823 int normalDir = gridRegion.normalDirection;
825 std::vector<size_t> regionSizes(regionInterval.
Sizes());
826 size_t numPointsRegion = regionInterval.
NNodes();
832 StateType &bcParamState(boundaryCondition.Params());
833 double *
bcParams = bcParamState.template GetFieldBuffer<double>(
"Numbers");
834 int *bcFlags = bcParamState.template GetFieldBuffer<int>(
"Flags");
836 const std::string &bcName(boundaryCondition.BCName());
837 const std::string ®ionName(gridRegion.regionName);
838 const std::string &boundaryName(domainBoundary.
Name());
857 regionOpBufferInterval,bcParams,bcFlags,
858 numEquations,myStateBuffers,targetStateBuffers,myRHSBuffers);
863 regionOpBufferInterval,bcParams,bcFlags,iMask,
holeMask,
864 numEquations,myStateBuffers,targetStateBuffers,myRHSBuffers);
871 bcParams[2] = boundaryWeight;
875 &numPointsRegion,&numBoundaryPoints,&threadPartBufferIndices[0],
877 rhoPtr,rhoVPtr,rhoEPtr,viscidFluxBuffer,&
numScalar,scalarPtr,
878 rhoRHSPtr,rhoVRHSPtr,rhoERHSPtr,scalarRHSPtr,
879 rhoTargetPtr,rhoVTargetPtr,rhoETargetPtr,scalarTargetPtr);
887 std::ostringstream errMsg;
888 errMsg <<
"NavierStokesRHS::HandleBoundaryConditions: Error SAT_NOSLIP_ISOTHERMAL " 889 <<
"is not implemented for scalar transport. Aborting." << std::endl;
891 globalPtr->StdOut(errMsg.str());
893 std::cout << errMsg.str();
897 bcParams[3] = boundaryWeight;
902 &numPointsRegion,&numBoundaryPoints,&threadPartBufferIndices[0],&
gridType,
903 &gridMetric[0],&gridJacobian[0],&bcParams[0],&gasParams[0],rhoPtr,rhoVPtr,
904 rhoEPtr,&
numScalar,scalarPtr,rhoRHSPtr,rhoVRHSPtr,rhoERHSPtr,
905 scalarRHSPtr,rhoTargetPtr,rhoVTargetPtr,rhoETargetPtr,
906 scalarTargetPtr,muPtr,lambdaPtr);
910 std::ostringstream errMsg;
911 errMsg <<
"NavierStokesRHS::HandleBoundaryConditions:Error: SAT_NOSLIP_ISOTHERMAL " 912 <<
"not a valid BC for inviscid (Re = 0) flows. Aborting." << std::endl;
914 globalPtr->StdOut(errMsg.str());
916 std::cout << errMsg.str();
923 bcParams[1] = boundaryWeight;
927 &numPointsRegion,&numBoundaryPoints,&threadPartBufferIndices[0],&
gridType,
928 &gridMetric[0],&gridJacobian[0],
bcParams,&gasParams[0],rhoPtr,rhoVPtr,
929 rhoEPtr,&
numScalar,scalarPtr,rhoRHSPtr,rhoVRHSPtr,rhoERHSPtr,
930 scalarRHSPtr,rhoTargetPtr,rhoVTargetPtr,rhoETargetPtr,scalarTargetPtr);
934 std::ostringstream errMsg;
935 errMsg <<
"NavierStokesRHS::HandleBoundaryConditions:Error: Could not resolve BC(" 936 << bcName <<
") = " << bcType <<
". Aborting." << std::endl;
938 globalPtr->StdOut(errMsg.str());
940 std::cout << errMsg.str() << std::endl;
959 size_t *dOpInterval = &fortranPartitionBufferIntervals[threadId][0];
968 globalPtr->FunctionEntry(
"ApplyOperator");
980 &uBuffer[compOffset],&duBuffer[compOffset]);
995 globalPtr->FunctionExit(
"ApplyOperator");
1003 bool split =
false,
int threadId = 0)
1007 size_t *dOpInterval = &fortranPartitionBufferIntervals[threadId][0];
1016 globalPtr->FunctionEntry(
"ApplyDissOperator");
1020 int numStencilsDiss = 0;
1021 int *stencilSizesDiss = NULL;
1022 int *stencilStartsDiss = NULL;
1023 int numStencilValuesDiss = 0;
1024 double *stencilWeightsDiss = NULL;
1025 int *stencilOffsetsDiss = NULL;
1026 int *stencilConnDiss = NULL;
1030 numStencilsDiss = artDissOperatorSplit.numStencils;
1031 stencilSizesDiss = artDissOperatorSplit.stencilSizes;
1032 stencilStartsDiss = artDissOperatorSplit.stencilStarts;
1033 stencilWeightsDiss = artDissOperatorSplit.stencilWeights;
1034 stencilOffsetsDiss = artDissOperatorSplit.stencilOffsets;
1035 stencilConnDiss = artDissConnSplit;
1037 numStencilsDiss = artDissOperator.numStencils;
1038 stencilSizesDiss = artDissOperator.stencilSizes;
1039 stencilStartsDiss = artDissOperator.stencilStarts;
1040 stencilWeightsDiss = artDissOperator.stencilWeights;
1041 stencilOffsetsDiss = artDissOperator.stencilOffsets;
1042 stencilConnDiss = artDissConn;
1052 &dDir,dOpInterval,&numStencilsDiss,stencilSizesDiss,
1053 stencilStartsDiss,&numStencilValuesDiss,stencilWeightsDiss,
1054 stencilOffsetsDiss,&stencilConnDiss[dirOffset],
1055 &uBuffer[compOffset],&duBuffer[compOffset]);
1070 globalPtr->FunctionExit(
"ApplyDissOperator");
1083 globalPtr->FunctionEntry(
"GridScaleRHS");
1087 size_t *dOpInterval = &fortranPartitionBufferIntervals[threadId][0];
1088 std::vector<double> &
gridJacobian(gridPtr->Jacobian());
1091 for(
int iEquation = 0;iEquation < numEquations;iEquation++){
1092 double *rhsPtr = myRHSBuffers[iEquation];
1108 globalPtr->FunctionExit(
"GridScaleRHS");
1120 globalPtr->FunctionEntry(
"AccumulateRHS");
1124 size_t *dOpInterval = &fortranPartitionBufferIntervals[threadId][0];
1128 for(
int iComp = 0;iComp < numEquations;iComp++){
1139 globalPtr->FunctionExit(
"AccumulateRHS");
1153 globalPtr->FunctionEntry(
"ComputeDV");
1158 const std::vector<pcpp::IndexIntervalType>
1159 &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
1162 std::vector<size_t> &fortranOpInterval(fortranPartitionBufferIntervals[threadId]);
1169 myStateBuffers,dvBuffers);
1172 returnCode = eosPtr->ComputePressureBuffer(regionInterval,
bufferSizes);
1173 returnCode = eosPtr->ComputeTemperatureBuffer(regionInterval,
bufferSizes);
1175 if(!returnCode && exchangeDV){
1176 returnCode = ExchangeDV(threadId);
1183 globalPtr->FunctionExit(
"ComputeDV");
1197 globalPtr->FunctionEntry(
"ComputeTV");
1202 const std::vector<pcpp::IndexIntervalType>
1203 &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
1206 std::vector<size_t> &fortranOpInterval(fortranPartitionBufferIntervals[threadId]);
1210 if(transportModel ==
POWER) {
1211 double heatCapacityCp = eosPtr->GetHeatCapacityCp();
1213 temperaturePtr,tvBuffers,
1214 beta,power,bulkViscFac,heatCapacityCp,refPr);
1216 std::ostringstream errMsg;
1217 errMsg <<
"NavierStokesRHS::RHS:Error: Unknown transport model in NavierStokesRHS. Aborting." 1220 globalPtr->StdOut(errMsg.str());
1221 globalPtr->FunctionExit(
"ComputeTV");
1223 std::cout << errMsg.str();
1233 globalPtr->FunctionExit(
"ComputeTV");
1248 globalPtr->FunctionEntry(
"InviscidFlux");
1252 size_t *fluxOpInterval = NULL;
1254 fluxOpInterval = &fortranPartitionBufferIntervals[threadId][0];
1256 fluxOpInterval = &fortranBufferIntervals[threadId][0];
1263 globalPtr->FunctionEntry(
"ComputeVHat");
1268 std::vector<double> &
gridMetric(gridPtr->Metric());
1269 double *
velocity = &velocityPtr[0];
1270 int fluxDim = iDim + 1;
1281 globalPtr->FunctionExit(
"ComputeVHat");
1282 globalPtr->FunctionEntry(
"Flux1D");
1291 &fluxDim,&
gridType,&gridMetric[0],rhoPtr,rhoVPtr,rhoEPtr,
1292 &
velHat[0],pressurePtr,inviscidFluxBuffer);
1314 globalPtr->FunctionExit(
"Flux1D");
1324 globalPtr->FunctionEntry(
"ScalarFlux1D");
1331 &
numScalar,scalarPtr,&velHat[0],inviscidFluxBuffers[numDim+2]);
1337 globalPtr->FunctionExit(
"ScalarFlux1D");
1347 globalPtr->FunctionExit(
"InviscidFlux");
1360 globalPtr->FunctionEntry(
"ComputeVelDiv");
1364 const std::vector<pcpp::IndexIntervalType>
1365 &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
1367 std::vector<size_t> &fortranOpInterval(fortranPartitionBufferIntervals[threadId]);
1368 std::vector<double> &
gridMetric(gridPtr->Metric());
1369 std::vector<double> &
gridJacobian(gridPtr->Jacobian());
1375 &fortranOpInterval[0],&zero,velDiv);
1386 globalPtr->FunctionExit(
"ComputeVelDiv");
1401 globalPtr->FunctionEntry(
"ComputeVelGrad");
1407 for(
int i=0;i<numDim;i++){
1408 ApplyOperator(i,numDim,velocityPtr,velGradBuffers[i],threadId);
1415 globalPtr->FunctionExit(
"ComputeVelGrad");
1428 globalPtr->FunctionEntry(
"ComputeTempGrad");
1433 const std::vector<pcpp::IndexIntervalType>
1434 &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
1436 std::vector<size_t> &fortranOpInterval(fortranPartitionBufferIntervals[threadId]);
1439 for(
int i=0;i<numDim;i++){
1440 ApplyOperator(i,1,temperaturePtr,temperatureGradBuffers[i],threadId);
1447 globalPtr->FunctionExit(
"ComputeTempGrad");
1460 globalPtr->FunctionEntry(
"ComputeStressTensor");
1464 const std::vector<pcpp::IndexIntervalType>
1465 &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
1467 std::vector<double> &
gridMetric(gridPtr->Metric());
1468 std::vector<double> &
gridJacobian(gridPtr->Jacobian());
1472 tvBuffers, velGradBuffers, tauBuffers);
1495 globalPtr->FunctionExit(
"ComputeStressTensor");
1509 globalPtr->FunctionEntry(
"ComputeHeatFlux");
1513 const std::vector<pcpp::IndexIntervalType>
1514 &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
1516 std::vector<double> &
gridMetric(gridPtr->Metric());
1517 std::vector<double> &
gridJacobian(gridPtr->Jacobian());
1522 tvBuffers, temperatureGradBuffers, heatFluxBuffers);
1546 globalPtr->FunctionExit(
"ComputeHeatFlux");
1559 globalPtr->FunctionEntry(
"ViscidEnergyFlux");
1563 const std::vector<pcpp::IndexIntervalType>
1564 &threadIntervals(gridPtr->ThreadPartitionBufferIntervals());
1566 const size_t *
opInterval = &fortranPartitionBufferIntervals[threadId][0];
1568 for(
int iDim = 0;iDim < numDim;iDim++){
1571 for(
int iVel = 0;iVel < numDim;iVel++){
1572 int tensorIndex = iDim*numDim + iVel;
1577 &velocityPtr[dirOffset],tauBuffers[tensorIndex],heatFluxBuffers[iDim]);
1586 globalPtr->FunctionExit(
"ViscidEnergyFlux");
1598 globalPtr->FunctionEntry(
"ViscidFlux");
1602 size_t *fluxOpInterval = NULL;
1604 fluxOpInterval = &fortranPartitionBufferIntervals[threadId][0];
1606 fluxOpInterval = &fortranBufferIntervals[threadId][0];
1608 std::vector<double> &
gridMetric(gridPtr->Metric());
1616 int fluxDim = iDim + 1;
1617 int tau1=iDim*numDim;
1648 globalPtr->FunctionExit(
"ViscidFlux");
1658 int threadId,
bool fullInterval =
false)
1661 size_t *bufferOpInterval = NULL;
1663 bufferOpInterval = &fortranPartitionBufferIntervals[threadId][0];
1665 bufferOpInterval = &fortranBufferIntervals[threadId][0];
1668 for(
int iComp = 1;iComp < numDim+2;iComp++){
1679 for(
int iDim = 0;iDim < numDV;iDim++)
1680 haloPtr->PackMessageBuffers(dvMessageId,iDim,1,dvBuffers[iDim]);
1688 haloPtr->PackMessageBuffers(messageId,0,numEquations,dataBuffer);
1694 haloPtr->UnPackMessageBuffers(messageId,0,numEquations,dataBuffer);
1699 haloPtr->SendMessage(messageId,gridNeighbors,*gridCommPtr);
1700 haloPtr->ReceiveMessage(messageId,gridNeighbors,*gridCommPtr);
1706 for(
int iDim = 0;iDim < numDV;iDim++)
1707 haloPtr->UnPackMessageBuffers(dvMessageId,iDim,1,dvBuffers[iDim]);
1714 haloPtr->SendMessage(dvMessageId,gridNeighbors,*gridCommPtr);
1715 haloPtr->ReceiveMessage(dvMessageId,gridNeighbors,*gridCommPtr);
1724 globalPtr->FunctionEntry(
"ExchangeDV");
1728 PackDVMessage(threadId);
1732 ExchangeDVMessage();
1735 UnpackDVMessage(threadId);
1740 globalPtr->FunctionExit(
"ExchangeDV");
1753 globalPtr->FunctionEntry(
"ExchangeBuffer");
1756 PackBufferMessage(messageId,numEquations,dataBuffer,threadId);
1761 ExchangeBufferMessage(messageId);
1766 UnpackBufferMessage(messageId,numEquations,dataBuffer,threadId);
1772 globalPtr->FunctionExit(
"ExchangeBuffer");
1783 return(ExchangeBuffer(fluxMessageId,numEquations,inviscidFluxBuffer,threadId));
1792 globalPtr->FunctionEntry(
"ExchangeState");
1798 haloPtr->PackMessageBuffers(fluxMessageId,0,1,rhoPtr);
1799 haloPtr->PackMessageBuffers(fluxMessageId,1,numDim,rhoVPtr);
1800 haloPtr->PackMessageBuffers(fluxMessageId,numDim+1,1,rhoEPtr);
1802 haloPtr->SendMessage(fluxMessageId,gridNeighbors,gridPtr->Communicator());
1803 haloPtr->ReceiveMessage(fluxMessageId,gridNeighbors,gridPtr->Communicator());
1805 haloPtr->UnPackMessageBuffers(fluxMessageId,0,1,rhoPtr);
1806 haloPtr->UnPackMessageBuffers(fluxMessageId,1,numDim,rhoVPtr);
1807 haloPtr->UnPackMessageBuffers(fluxMessageId,numDim+1,1,rhoEPtr);
1814 globalPtr->FunctionExit(
"ExchangeState");
1831 globalPtr->FunctionEntry(
"SetParamState");
1840 int paramHandle = paramState.GetDataIndex(
"refRe");
1841 if(paramHandle >= 0){
1843 const double *paramPtr = paramData.
Data<
double>();
1853 int paramHandle = paramState.GetDataIndex(
"refPr");
1854 if(paramHandle >= 0){
1856 const double *paramPtr = paramData.
Data<
double>();
1880 int paramHandle = paramState.GetDataIndex(
"nonDimensional");
1881 std::ostringstream statusMsg;
1882 statusMsg << (paramHandle >= 0 ?
"User chose" :
"Defaulting to");
1883 if(paramHandle >= 0){
1885 const int *paramPtr = paramData.
Data<
int>();
1886 nonDimensionalMode = *paramPtr;
1887 statusMsg << (nonDimensionalMode==1 ?
" non-dimensional " :
1888 nonDimensionalMode==2 ?
" PlasComCM-non-dimensional " :
1891 nonDimensionalMode = 1;
1892 statusMsg <<
" non-dimensional ";
1894 statusMsg <<
"mode." << std::endl;
1896 globalPtr->StdOut(statusMsg.str());
1898 std::cout << statusMsg.str();
1907 bool haveRefRho =
false;
1908 bool haveRefPressure =
false;
1909 bool haveRefTemperature =
false;
1911 int paramHandle = paramState.GetDataIndex(
"refRho");
1912 if(paramHandle >= 0){
1914 ¶mData(paramState.Field(paramHandle));
1915 const double *paramPtr = paramData.
Data<
double>();
1922 int paramHandle = paramState.GetDataIndex(
"refPressure");
1923 if(paramHandle >= 0){
1925 ¶mData(paramState.Field(paramHandle));
1926 const double *paramPtr = paramData.
Data<
double>();
1927 refPressure = *paramPtr;
1928 haveRefPressure =
true;
1933 int paramHandle = paramState.GetDataIndex(
"refTemperature");
1934 if(paramHandle >= 0){
1936 ¶mData(paramState.Field(paramHandle));
1937 const double *paramPtr = paramData.
Data<
double>();
1938 refTemperature = *paramPtr;
1939 haveRefTemperature =
true;
1943 if(!(haveRefRho && haveRefPressure && haveRefTemperature)){
1945 refPressure = 101325;
1946 refTemperature = 288.15;
1947 std::ostringstream statusMessage;
1948 statusMessage <<
"Warning: EOS reference state missing or incomplete in input. " << std::endl
1949 <<
"User must specify refRho, refPressure, and refTemperature. " << std::endl
1950 <<
"Setting to default values: " << std::endl
1951 <<
"\t refRho=" << refRho << std::endl
1952 <<
"\t refPressure=" << refPressure << std::endl
1953 <<
"\t refTemperature=" << refTemperature << std::endl
1956 globalPtr->StdOut(statusMessage.str());
1958 std::cout << statusMessage.str();
1963 int paramHandle = paramState.GetDataIndex(
"refLength");
1964 if(paramHandle >= 0){
1966 ¶mData(paramState.Field(paramHandle));
1967 const double *paramPtr = paramData.
Data<
double>();
1968 refLength = *paramPtr;
1971 std::ostringstream statusMessage;
1972 statusMessage <<
"Warning: Could not find refLength in input" << std::endl
1973 <<
"Setting to default value " << refLength << std::endl;
1975 globalPtr->StdOut(statusMessage.str());
1977 std::cout << statusMessage.str();
1987 if(gasModel ==
IDEAL) {
1990 gammaHandle = paramState.GetDataIndex(
"gamma");
1991 if(gammaHandle < 0) {
1992 std::ostringstream statusMessage;
1993 statusMessage <<
"Gamma does not exist in input parameter state." 1994 <<
" Exiting." << std::endl;
1996 globalPtr->StdOut(statusMessage.str());
1997 globalPtr->FunctionExit(
"SetParamState");
1999 std::cout << statusMessage.str();
2005 const double *gammaPtr = gammaData.
Data<
double>();
2007 std::ostringstream statusMessage;
2008 statusMessage <<
"Gamma not defined in input parameter state." 2009 <<
" Exiting." << std::endl;
2011 globalPtr->StdOut(statusMessage.str());
2012 globalPtr->FunctionExit(
"SetParamState");
2014 std::cout << statusMessage.str();
2019 eosPtr->SetGamma(*gammaPtr);
2021 std::ostringstream statusMessage;
2022 statusMessage <<
"Error: Unknown gas model selected." << gasModel
2023 <<
" Exiting." << std::endl;
2025 globalPtr->StdOut(statusMessage.str());
2026 globalPtr->FunctionExit(
"SetParamState");
2028 std::cout << statusMessage.str();
2034 if(transportModel ==
POWER) {
2037 powerHandle = paramState.GetDataIndex(
"power");
2038 if(powerHandle < 0) {
2039 std::ostringstream statusMessage;
2040 statusMessage <<
"Power does not exist in input parameter state." << std::endl
2041 <<
" Required for the Power Law transport model." << std::endl
2042 <<
" Exiting." << std::endl;
2044 globalPtr->StdOut(statusMessage.str());
2045 globalPtr->FunctionExit(
"SetParamState");
2047 std::cout << statusMessage.str();
2053 &powerData(paramState.Field(powerHandle));
2054 powerPtr = powerData.
Data<
double>();
2056 std::ostringstream statusMessage;
2057 statusMessage <<
"Power not defined in input parameter state." << std::endl
2058 <<
"Required for the Power Law transport model." << std::endl
2059 <<
"Example: power=2/3 for air." << std::endl
2060 <<
"Exiting." << std::endl;
2062 globalPtr->StdOut(statusMessage.str());
2063 globalPtr->FunctionExit(
"SetParamState");
2065 std::cout << statusMessage.str();
2072 betaHandle = paramState.GetDataIndex(
"beta");
2073 if(betaHandle < 0) {
2074 std::ostringstream statusMessage;
2075 statusMessage <<
"beta does not exist in input parameter state." << std::endl
2076 <<
"Required for the beta Law transport model." << std::endl
2077 <<
"Exiting." << std::endl;
2079 globalPtr->StdOut(statusMessage.str());
2080 globalPtr->FunctionExit(
"SetParamState");
2082 std::cout << statusMessage.str();
2088 &betaData(paramState.Field(betaHandle));
2089 betaPtr = betaData.
Data<
double>();
2091 std::ostringstream statusMessage;
2092 statusMessage <<
"beta not defined in input parameter state." << std::endl
2093 <<
"Required for the beta Law transport model." << std::endl
2094 <<
"Example: beta=4.093*10^-7 for air." << std::endl
2095 <<
"Exiting." << std::endl;
2097 globalPtr->StdOut(statusMessage.str());
2098 globalPtr->FunctionExit(
"SetParamState");
2100 std::cout << statusMessage.str();
2107 bulkViscFacHandle = paramState.GetDataIndex(
"bulkViscFac");
2108 if(bulkViscFacHandle < 0) {
2109 std::ostringstream statusMessage;
2110 statusMessage <<
"bulkViscFac does not exist in input parameter state." << std::endl
2111 <<
"Required for the bulkViscFac Law transport model."<< std::endl
2112 <<
"Exiting." << std::endl;
2114 globalPtr->StdOut(statusMessage.str());
2115 globalPtr->FunctionExit(
"SetParamState");
2117 std::cout << statusMessage.str();
2123 &bulkViscFacData(paramState.Field(bulkViscFacHandle));
2124 bulkViscFacPtr = bulkViscFacData.
Data<
double>();
2125 if(!bulkViscFacPtr) {
2126 std::ostringstream statusMessage;
2127 statusMessage <<
"bulkViscFac not defined in input parameter state." << std::endl
2128 <<
"Required for the bulkViscFac Law transport model."<< std::endl
2129 <<
"Example: bulkViscFac=0.6 for air." << std::endl
2130 <<
"Exiting." << std::endl;
2132 globalPtr->StdOut(statusMessage.str());
2133 globalPtr->FunctionExit(
"SetParamState");
2135 std::cout << statusMessage.str();
2139 bulkViscFac = *bulkViscFacPtr;
2142 std::ostringstream statusMessage;
2143 statusMessage <<
"Error: Unknown transport model selected." 2144 << transportModel <<
" Exiting." << std::endl;
2146 globalPtr->StdOut(statusMessage.str());
2147 globalPtr->FunctionExit(
"SetParamState");
2149 std::cout << statusMessage.str();
2158 if(gasModel ==
IDEAL) {
2159 double gamma = eosPtr->GetGamma();
2160 if(nonDimensionalMode == 1){
2162 refSpecGasConst = refPressure/refRho/refTemperature;
2163 refVelocity = sqrt(refPressure/refRho);
2164 refSndSpd = refVelocity*sqrt(gamma);
2165 refEnergy = refPressure/(gamma-1.0)/refRho;
2166 refCp = refSpecGasConst*gamma/(gamma-1.0);
2167 refCv = refCp-refSpecGasConst;
2170 pressureScale = refPressure;
2171 temperatureScale = refTemperature;
2172 velocityScale = refVelocity;
2174 }
else if(nonDimensionalMode == 2){
2176 specGasConst = (gamma-1)/gamma;
2177 refSpecGasConst = refPressure/refRho/refTemperature;
2178 refVelocity = sqrt(gamma*refPressure/refRho);
2179 refSndSpd = refVelocity;
2180 refEnergy = refPressure/(gamma-1.0)/refRho;
2181 refCp = refSpecGasConst*gamma/(gamma-1.0);
2182 refCv = refCp-refSpecGasConst;
2184 pressureScale = refRho*refSndSpd*refSndSpd;
2185 temperatureScale = (gamma-1)*refTemperature;
2186 velocityScale = refSndSpd;
2188 std::ostringstream statusMessage;
2189 statusMessage <<
"Enabling legacy PlasComCM nonDimensionalization mode." << std::endl;
2191 globalPtr->StdOut(statusMessage.str());
2193 std::cout << statusMessage.str();
2198 refVelocity = sqrt(refPressure/refRho);
2199 specGasConst = refPressure/refRho/refTemperature;
2200 refSndSpd = refVelocity*sqrt(gamma);
2201 refCp = specGasConst*gamma/(gamma-1.0);
2202 refCv = refCp-specGasConst;
2203 refEnergy = refPressure/(gamma-1.0)/refRho;
2207 refTemperature = 1.0;
2209 refSpecGasConst = 1.0;
2211 pressureScale = 1.0;
2212 temperatureScale = 1.0;
2213 velocityScale = 1.0;
2216 eosPtr->SetSpecificGasConstant(specGasConst);
2217 eosPtr->InitializeMaterialProperties();
2225 if(nonDimensionalMode != 0){
2226 refMu = refRho*refVelocity*refLength/refRe;
2227 refKappa = refCp*refMu/refPr;
2238 if(nonDimensionalMode == 1){
2239 if(transportModel ==
POWER){
2242 }
else if(nonDimensionalMode == 2){
2243 double gamma = eosPtr->GetGamma();
2244 if(transportModel ==
POWER){
2246 beta = pow(gamma-1,power);
2254 std::ostringstream refStateMsg;
2255 refStateMsg <<
"NavierStokesRHS Reference State: " << std::endl
2256 <<
"Reynolds Number: \t" << refRe << std::endl
2257 <<
"Prandtl Number: \t" << refPr << std::endl
2259 <<
"Length Reference: \t" << refLength << std::endl
2260 <<
"Density Reference: \t" << refRho << std::endl
2261 <<
"Pressure Reference: \t" << refPressure << std::endl
2262 <<
"Temperature Reference: \t" << refTemperature << std::endl
2264 <<
"Specific Gas Constant Reference: \t" << refSpecGasConst << std::endl
2265 <<
"Energy Reference: \t" << refEnergy << std::endl
2266 <<
"Velocity Reference: \t" << refVelocity << std::endl
2267 <<
"Speed of Sound Reference: \t" << refSndSpd << std::endl
2268 <<
"Cp Reference: \t" << refCp << std::endl
2269 <<
"Cv Reference: \t" << refCv << std::endl;
2271 if(nonDimensionalMode == 0) {
2272 refStateMsg <<
"NavierStokesRHS set to run in dimensional mode." << std::endl
2273 <<
"Specific Gas Constant: \t" << specGasConst << std::endl
2274 <<
"Speed of sound: \t" << refSndSpd << std::endl
2275 <<
"Cp: \t" << eosPtr->GetHeatCapacityCp() << std::endl
2276 <<
"Cv: \t" << eosPtr->GetHeatCapacityCv() << std::endl;
2280 globalPtr->StdOut(refStateMsg.str());
2282 std::cout << refStateMsg.str();
2285 std::ostringstream paramStatusMsg;
2291 int inputCFLHandle = paramState.GetDataIndex(
"inputCFL");
2292 if(inputCFLHandle >= 0){
2294 const double *inCFLPtr = cflData.
Data<
double>();
2295 inputCFL = *inCFLPtr;
2296 paramStatusMsg <<
"Input CFL = " << inputCFL << std::endl;
2300 paramStatusMsg <<
"Defaulting CFL = " << inputCFL << std::endl;
2305 int inputDTHandle = paramState.GetDataIndex(
"inputDT");
2306 if(inputDTHandle >= 0){
2308 const double *inputDTPtr = inputDTData.
Data<
double>();
2309 inputDT = *inputDTPtr;
2310 paramStatusMsg <<
"Input DT = " << inputDT << std::endl;
2314 paramStatusMsg <<
"User-specified DT is invalid, defaulting to CONSTANT_CFL mode." 2318 sigmaHandle = paramState.GetDataIndex(
"sigma");
2319 if(sigmaHandle >= 0){
2321 sigmaPtr = sigmaData.
Data<
double>();
2322 inputSigma = *sigmaPtr;
2328 int sigmaDilHandle = paramState.GetDataIndex(
"sigmaDilatation");
2329 if(sigmaDilHandle >= 0){
2331 sigmaDilPtr = sigmaData.
Data<
double>();
2332 inputSigmaDil = *sigmaDilPtr;
2335 inputSigmaDil = 0.0;
2339 if(inputSigma > 0.0 || inputSigmaDil > 0.0){
2340 artificialDissipation =
true;
2342 if(artificialDissipation){
2343 if(inputSigma < 0 || inputSigmaDil < 0){
2351 int paramHandle = paramState.GetDataIndex(
"Flags");
2352 if(paramHandle >= 0){
2354 const int *paramPtr = paramData.
Data<
int>();
2355 int optionFlags = *paramPtr;
2360 useWENO = optionFlags&
USEWENO;
2364 paramStatusMsg <<
"NavierStokesRHS: Viscous flow, enabling DV exchange." << std::endl;
2367 if(artificialDissipation){
2368 paramStatusMsg <<
"NavierStokesRHS: Artificial dissipation ON.";
2370 paramStatusMsg <<
" Enabling DV exchange.";
2371 paramStatusMsg << std::endl;
2374 if(useAlternateRHS){
2375 exchangeFlux =
true;
2376 paramStatusMsg <<
"NavierStokesRHS: Enabling flux exchange." << std::endl;
2384 globalPtr->StdOut(paramStatusMsg.str());
2385 globalPtr->FunctionExit(
"SetParamState");
2388 std::cout << paramStatusMsg.str();
2398 if(!gridInitialized)
2401 timeHandle = inState.GetDataIndex(
"simTime");
2405 timePtr = timeData.
Data<
double>();
2410 cflHandle = inState.GetDataIndex(
"simCFL");
2413 cflPtr = cflData.
Data<
double>();
2416 dtHandle = inState.GetDataIndex(
"simDT");
2419 dtPtr = dtData.
Data<
double>();
2422 rhoHandle = inState.GetDataIndex(
"rho");
2426 rhoPtr = rhoData.
Data<
double>();
2430 rhoVHandle = inState.GetDataIndex(
"rhoV");
2434 rhoVPtr = rhoVData.
Data<
double>();
2438 rhoEHandle = inState.GetDataIndex(
"rhoE");
2442 rhoEPtr = rhoEData.
Data<
double>();
2446 pressureHandle = inState.GetDataIndex(
"pressure");
2447 if(pressureHandle >= 0){
2449 pressurePtr = pressureData.
Data<
double>();
2452 temperatureHandle = inState.GetDataIndex(
"temperature");
2453 if(temperatureHandle >= 0){
2455 temperaturePtr = temperatureData.
Data<
double>();
2459 rhom1Handle = inState.GetDataIndex(
"rhom1");
2460 if(rhom1Handle >= 0){
2462 rhom1Ptr = rhom1Data.
Data<
double>();
2465 velocityHandle = inState.GetDataIndex(
"velocity");
2466 if(velocityHandle >= 0){
2468 velocityPtr = velocityData.
Data<
double>();
2471 internalEnergyHandle = inState.GetDataIndex(
"internalEnergy");
2472 if(internalEnergyHandle >= 0){
2474 internalEnergyPtr = internalEnergyData.
Data<
double>();
2477 scalarHandle = inState.GetDataIndex(
"scalarVars");
2478 if(scalarHandle >= 0){
2480 numScalar = inState.Meta()[scalarHandle].ncomp;
2481 scalarPtr = scalarData.Data<
double>();
2490 int stencilConnHandle = inState.GetDataIndex(
"sconn");
2491 if(stencilConnHandle >= 0){
2493 stencilConn = scData.
Data<
int>();
2500 int iMaskHandle = inState.GetDataIndex(
"imask");
2501 if(iMaskHandle >= 0){
2503 iMask = imData.
Data<
int>();
2510 int dilatationHandle = inState.GetDataIndex(
"divV");
2511 if(dilatationHandle >= 0){
2513 velDiv = dilData.
Data<
double>();
2521 muHandle = inState.GetDataIndex(
"mu");
2524 muPtr = muData.
Data<
double>();
2530 lambdaHandle = inState.GetDataIndex(
"lambda");
2531 if(lambdaHandle >= 0 ){
2533 lambdaPtr = lambdaData.
Data<
double>();
2539 kappaHandle = inState.GetDataIndex(
"kappa");
2540 if(kappaHandle >= 0){
2542 kappaPtr = kappaData.
Data<
double>();
2549 myStateBuffers.resize(numEquations+1,NULL);
2553 myStateBuffers[iEquation++] = rhoPtr;
2554 for(
int iDim = 0;iDim < numDim;iDim++)
2556 myStateBuffers[iEquation++] = rhoEPtr;
2557 for(
int iScalar = 0;iScalar <
numScalar;iScalar++)
2559 myRHSBuffers.resize(numEquations+1,NULL);
2561 if(dvBuffers.empty()){
2571 dvBuffers.resize(numDV);
2577 dvBuffers[0] = pressurePtr;
2580 pressurePtr = dvBuffers[0];
2583 dvBuffers[1] = temperaturePtr;
2586 temperaturePtr = dvBuffers[1];
2589 dvBuffers[2] = rhom1Ptr;
2592 rhom1Ptr = dvBuffers[2];
2596 for(
int iDim = 0;iDim < numDim;iDim++)
2598 if(internalEnergyPtr){
2599 dvBuffers[3+numDim] = internalEnergyPtr;
2602 internalEnergyPtr = dvBuffers[3+numDim];
2606 if(tvBuffers.empty()){
2609 tvBuffers.resize(numTV);
2613 tvBuffers[0] = muPtr;
2616 muPtr = tvBuffers[0];
2619 tvBuffers[1] = lambdaPtr;
2622 lambdaPtr = tvBuffers[1];
2625 tvBuffers[2] = kappaPtr;
2628 kappaPtr = tvBuffers[2];
2633 rhsState.AddField(
"vHat",
'n',1,8,
"");
2634 rhsState.AddField(
"viscidEnergy",
'n',3,8,
"");
2635 rhsState.AddField(
"inviscidFlux",
'n',numEquations,8,
"");
2636 rhsState.AddField(
"viscidFlux",
'n',numEquations,8,
"");
2637 rhsState.AddField(
"dFlux",
'n',numEquations,8,
"");
2638 if(refRe > 0 || artificialDissipation){
2639 rhsState.AddField(
"velGrad",
'n',numDim*numDim,8,
"");
2640 if(artificialDissipation){
2641 rhsState.AddField(
"velDiv",
'n',1,8,
"");
2642 rhsState.AddField(
"sigmaDiss",
'n',1,8,
"");
2643 rhsState.AddField(
"alphaDiss",
'n',1,8,
"");
2647 rhsState.AddField(
"tau",
'n',(numDim*(numDim+1))/2,8,
"");
2648 rhsState.AddField(
"temperatureGrad",
'n',numDim,8,
"");
2649 rhsState.AddField(
"heatFlux",
'n',numDim,8,
"");
2652 velHat = rhsState.template GetFieldBuffer<double>(
"vHat");
2653 viscidEnergy = rhsState.template GetFieldBuffer<double>(
"viscidEnergy");
2654 if(artificialDissipation){
2656 rhsState.SetFieldBuffer(
"velDiv",velDiv);
2658 velDiv = rhsState.template GetFieldBuffer<double>(
"velDiv");
2661 alphaDissipation = rhsState.template GetFieldBuffer<double>(
"alphaDiss");
2664 alphaDissipation = NULL;
2667 velocityHandle = inState.GetDataIndex(
"velocity");
2668 if(velocityHandle >= 0){
2670 velocityPtr = velocityData.
Data<
double>();
2674 if(inviscidFluxBuffers.empty()){
2675 inviscidFluxBuffers.resize(numEquations);
2677 int inviscidFluxHandle = rhsState.GetDataIndex(
"inviscidFlux");
2678 inviscidFluxBuffer = rhsState.GetRealFieldBuffer(inviscidFluxHandle);
2679 for(
int i=0;i<numEquations;i++){
2684 if(viscidFluxBuffers.empty()){
2685 viscidFluxBuffers.resize(numEquations);
2687 int viscidFluxHandle = rhsState.GetDataIndex(
"viscidFlux");
2688 viscidFluxBuffer = rhsState.GetRealFieldBuffer(viscidFluxHandle);
2689 for(
int i=0;i<numEquations;i++){
2694 if(dFluxBuffers.empty()){
2695 dFluxBuffers.resize(numEquations);
2697 int dFluxHandle = rhsState.GetDataIndex(
"dFlux");
2698 dFluxBuffer = rhsState.GetRealFieldBuffer(dFluxHandle);
2699 for(
int i=0;i<numEquations;i++){
2704 if(velGradBuffers.empty()){
2705 velGradBuffers.resize(numDim);
2707 int velGradHandle = rhsState.GetDataIndex(
"velGrad");
2708 velGradBufferPtr = rhsState.GetRealFieldBuffer(velGradHandle);
2709 for(
int i=0;i<numDim;i++){
2713 if(tauBuffers.empty()){
2714 tauBuffers.resize(numDim*numDim);
2716 int tauHandle = rhsState.GetDataIndex(
"tau");
2717 tauBufferPtr = rhsState.GetRealFieldBuffer(tauHandle);
2722 size_t bufferIndex = 0;
2723 for(
int i=0;i<numDim;i++){
2724 for(
int j=0;j<numDim;j++){
2725 int index = i*numDim + j;
2727 tauBuffers[index] = &tauBufferPtr[bufferIndex];
2730 tauBuffers[index] = tauBuffers[j*numDim+i];
2736 if(temperatureGradBuffers.empty()){
2737 temperatureGradBuffers.resize(numDim);
2739 int temperatureGradHandle = rhsState.GetDataIndex(
"temperatureGrad");
2740 temperatureGradBufferPtr = rhsState.GetRealFieldBuffer(temperatureGradHandle);
2741 for(
int i=0;i<numDim;i++){
2742 temperatureGradBuffers[i] = &temperatureGradBufferPtr[i*
numPointsBuffer];
2745 if(heatFluxBuffers.empty()){
2746 heatFluxBuffers.resize(numDim);
2748 int heatFluxHandle = rhsState.GetDataIndex(
"heatFlux");
2749 heatFluxBufferPtr = rhsState.GetRealFieldBuffer(heatFluxHandle);
2751 for(
int i=0;i<numDim;i++){
2760 eosPtr->SetupPressureBuffer(pressurePtr);
2761 eosPtr->SetupTemperatureBuffer(temperaturePtr);
2762 eosPtr->SetupSpecificVolumeBuffer(rhom1Ptr);
2763 eosPtr->SetupInternalEnergyBuffer(internalEnergyPtr);
2765 dScalarHandle = scalarHandle;
2766 stateInitialized =
true;
2768 SetTargetState(inState);
2776 if(!stateInitialized)
2779 targetState.Copy(inState);
2781 targetStateBuffers.resize(0);
2782 int numTargetVars = numEquations+1;
2783 targetStateBuffers.resize(numTargetVars,NULL);
2785 rhoTargetPtr = targetState.GetRealFieldBuffer(rhoHandle);
2786 rhoVTargetPtr = targetState.GetRealFieldBuffer(rhoVHandle);
2787 rhoETargetPtr = targetState.GetRealFieldBuffer(rhoEHandle);
2790 targetStateBuffers[iEquation++] = rhoTargetPtr;
2791 for(
int iDim = 0;iDim < numDim;iDim++)
2792 targetStateBuffers[iEquation++] = &rhoVTargetPtr[iDim*
numPointsBuffer];
2793 targetStateBuffers[iEquation++] = rhoETargetPtr;
2796 scalarTargetPtr = targetState.GetRealFieldBuffer(scalarHandle);
2797 for(
int iScalar = 0;iScalar <
numScalar;iScalar++)
2798 targetStateBuffers[iEquation++] = &scalarTargetPtr[iScalar*
numPointsBuffer];
2809 if(!stateInitialized)
2812 rhoPtr = inState.GetRealFieldBuffer(rhoHandle);
2813 rhoVPtr = inState.GetRealFieldBuffer(rhoVHandle);
2814 rhoEPtr = inState.GetRealFieldBuffer(rhoEHandle);
2817 myStateBuffers[iEquation++] = rhoPtr;
2818 for(
int iDim = 0;iDim < numDim;iDim++)
2820 myStateBuffers[iEquation++] = rhoEPtr;
2823 scalarPtr = inState.GetRealFieldBuffer(scalarHandle);
2824 for(
int iScalar = 0;iScalar <
numScalar;iScalar++)
2829 timePtr = inState.GetRealFieldBuffer(timeHandle);
2833 dtPtr = inState.GetRealFieldBuffer(dtHandle);
2836 cflPtr = inState.GetRealFieldBuffer(cflHandle);
2847 dRhoHandle = rhsState.GetDataIndex(
"rho");
2848 dRhoVHandle = rhsState.GetDataIndex(
"rhoV");
2849 dRhoEHandle = rhsState.GetDataIndex(
"rhoE");
2851 dScalarHandle = rhsState.GetDataIndex(
"scalarVars");
2854 rhoRHSPtr = rhsState.GetRealFieldBuffer(dRhoHandle);
2855 rhoVRHSPtr = rhsState.GetRealFieldBuffer(dRhoVHandle);
2856 rhoERHSPtr = rhsState.GetRealFieldBuffer(dRhoEHandle);
2859 scalarRHSPtr = rhsState.GetRealFieldBuffer(dScalarHandle);
2861 scalarRHSPtr = NULL;
2864 myRHSBuffers[iEquation++] = rhoRHSPtr;
2865 for(
int iDim = 0;iDim < numDim;iDim++)
2867 myRHSBuffers[iEquation++] = rhoERHSPtr;
2868 for(
int iScalar = 0;iScalar <
numScalar;iScalar++)
2883 if(!operatorInitialized)
2902 if(artificialDissipation){
2917 boundaryDepth,&periodicDirs[0],
2920 if(artificialDissipation){
2922 artDissOperator.boundaryDepth,&periodicDirs[0],
2926 artDissOperatorSplit.boundaryDepth,&periodicDirs[0],
2927 artDissConnSplit,
true);
2969 if(!operatorInitialized){
2970 myOperator.Copy(inOperator);
2977 boundaryDepth = myOperator.boundaryDepth;
2978 boundaryWeight = myOperator.boundaryWeight;
2979 numStencilValues = myOperator.numValues;
2982 for(
int iStencil = 0;iStencil <
numStencils;iStencil++)
2986 if(artificialDissipation){
2987 int interiorOrderAD = (myOperator.overallOrder-1)*2;
2989 int *startPtr = artDissOperator.stencilStarts;
2990 int numStencilsAD = artDissOperator.numStencils;
2992 for(
int iStencil = 0;iStencil < numStencilsAD;iStencil++)
2993 startPtr[iStencil]++;
2994 startPtr = artDissOperatorSplit.stencilStarts;
2995 numStencilsAD = artDissOperator.numStencils;
2996 for(
int iStencil = 0;iStencil < numStencilsAD;iStencil++)
2997 startPtr[iStencil]++;
3000 operatorInitialized =
true;
3004 if(SetupStencilConnectivities())
3007 gridPtr->SetupDifferentialOperator(&myOperator,&stencilConn[0]);
3008 std::ostringstream metricMessages;
3009 if(gridPtr->ComputeMetrics(metricMessages)){
3010 std::ostringstream failMessage;
3011 failMessage <<
"navierstokes::rhs::SetupOperators: Error: grid::ComputeMetrics failed " 3012 <<
"with messages:" << std::endl << metricMessages.str() << std::endl;
3014 globalPtr->ErrOut(failMessage.str());
3016 std::cout << failMessage.str();
3053 bool inParallel =
false;
3056 inParallel = omp_in_parallel();
3058 myThreadId = omp_get_thread_num();
3059 numThreads = omp_get_num_threads();
3065 fortranBufferIntervals.resize(numThreads);
3066 fortranPartitionBufferIntervals.resize(numThreads);
3073 std::vector<pcpp::IndexIntervalType> &threadPartitionBufferIntervals(gridPtr->ThreadPartitionBufferIntervals());
3074 if(threadPartitionBufferIntervals.size() != numThreads)
3077 threadPartitionBufferIntervals[myThreadId].Flatten(fortranPartitionBufferIntervals[myThreadId]);
3079 std::vector<pcpp::IndexIntervalType> &threadBufferIntervals(gridPtr->ThreadBufferIntervals());
3080 if(threadBufferIntervals.size() != numThreads)
3082 threadBufferIntervals[myThreadId].Flatten(fortranBufferIntervals[myThreadId]);
3085 std::vector<size_t>::iterator opIt = fortranPartitionBufferIntervals[myThreadId].begin();
3086 while(opIt != fortranPartitionBufferIntervals[myThreadId].end()){
3091 std::vector<size_t>::iterator bufIt = fortranBufferIntervals[myThreadId].begin();
3092 while(bufIt != fortranBufferIntervals[myThreadId].end()){
3109 numThreads = omp_get_num_threads();
3112 if(gridPtr != &inGrid){
3117 myThreadId = omp_get_thread_num();
3127 partitionInterval = gridPtr->PartitionInterval();
3128 partitionBufferInterval = gridPtr->PartitionBufferInterval();
3129 numPointsPartition = partitionInterval.NNodes();
3132 gridCommPtr = &gridPtr->Communicator();
3133 subRegionsPtr = &gridPtr->SubRegions();
3134 gridCommPtr->CartNeighbors(gridNeighbors);
3145 gridInitialized =
true;
3146 gridDX.resize(numDim,1.0);
3147 gridDXM1.resize(numDim,1.0);
3150 std::vector<double> &dX(inGrid.DX());
3151 std::vector<double>::iterator dXIt = dX.begin();
3152 std::vector<double>::iterator gridDXIt = gridDX.begin();
3153 std::vector<double>::iterator gridDXM1It = gridDXM1.begin();
3154 while(dXIt != dX.end()){
3155 double dx = *dXIt++;
3157 *gridDXM1It++ = 1/dx;
3163 while(dXIt != dX.end()){
3164 double dx = *dXIt++;
3171 std::vector<size_t>::iterator opIt =
opInterval.begin();
3178 InitThreadIntervals();
3195 if(periodicDirs.empty()){
3196 periodicDirs.resize(numDim,1);
3222 void SetDomainBCs(std::vector<BoundaryType> &domainBoundaries,std::vector<BCType> &domainBCs)
3225 domainBCsPtr = &domainBCs;
3226 domainBoundariesPtr = &domainBoundaries;
3227 assert(subRegionsPtr != NULL);
3229 std::vector<GridRegionType> &gridRegions(*subRegionsPtr);
3231 int numBCs = domainBCs.size();
3232 int numGridRegions = gridRegions.size();
3233 int numDomainBoundaries = domainBoundaries.size();
3234 int holeStencil = myOperator.numStencils;
3235 int adHole = artDissOperator.numStencils;
3236 int adSplitHole = artDissOperatorSplit.numStencils;
3237 int adBoundaryDepth = artDissOperator.boundaryDepth;
3238 int adSplitBoundaryDepth = artDissOperatorSplit.boundaryDepth;
3249 if(numDomainBoundaries == 0)
3254 std::vector<size_t> partitionBufferExtent;
3255 partitionBufferInterval.Flatten(partitionBufferExtent);
3257 typename std::vector<BoundaryType>::iterator dbIt = domainBoundaries.begin();
3264 dbIt = domainBoundaries.begin();
3265 while(dbIt != domainBoundaries.end()){
3266 BoundaryType &domainBoundary(*dbIt++);
3267 GridRegionType &gridRegion(gridRegions[domainBoundary.
GridRegionID()]);
3269 int boundaryNormalDirection = gridRegion.normalDirection;
3273 BCType &boundaryCondition(domainBoundary.
BC());
3275 int bcType = boundaryCondition.BCType();
3281 partitionBufferInterval,
3283 if(opInterval.
NNodes() > 0) {
3284 std::vector<size_t> opExtents;
3285 opInterval.
Flatten(opExtents);
3295 haloPtr->PackMessageBuffers(maskMessageId,0,1,iMask);
3296 haloPtr->SendMessage(maskMessageId,gridNeighbors,*gridCommPtr);
3297 haloPtr->ReceiveMessage(maskMessageId,gridNeighbors,*gridCommPtr);
3298 haloPtr->UnPackMessageBuffers(maskMessageId,0,1,iMask);
3302 dbIt = domainBoundaries.begin();
3303 while(dbIt != domainBoundaries.end()){
3304 BoundaryType &domainBoundary(*dbIt++);
3305 GridRegionType &gridRegion(gridRegions[domainBoundary.
GridRegionID()]);
3307 int boundaryNormalDirection = gridRegion.normalDirection;
3311 BCType &boundaryCondition(domainBoundary.
BC());
3312 int bcType = boundaryCondition.BCType();
3316 boundaryDepth,
holeMask,iMask,stencilConn);
3319 if(artificialDissipation){
3321 adBoundaryDepth,
holeMask,iMask,artDissConn);
3323 adSplitBoundaryDepth,
holeMask,iMask,artDissConnSplit);
3336 virtual int ProcessHalo(HaloType &inHalo,GridType &inGrid,StateType &inState,StateType &rhsState)
3345 globalPtr = &inGlobal;
3349 {
return(dvBuffers); };
3351 {
return(tvBuffers); };
3353 {
return(myStateBuffers); };
3355 {
return(myRHSBuffers); };
3357 {
return(inviscidFluxBuffer); };
3359 {
return(dFluxBuffer); };
3361 {
return(numEquations); };
3365 {
return(coeffsWENO); };
3367 {
return(velGradBuffers); };
3369 {
return(temperatureGradBuffers); };
3375 minSpacing.resize(numDim,std::numeric_limits<double>::max());
3378 for(
int iDim = 0;iDim < numDim;iDim++){
3379 minSpacing[iDim] = gridDX[iDim];
3407 if(minSpacing.empty())
3409 double spaceFactor = 0.0;
3410 for(
int iDim = 0;iDim < numDim;iDim++)
3411 spaceFactor += 1.0/(minSpacing[iDim]);
3417 if(*dtPtr != inputDT){
3419 *cflPtr = refSndSpd*(*dtPtr)*spaceFactor;
3430 if(*cflPtr != inputCFL){
3432 *dtPtr = *cflPtr/(refSndSpd*spaceFactor);
3449 myThreadId = omp_get_thread_num();
3450 numThreads = omp_get_num_threads();
3470 temperatureHandle = -1;
3471 pressureHandle = -1;
3489 velocityHandle = -1;
3491 internalEnergyHandle = -1;
3494 numPointsPartition = 0;
3509 rhoTargetPtr = NULL;
3510 rhoVTargetPtr = NULL;
3511 rhoETargetPtr = NULL;
3512 scalarTargetPtr = NULL;
3514 temperaturePtr = NULL;
3518 scalarRHSPtr = NULL;
3520 momentumFluxPtr = NULL;
3521 energyFluxPtr = NULL;
3522 scalarFluxPtr = NULL;
3523 fluxComponentPtr = NULL;
3526 internalEnergyPtr = NULL;
3527 massSourcePtr = NULL;
3528 momentumSourcePtr = NULL;
3529 energySourcePtr = NULL;
3530 scalarSourcePtr = NULL;
3532 viscidFluxBuffer = NULL;
3533 inviscidFluxBuffer = NULL;
3535 scaledPressure = NULL;
3538 viscidEnergy = NULL;
3539 velGradBufferPtr = NULL;
3542 alphaDissipation = NULL;
3547 temperatureGradBufferPtr = NULL;
3553 subRegionsPtr = NULL;
3554 domainBoundariesPtr = NULL;
3555 domainBCsPtr = NULL;
3560 artDissConnSplit = NULL;
3565 stateInitialized =
false;
3566 gridInitialized =
false;
3567 threadsInitialized =
false;
3568 operatorInitialized =
false;
double * temperatureGradBufferPtr
void AccumulateViscousFluxes(double a, std::vector< double *> &X, std::vector< double *> &Y, int threadId, bool fullInterval=false)
double * heatFluxBufferPtr
OperatorType artDissOperatorSplit
void const size_t const size_t const size_t const int const int const double * gridMetric
size_t * numPointsStencil
bool artificialDissipation
std::vector< double * > myStateBuffers
subroutine vhatcomponent(numDim, numPointsBuffer, bufferSizes, opInterval, gridType, gridMetric, velDir, velocity, velHatComponent)
int ComputeDVBuffer(const pcpp::IndexIntervalType ®ionInterval, const std::vector< size_t > &bufferSizes, const std::vector< double *> &stateBuffers, std::vector< double *> &dvBuffers)
subroutine assignmentyx(numDim, numPoints, bufferSize, bufferInterval, X, Y)
ASSIGNMENTYX point-wise operator performing Y = X.
static const int holeMask
int ComputeTauBuffer(const pcpp::IndexIntervalType ®ionInterval, const std::vector< size_t > &bufferSize, int gridType, const std::vector< double > &gridMetrics, const std::vector< double > &gridJacobian, const double &Re, const std::vector< double *> &tvBuffers, const std::vector< double *> &velGradBuffers, std::vector< double *> &tauBuffers)
subroutine noslip_isothermal(numDim, bufferSizes, numPointsBuffer, patchNormalDir, patchSizes, numPointsPatch, numPatchPointsOp, patchPointsOp, gridType, gridMetric, jacobianDeterminant, bcParams, gasParams, rhoBuffer, rhoVBuffer, rhoEBuffer, numscalar, scalarBuffer, rhoRHS, rhoVRHS, rhoERHS, scalarRHS, rhoTarget, rhoVTarget, rhoETarget, scalarTarget, muBuffer, lambdaBuffer)
subroutine dissipationweight(numDim, bufferSize, numPoints, bufferInterval, sigmaDissipation, sigmaDilatation, dilatationRamp, dilatationCutoff, divV, sigmaDiss)
subroutine ywxpy(numDim, numPoints, bufferSize, bufferInterval, W, X, Y)
YWXPY point-wise operator performing Y = WX + Y, where all are vectors.
void SetDomainBCs(std::vector< BoundaryType > &domainBoundaries, std::vector< BCType > &domainBCs)
std::vector< size_t > opInterval
rhs_base_t::BCType BCType
virtual int ProcessHalo(HaloType &inHalo, GridType &inGrid, StateType &inState, StateType &rhsState)
std::vector< std::vector< size_t > > fortranPartitionBufferIntervals
double * velGradBufferPtr
void PackBufferMessage(int messageId, int numEquations, double *dataBuffer, int threadId)
std::vector< double * > & TemperatureGradBuffers()
const std::string & Name() const
std::vector< BoundaryType > * domainBoundariesPtr
int HaloSetup(HaloType &inHalo)
int HandleBoundaryConditions(int threadId)
void size_t int size_t int size_t int int * stencilSizes
WENO::CoeffsWENO & WENOCoefficients()
int ExchangeDV(int threadId)
int SetGrid(GridType &inGrid)
void GetFlatIndices(const sizeextent &extent, ContainerType &indices) const
int InitStep(GridType &inGrid, StateType &inState)
std::vector< int > gridNeighbors
std::vector< int > & PeriodicDirs()
void Flatten(ContainerType &output) const
double * internalEnergyPtr
void size_t int size_t int size_t int int int int double int int double double *void size_t int size_t int int int int int double int size_t size_t size_t double double *void size_t int size_t int size_t size_t int double int double double *void size_t size_t size_t double double * X
int ComputeStressTensor(int threadId)
OperatorType artDissOperator
std::vector< double > gridDXM1
std::vector< double * > heatFluxBuffers
int ComputeDV(int threadId)
std::vector< double * > dvBuffers
void double double double double double double double *void const size_t const size_t const size_t const double const double const double const double const double double * sigmaDiss
rhs_base_t::BoundaryType BoundaryType
double * viscidFluxBuffer
void size_t int size_t int size_t int int int int double int int double * uBuffer
int InitThreadIntervals()
eos::perfect_gas * perfectGasPtr
int ComputeTV(int threadId)
int ExchangeInviscidFlux(int threadId)
subroutine applyoperator(numDim, dimSizes, numComponents, numPoints, opDir, opInterval, numStencils, stencilSizes, stencilStarts, numValues, stencilWeights, stencilOffsets, stencilID, U, dU)
applyoperator applies an operator specified as a stencil set to the provided state data ...
std::vector< BCType > * domainBCsPtr
void AccumulateToRHS(double a, double *dfBuffer, int threadId)
void UnpackDVMessage(int threadId)
void ReconstPointVal(double vals[], CoeffsWENO &coeffs, int group, double &pVal)
int ComputeHeatFluxBuffer(const pcpp::IndexIntervalType ®ionInterval, const std::vector< size_t > &bufferSize, int gridType, const std::vector< double > &gridMetrics, const std::vector< double > &gridJacobian, const double &Re, const std::vector< double *> &tvBuffers, const std::vector< double *> &temperatureGradBuffers, std::vector< double *> &heatFluxBuffers)
int ExchangeBuffer(int messageId, int numEquations, double *dataBuffer, int threadId)
void size_t int size_t int size_t int int int * stencilStarts
void const size_t const size_t * gridSizes
double EntropyFixEta(int n, double lambdaL[], double lambdaR[])
void RenewStream(std::ostringstream &outStream)
std::vector< double > gridMetric
WENO::CoeffsWENO coeffsWENO
virtual int ProcessInterior(GridType &inGrid, StateType &inState, StateType &rhsState)
void const size_t const size_t const size_t const int const int const double const double const double const double const double * velHat
std::vector< double * > & TVBuffers()
int SetState(StateType &inState)
void const size_t const size_t const size_t const int const int * gridType
std::vector< size_t > bufferSizes
int ApplyDissipationOperator(int dDir, int numComponents, double *uBuffer, double *duBuffer, bool split=false, int threadId=0)
void SpongeZone(const std::vector< size_t > &bufferSizes, int normalDir, const pcpp::IndexIntervalType &spongeRegion, const pcpp::IndexIntervalType &partitionOpInterval, const pcpp::IndexIntervalType &opBufferInterval, const double *bcParams, const int *bcFlags, const int *iMask, const int disabledMask, const int numEquations, const std::vector< double *> &stateBuffers, const std::vector< double *> &targetBuffers, std::vector< double *> &rhsBuffers)
void const size_t const size_t const size_t * opInterval
subroutine yxy(numDim, numPoints, bufferSize, bufferInterval, X, Y)
YXY point-wise operator performing Y = XY (all vectors)
subroutine scalarflux1d(numDim, numPoints, gridSizes, opInterval, numScalars, scalarBuffer, velHat, fluxBuffer)
Flux for scalar transport.
std::vector< double * > targetStateBuffers
int PrepareBuffers(int threadId)
static const int EXCHANGEFLUX
std::vector< double > minSpacing
void const size_t const size_t const size_t const int const double const int const double * velocity
simulation::grid::subregion GridRegionType
subroutine assignmentxa(numDim, numPoints, bufferSize, bufferInterval, a, X)
ASSIGNMENTXA point-wise operator performing X = scalar a.
virtual int ProcessBoundary(GridType &inGrid, StateType &inState, StateType &rhsState)
std::vector< double * > & StateBuffers()
std::vector< double * > viscidFluxBuffers
subroutine strongflux1d(numDim, fluxDir, gridSizes, numPoints, opInterval, gridType, gridMetric, tauBuffer, energyBuffer, fluxBuffer)
Compute the curvilinear cartesian viscous fluxes in 1 dimension.
void UnpackBufferMessage(int messageId, int numEquations, double *dataBuffer, int threadId)
std::vector< double * > & VelGradBuffers()
int InviscidFlux(int iDim, int threadId)
int DetectHoles(int numDim, size_t *dimSizes, size_t *opInterval, int boundaryDepth, int holeBit, int *inMask, int *stencilID)
Detect unstructured holes and set up boundary stencil id's accordingly.
fixtures::CommunicatorType * gridCommPtr
double * InviscidFluxBuffer()
int SetupRHS(GridType &inGrid, StateType &inState, StateType &rhsState, int myThreadId=0)
virtual int ProcessInterior2(GridType &inGrid, StateType &inState, StateType &rhsState)
void ExchangeBufferMessage(int messageId)
void GridScaleRHS(int threadId)
void size_t int size_t int size_t int int int int double int * stencilOffsets
std::vector< int > periodicDirs
int HandleSources(const pcpp::IndexIntervalType ®ionInterval)
std::vector< size_t > gridSizes
subroutine ijkgradtoxyzdiv(numDim, numPoints, gridSizes, opInterval, gridType, gridJacobian, gridMetric, gradVBuffer, divBuffer)
int CreateMessageBuffers(int numComponents, int componentSize=8)
subroutine farfield(numDim, bufferSizes, numPointsBuffer, patchNormalDir, patchSizes, numPointsPatch, numPatchPointsOp, patchPointsOp, gridType, gridMetric, jacobianDeterminant, bcParams, gasParams, rhoBuffer, rhoVBuffer, rhoEBuffer, viscousFluxBuffer, numscalar, scalarBuffer, rhoRHS, rhoVRHS, rhoERHS, scalarRHS, rhoTarget, rhoVTarget, rhoETarget, scalarTarget)
int ComputeVelGrad(int threadId)
void const size_t const size_t const int const size_t const size_t const size_t const size_t const int const double const double const double const double const double const double const double const double const int * numScalar
std::vector< double > gridDX
int ArtificialDissipation(int threadId)
const double * sigmaDilPtr
std::vector< double * > velGradBuffers
subroutine flux1d(numDim, numPoints, gridSizes, opInterval, fluxDir, gridType, gridMetric, rhoBuffer, rhoVBuffer, rhoEBuffer, velHat, pressureBuffer, fluxBuffer)
int Initialize(GridType &inGrid, StateType &inState, StateType &inParam, OperatorType &inOperator)
static const int EXCHANGEDV
std::vector< std::vector< double * > > haloFluxBuffers
int Initialize(stencilset &stencilSet1, stencilset &stencilSet2, int interiorOrder)
std::vector< std::vector< double * > > haloDVBuffers
std::vector< double * > & RHSBuffers()
int HandleSources(int myThreadId)
int SetupStencilConnectivities()
int ComputeTVBufferPower(const pcpp::IndexIntervalType ®ionInterval, const std::vector< size_t > &bufferSize, const double *temperatureBuffer, std::vector< double *> &tvBuffer, const double beta, const double power, const double bulkViscFac, const double specificHeat, const double prandtlNumber)
Compute transport coefficients using the power law.
void Project(int n, int m, double T[], double vals[], double res[])
double * momentumSourcePtr
subroutine slip_adiabatic(numDim, bufferSizes, numPointsBuffer, patchNormalDir, patchSizes, numPointsPatch, numPatchPointsOp, patchPointsOp, gridType, gridMetric, jacobianDeterminant, bcParams, gasParams, rhoBuffer, rhoVBuffer, rhoEBuffer, numscalar, scalarBuffer, rhoRHS, rhoVRHS, rhoERHS, scalarRHS, rhoTarget, rhoVTarget, rhoETarget, scalarTarget)
Main encapsulation of MPI.
void LinearSponge(const std::vector< size_t > &bufferSizes, int normalDir, const pcpp::IndexIntervalType &spongeRegion, const pcpp::IndexIntervalType &partitionOpInterval, const pcpp::IndexIntervalType &opBufferInterval, const double *bcParams, const int *bcFlags, const int numEquations, const std::vector< double *> &stateBuffers, const std::vector< double *> &targetBuffers, std::vector< double *> &rhsBuffers)
base API for equations of state
int CreateStencilConnectivity(int numDim, size_t *dimSizes, size_t *opInterval, int boundaryDepth, int *stencilID, bool fortranInterval=false)
Creates simple stencil connectivity assuming all domain boundaries are physical boundaries.
void size_t int size_t int size_t int int int int double int int double double *void size_t int size_t int int int int int double int size_t size_t size_t double double *void size_t int size_t int size_t size_t int double int double double *void size_t size_t size_t double * a
size_t numPointsPartition
simulation::grid::halo HaloType
sizeextent RelativeTranslation(const sizeextent &inOrigin, const sizeextent &inDestination) const
int ComputeViscidEnergyFlux(int threadId)
double * alphaDissipation
static const int ALTERNATERHS
void SetGlobal(pcpp::ParallelGlobalType &inGlobal)
subroutine alphaweight(numDim, numPointsBuffer, bufferSizes, opInterval, gridType, gridMetric, alphaDir, alphaW)
int ExchangeState(int threadId)
int ComputeVelDiv(int threadId)
void size_t int size_t int size_t int int int int double int int double double * duBuffer
void const size_t const size_t const size_t const int const double * gridJacobian
int ComputeHeatFlux(int threadId)
pcpp::IndexIntervalType partitionBufferInterval
const double * bulkViscFacPtr
void const size_t const size_t * bufferSizes
int ComputeTemperatureGrad(int threadId)
void size_t int size_t int size_t int int int int double int int double double *void size_t int size_t int int int int int double int size_t size_t size_t double double *void size_t int size_t int size_t size_t int double int double double *void size_t size_t size_t double double double * Y
void size_t int size_t int size_t int int int int double int int double double *void size_t int size_t int int int int int double int size_t size_t size_t double double *void size_t int size_t int size_t size_t int double int double double *void size_t size_t size_t * bufferInterval
Perfect Gas Equation of State.
double TimeStep(double *outCFL=NULL)
void const size_t const size_t const int const size_t const size_t const size_t const size_t const int const double const double const double const double * gasParams
std::vector< std::vector< size_t > > fortranBufferIntervals
simulation::rhs::domain_base< GridType, StateType, OperatorType > rhs_base_t
subroutine zxy(numDim, numPoints, bufferSize, bufferInterval, X, Y, Z)
ZXY point-wise operator performing Z = XY (all vectors)
double * fluxComponentPtr
double * inviscidFluxBuffer
void double double double double double double double *void const size_t const size_t const size_t const double * sigmaDissipation
std::vector< double * > inviscidFluxBuffers
void size_t int size_t int size_t int int int int double * stencilWeights
void size_t int size_t int size_t int * numStencils
int SetRHSState(StateType &rhsState)
int SetHalo(HaloType &inHalo)
int SetParamState(const StateType ¶mState)
std::vector< size_t > Sizes() const
std::vector< GridRegionType > * subRegionsPtr
std::vector< double * > myRHSBuffers
subroutine yaxpy(N1, N2, N3, localInterval, a, X, Y)
std::vector< double * > tvBuffers
int SetupOperators(OperatorType &inOperator)
Simple Block Structured Mesh object.
double * sigmaDissipation
std::vector< bool > haveHalo
int ApplyOperator(int dDir, int numComponents, double *uBuffer, double *duBuffer, int threadId)
void size_t int * numComponents
int SetTargetState(StateType &inState)
std::vector< double * > dFluxBuffers
int ViscidFlux(int iDim, int threadId)
std::vector< bool > dvBuffersMine
void SetMask(const std::vector< size_t > &opIndices, int maskBits, int *inMask)
std::vector< double * > temperatureGradBuffers
subroutine sat_form_roe_matrices(ND, u_in, u_out, gamma, norm, tmat, tinv, lambda)
std::vector< std::vector< double * > > haloStateBuffers
pcpp::IndexIntervalType partitionInterval
void InitSimple(const ContainerType &inSize)
void const size_t * numPointsBuffer
void const size_t const size_t const int const size_t const size_t const size_t const size_t const int const double const double const double * bcParams
void FC_MODULE(euler, flux1d, EULER, FLUX1D)(const int *numDim
int ApplyWENO(int iDim, int threadId)
void PackDVMessage(int threadId)
int InitializeState(StateType &inState)
std::vector< double * > & DVBuffers()
std::vector< double * > tauBuffers