PlasCom2  1.0
XPACC Multi-physics simluation application
TestWENOParallel.C
Go to the documentation of this file.
1 // Parallel unit tests for WENO routines and classes
2 
3 #include "Testing.H"
4 #include "Simulation.H"
5 #include "RKTestFixtures.H"
6 #include "EulerRHS.H"
7 #include "Stencil.H"
8 #include "PCPPCommUtil.H"
9 #include "PCPPReport.H"
10 #include "PCPPIntervalUtils.H"
11 #include "TestFixtures.H"
12 #include "EulerTestFixtures.H"
13 #include "SATKernels.H"
14 #include "WENO.H"
15 
16 
17 void TestWENO_ApplyWENO(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
18 {
19  // Test ApplyWENO routine
20 
21  double eps = 1e-12;
22  double epsLarge = 1e-3; // right near jumps, WENO reports values that are not quite zero
23  int stenWidth = 3; // max length of one WENO stencil
24 
25  typedef plascom2::operators::sbp::base operator_t;
26  typedef simulation::grid::halo halo_t;
27  typedef simulation::state::base state_t;
28  typedef pcpp::ParallelGlobalType global_t;
29  typedef pcpp::CommunicatorType comm_t;
30  typedef pcpp::IndexIntervalType interval_t;
31 
32  typedef simulation::grid::parallel_blockstructured grid_t;
33  typedef simulation::domain::base<grid_t,state_t> domain_t;
34  typedef euler::rhs<grid_t,state_t,operator_t> rhs_t;
35  typedef WENO::CoeffsWENO weno_t;
36 
37  std::ostringstream messageStream;
38 
39  // --- Equal ---
40  // Tests that when states are initially equal, result dFluxBuffer should be 0
41  for(int numDim = 2;numDim < 4;numDim++){
42 
43  grid_t simGrid;
44  operator_t sbpOperator;
45  halo_t &simHalo(simGrid.Halo());
46  state_t simState;
47  state_t paramState;
48  rhs_t rhsWENO;
49  weno_t &coeffsWENO(rhsWENO.WENOCoefficients());
50 
51  bool parTestWENO = true;
52  std::vector<size_t> gridSizes(numDim,21);
53 
54  if(testfixtures::CreateSimulationFixtures(simGrid,simHalo,sbpOperator,coeffsWENO,
55  gridSizes,2,testComm,&messageStream)){
56  std::cout << "CreateSimulationFixtures failed: " << messageStream.str() << std::endl;
57  parTestWENO = false;
58  }
59 
60  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Equal:InitFixtures",parTestWENO);
61 
62  if(!parTestWENO)
63  return;
64 
65  pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());
66  pcpp::IndexIntervalType &partitionInterval(simGrid.PartitionInterval());
67  std::vector<size_t> &bufferSizes(simGrid.BufferSizes());
68 
69  if(testfixtures::euler::SetupEulerState(simGrid,simState,paramState,0)){
70  parTestWENO = false;
71  }
72 
73  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Equal:SetupState",parTestWENO);
74 
75  if(!parTestWENO)
76  return;
77 
78  std::cout << messageStream.str() << std::endl;
79  std::cout << "WENO:ApplyWENO:Equal:Testing WENO in " << numDim << " dimensions." << std::endl;
80 
81  size_t numPointsBuffer = simGrid.BufferSize();
82 
83  double t = 0.0;
84  double gamma = 1.4;
85  double dt = .001;
86  double cfl = 1.0;
87  double Re = 0.0;
88  std::vector<double> numbers(4,0.0);
89  int flags = USEWENO;
90 
91  // Add some fields that are not usual
92  paramState.Create(numPointsBuffer,0);
93 
94  paramState.SetFieldBuffer("gamma",&gamma);
95  paramState.SetFieldBuffer("inputDT",&dt);
96  paramState.SetFieldBuffer("inputCFL",&cfl);
97  paramState.SetFieldBuffer("refRe",&Re);
98  paramState.SetFieldBuffer("Flags",&flags);
99  paramState.SetFieldBuffer("Numbers",&numbers[0]);
100 
101  std::vector<double> rho(numPointsBuffer,-1000.0);
102  std::vector<double> rhoV(numDim*numPointsBuffer,-1000.0);
103  std::vector<double> rhoE(numPointsBuffer,-1000.0);
104 
105  size_t iStart = partitionBufferInterval[0].first;
106  size_t iEnd = partitionBufferInterval[0].second;
107  size_t jStart = partitionBufferInterval[1].first;
108  size_t jEnd = partitionBufferInterval[1].second;
109 
110  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
111  // computed correctly
112  size_t kStart = 0;
113  size_t kEnd = 0;
114  if (numDim == 3) {
115  kStart = partitionBufferInterval[2].first;
116  kEnd = partitionBufferInterval[2].second;
117  }
118 
119  for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
120  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
121  for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
122  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
123  for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
124  size_t bufferIndex = jBufferIndex + iIndex;
125  rho[bufferIndex] = 1.0;
126  for (int vDim=0; vDim<numDim; vDim++) {
127  rhoV[bufferIndex + vDim*numPointsBuffer] = 1.0;
128  }
129  rhoE[bufferIndex] = 4.0 + 0.5*numDim;
130  }
131  }
132  }
133 
134 
135  std::vector<double> dRho(numPointsBuffer,-1000.0);
136  std::vector<double> dRhoV(numDim*numPointsBuffer,-1000.0);
137  std::vector<double> dRhoE(numPointsBuffer,-1000.0);
138 
139  std::vector<double> T(numPointsBuffer,0.0);
140  std::vector<double> p(numPointsBuffer,0.0);
141  std::vector<double> V(numDim*numPointsBuffer,0.0);
142 
143  std::vector<double> rhom1(numPointsBuffer,0.0);
144 
145 
146  simState.SetFieldBuffer("simTime",&t);
147  simState.SetFieldBuffer("rho",&rho[0]);
148  simState.SetFieldBuffer("rhoV",&rhoV[0]);
149  simState.SetFieldBuffer("rhoE",&rhoE[0]);
150 
151  simState.SetFieldBuffer("pressure",&p[0]);
152  simState.SetFieldBuffer("temperature",&T[0]);
153  simState.SetFieldBuffer("velocity",&V[0]);
154  simState.SetFieldBuffer("rhom1",&rhom1[0]);
155 
156  state_t rhsState(simState);
157 
158  rhsState.SetFieldBuffer("rho",&dRho[0]);
159  rhsState.SetFieldBuffer("rhoV",&dRhoV[0]);
160  rhsState.SetFieldBuffer("rhoE",&dRhoE[0]);
161 
162  // State/Soln initialization
163  std::vector<double> inParams;
164  std::vector<int> inFlags;
165 
166  // if(InitializeAcousticPulse(simGrid,simState,
167  // paramState,inParams,
168  // inFlags,0)){
169  // // then it failed.
170  // }
171 
172  std::cout << "Calling rhsWENO initialize." << std::endl;
173  if(rhsWENO.Initialize(simGrid,simState,paramState,sbpOperator)){
174  parTestWENO = false;
175  }
176  std::cout << "Calling rhsWENO initialize DONE!!." << std::endl;
177  std::vector<int> numThreadsDim;
178  std::cout << "MADE IT HERE" << std::endl;
179  if(simGrid.SetupThreads(numThreadsDim)){
180  std::cout << "Setting up threads for grid object failed." << std::endl;
181  return;
182  }
183  std::cout << "MADE IT HERE2" << std::endl;
184  rhsWENO.InitThreadIntervals();
185  std::cout << "MADE IT HERE3" << std::endl;
186 
187  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Equal:InitRHS",parTestWENO);
188  if(!parTestWENO)
189  return;
190 
191  if(rhsWENO.SetupRHS(simGrid,simState,rhsState)){
192  parTestWENO = false;
193  }
194 
195  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Equal:SetupRHS",parTestWENO);
196 
197  if(!parTestWENO)
198  return;
199 
200  //std::cout << "RHO BEFORE: " << std::endl;
201  //pcpp::io::DumpContents(std::cout,rho," ");
202  //std::cout << std::endl;
203 
204  // Exchange the state
205  rhsWENO.ExchangeState(0);
206 
207  //std::cout << "RHO AFTER: " << std::endl;
208  //pcpp::io::DumpContents(std::cout,rho," ");
209  //std::cout << std::endl;
210 
211  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
212  // computed correctly
213  size_t kFinal = 1;
214  if (numDim == 3) {
215  kFinal = bufferSizes[2];
216  }
217 
218  for (size_t kIndex = 0; kIndex < kFinal && parTestWENO; kIndex++) {
219  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
220  for (size_t jIndex = 0; jIndex < bufferSizes[1] && parTestWENO; jIndex++) {
221  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
222  for (size_t iIndex = 0; iIndex < bufferSizes[0] && parTestWENO; iIndex++) {
223  size_t bufferIndex = jBufferIndex + iIndex;
224 
225  int c = 0;
226  if (iIndex < iStart || iIndex > iEnd) c++;
227  if (jIndex < jStart || jIndex > jEnd) c++;
228  if (numDim == 3 && (kIndex < kStart || kIndex > kEnd)) c++;
229  if (c != 1) continue; // then we are not in a halo region
230 
231  if (std::abs(rho[bufferIndex] - 1.0) > eps) {
232  parTestWENO = false;
233  std::cout << "Failed state exchange, rho " << rho[bufferIndex] << std::endl;
234  }
235  for (int vDim=0; vDim<numDim; vDim++) {
236  if (std::abs(rhoV[bufferIndex + vDim*numPointsBuffer] - 1.0) > eps) {
237  parTestWENO = false;
238  std::cout << "Failed state exchange, rhoV" << vDim << " "
239  << rhoV[bufferIndex + vDim*numPointsBuffer] << std::endl;
240  }
241  }
242  if (std::abs(rhoE[bufferIndex] - (4.0 + 0.5*numDim)) > eps) {
243  parTestWENO = false;
244  std::cout << "Failed state exchange, rhoE " << rhoE[bufferIndex] << std::endl;
245  }
246 
247  if (!parTestWENO) {
248  std::cout << "iIndex = " << iIndex << ", jIndex = " << jIndex
249  << ", kIndex = " << kIndex << std::endl;
250  std::cout << "iStart = " << iStart << ", jStart = " << jStart
251  << ", kStart = " << kStart << std::endl;
252  std::cout << "iEnd = " << iEnd << ", jEnd = " << jEnd
253  << ", kEnd = " << kEnd << std::endl;
254  }
255  }
256  }
257  }
258 
259  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Equal:ExchangeState",parTestWENO);
260 
261  if(!parTestWENO)
262  return;
263 
264  double *inviscidFluxBuffer(rhsWENO.InviscidFluxBuffer());
265  double *dFluxBuffer(rhsWENO.DFluxBuffer());
266 
267  for(int iDim = 0;iDim < numDim;iDim++){
268 
269  // Running WENO in dimension (iDim)
270  std::cout << "Checking WENO in dimension(" << iDim << ")"
271  << std::endl;
272 
273  // initialize dFluxBuffer so test can't pass because data was initialized to zero
274  for (size_t i=0; i<numPointsBuffer; i++) {
275  for (int j=0; j<numDim+2; j++) {
276  size_t bufferIndex = j*numPointsBuffer + i;
277  dFluxBuffer[bufferIndex] = -1000;
278  }
279  }
280 
281  rhsWENO.ComputeDV(0);
282  rhsWENO.InviscidFlux(iDim,0);
283  rhsWENO.ExchangeInviscidFlux(0);
284 
285  //if(numDim == 2){
286  // std::cout << "Velocity = " << std::endl;
287  // for(int iVel = 0;iVel < numPointsBuffer;iVel++){
288  // std::cout << "Vel(" << iVel << ") = ("
289  // << V[iVel] << "," << V[iVel+numPointsBuffer]
290  // << ")" << std::endl;
291  // std::cout << "PRESSURE: " << p[iVel] << std::endl;
292  // std::cout << "FLUX_RHO: " << inviscidFluxBuffer[iVel] << std::endl;
293  // std::cout << "FLUX_RHOV1: " << inviscidFluxBuffer[iVel+numPointsBuffer] << std::endl;
294  // std::cout << "FLUX_RHOV2: " << inviscidFluxBuffer[iVel+2*numPointsBuffer] << std::endl;
295  // std::cout << "FLUX_E: " << inviscidFluxBuffer[iVel+3*numPointsBuffer] << std::endl;
296  // }
297  //}
298 
299  double expectedMassFlux = 1.0;
300  std::vector<double> expectedMomFlux(numDim,1.0);
301  expectedMomFlux[iDim] = 2.6;
302  double expectedEnergyFlux = 5.6+.5*numDim;
303 
304  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
305  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
306  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
307  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
308  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
309  size_t bufferIndex = jBufferIndex + iIndex;
310  double actualMassFlux = inviscidFluxBuffer[bufferIndex];
311  if(std::abs(actualMassFlux - expectedMassFlux) > eps){
312  parTestWENO = false;
313  std::cout << "MassFlux(actual,expected) = (" << actualMassFlux
314  << "," << expectedMassFlux << ")" << std::endl;
315  }
316  for(int vDim = 0;vDim < numDim;vDim++){
317  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
318  double actualMomFlux = inviscidFluxBuffer[vIndex];
319  if(std::abs(actualMomFlux - expectedMomFlux[vDim]) > eps) {
320  parTestWENO = false;
321  std::cout << "MomFlux(actual,expected) = (" << actualMomFlux
322  << "," << expectedMomFlux[vDim] << ")" << std::endl;
323  }
324  }
325  size_t bufferIndexE = bufferIndex + (numDim+1)*numPointsBuffer;
326  double actualEnergyFlux = inviscidFluxBuffer[bufferIndexE];
327  if(std::abs(actualEnergyFlux- expectedEnergyFlux) > eps){
328  parTestWENO = false;
329  std::cout << "EnergyFlux(actual,expected) = (" << actualEnergyFlux
330  << "," << expectedEnergyFlux << ")" << std::endl;
331  }
332  }
333  }
334  }
335 
336  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Equal:FluxCheck",parTestWENO);
337 
338  if(!parTestWENO)
339  return;
340 
341  int statusWENO = rhsWENO.ApplyWENO(iDim,0);
342  if(statusWENO) {
343  parTestWENO = false;
344  std::cout << "ApplyWENO returned error code " << statusWENO << std::endl;
345  }
346 
347  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Equal:ErrorApplyWENO",parTestWENO);
348 
349  if(!parTestWENO)
350  return;
351 
352  // Check dFluxBuffer -- should be zero everywhere since state was
353  // initialized constant everywhere
354  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
355  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
356  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
357  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
358  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
359  size_t bufferIndex = jBufferIndex + iIndex;
360  for(int eqn = 0; eqn < numDim+2; eqn++){
361  size_t eIndex = bufferIndex + eqn*numPointsBuffer;
362  if(std::abs(dFluxBuffer[eIndex]) > eps)
363  {
364  parTestWENO = false;
365  std::cout << "Failed dFluxBuffer " << dFluxBuffer[eIndex] << std::endl;
366  }
367  }
368  }
369  }
370  }
371 
372  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Equal:dFluxCheck",parTestWENO);
373 
374  if(!parTestWENO)
375  return;
376 
377  } // Loop over dimensions
378 
379  } // Loop over numDim = 2,3
380 
381  // --- Jump ---
382  // Tests that when states have initial discontinuity, dFluxBuffer is
383  // computed correctly
384 
385  for(int numDim = 2;numDim < 4;numDim++){
386 
387  std::cout << "TestWENO:ApplyWENO:Jump:Testing WENO in " << numDim << " dimensions." << std::endl;
388 
389  for (int jumpDim=0; jumpDim<numDim; jumpDim++) {
390  grid_t simGrid;
391  operator_t sbpOperator;
392  halo_t &simHalo(simGrid.Halo());
393  state_t simState;
394  state_t paramState;
395  rhs_t rhsWENO;
396  weno_t &coeffsWENO(rhsWENO.WENOCoefficients());
397 
398  bool parTestWENO = true;
399  int numPoints = 21;
400  double length = 20.0;
401  std::vector<size_t> gridSizes(numDim,numPoints);
402  double dx = length/(numPoints - 1);
403 
404  if(testfixtures::CreateSimulationFixtures(simGrid,simHalo,sbpOperator,coeffsWENO,
405  gridSizes,2,testComm,&messageStream)){
406  parTestWENO = false;
407  }
408 
409  if(simGrid.GenerateCoordinates(messageStream))
410  parTestWENO = false;
411 
412  std::vector<double> &gridCoords(simGrid.CoordinateData());
413 
414  fixtures::CommunicatorType *gridCommPtr;
415  std::vector<int> cartNeighbors;
416  gridCommPtr = &simGrid.Communicator();
417  gridCommPtr->CartNeighbors(cartNeighbors);
418 
419  // Fill halo regions with coordinates
420  int xyzMessageId = simHalo.CreateMessageBuffers(numDim);
421  simHalo.PackMessageBuffers(xyzMessageId,0,numDim,&gridCoords[0]);
422 
423  simHalo.SendMessage(0,cartNeighbors,testComm);
424  simHalo.ReceiveMessage(0,cartNeighbors,testComm);
425  simHalo.UnPackMessageBuffers(0,0,numDim,&gridCoords[0]);
426 
427  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Jump:InitFixtures",parTestWENO);
428 
429  if(!parTestWENO)
430  return;
431 
432  pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());
433  pcpp::IndexIntervalType &partitionInterval(simGrid.PartitionInterval());
434  std::vector<size_t> &bufferSizes(simGrid.BufferSizes());
435 
436  if(testfixtures::euler::SetupEulerState(simGrid,simState,paramState,0)){
437  parTestWENO = false;
438  }
439 
440  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Jump:SetupState",parTestWENO);
441 
442  if(!parTestWENO)
443  return;
444 
445  std::cout << messageStream.str() << std::endl;
446  std::cout << "Testing WENO jump in " << jumpDim << " dimension." << std::endl;
447 
448  size_t numPointsBuffer = simGrid.BufferSize();
449 
450  double t = 0.0;
451  double gamma = 1.4;
452  double dt = .001;
453  double cfl = 1.0;
454  double Re = 0.0;
455  std::vector<double> numbers(4,0.0);
456  int flags = USEWENO;
457 
458  // Add some fields that are not usual
459  paramState.Create(numPointsBuffer,0);
460 
461  paramState.SetFieldBuffer("gamma",&gamma);
462  paramState.SetFieldBuffer("inputDT",&dt);
463  paramState.SetFieldBuffer("inputCFL",&cfl);
464  paramState.SetFieldBuffer("refRe",&Re);
465  paramState.SetFieldBuffer("Flags",&flags);
466  paramState.SetFieldBuffer("Numbers",&numbers[0]);
467 
468  std::vector<double> rho(numPointsBuffer,-1000.0);
469  std::vector<double> rhoV(numDim*numPointsBuffer,-1000.0);
470  std::vector<double> rhoE(numPointsBuffer,-1000.0);
471 
472  // "left" state (edges of domain in direction jumpDim)
473  double rhoL = 1.0;
474  double rhoVL = 2.0;
475  double rhoEL = 4.0 + 2.0*numDim;
476 
477  // "right" state (center of domain in direction jumpDim)
478  double rhoR = 2.0;
479  double rhoVR = 6.0;
480  double rhoER = 8.0 + 9.0*numDim;
481 
482  // locations of jumps, using physical coordinates
483  double loc1 = 0.25*length;
484  double loc2 = 0.75*length;
485 
486  // locations of jumps, using physical coordinates
487  double checkLoc1 = loc1;
488  double checkLoc2 = loc2 + dx;
489 
490  size_t iStart = partitionBufferInterval[0].first;
491  size_t iEnd = partitionBufferInterval[0].second;
492  size_t jStart = partitionBufferInterval[1].first;
493  size_t jEnd = partitionBufferInterval[1].second;
494 
495  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
496  // computed correctly
497  size_t kStart = 0;
498  size_t kEnd = 0;
499  if (numDim == 3) {
500  kStart = partitionBufferInterval[2].first;
501  kEnd = partitionBufferInterval[2].second;
502  }
503 
504  for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
505  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
506  for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
507  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
508  for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
509  size_t bufferIndex = jBufferIndex + iIndex;
510 
511  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
512 
513  if (cloc < loc1 || cloc > loc2) {
514  rho[bufferIndex] = rhoL;
515  for (int vDim=0; vDim<numDim; vDim++) {
516  rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVL;
517  }
518  rhoE[bufferIndex] = rhoEL;
519  } else {
520  rho[bufferIndex] = rhoR;
521  for (int vDim=0; vDim<numDim; vDim++) {
522  rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVR;
523  }
524  rhoE[bufferIndex] = rhoER;
525  }
526  }
527  }
528  }
529 
530  std::vector<double> dRho(numPointsBuffer,-1000.0);
531  std::vector<double> dRhoV(numDim*numPointsBuffer,-1000.0);
532  std::vector<double> dRhoE(numPointsBuffer,-1000.0);
533 
534  std::vector<double> T(numPointsBuffer,0.0);
535  std::vector<double> p(numPointsBuffer,0.0);
536  std::vector<double> V(numDim*numPointsBuffer,0.0);
537 
538  std::vector<double> rhom1(numPointsBuffer,0.0);
539 
540  simState.SetFieldBuffer("simTime",&t);
541  simState.SetFieldBuffer("rho",&rho[0]);
542  simState.SetFieldBuffer("rhoV",&rhoV[0]);
543  simState.SetFieldBuffer("rhoE",&rhoE[0]);
544 
545  simState.SetFieldBuffer("pressure",&p[0]);
546  simState.SetFieldBuffer("temperature",&T[0]);
547  simState.SetFieldBuffer("velocity",&V[0]);
548  simState.SetFieldBuffer("rhom1",&rhom1[0]);
549 
550  state_t rhsState(simState);
551 
552  rhsState.SetFieldBuffer("rho",&dRho[0]);
553  rhsState.SetFieldBuffer("rhoV",&dRhoV[0]);
554  rhsState.SetFieldBuffer("rhoE",&dRhoE[0]);
555 
556  // State/Soln initialization
557  std::vector<double> inParams;
558  std::vector<int> inFlags;
559 
560  if(rhsWENO.Initialize(simGrid,simState,paramState,sbpOperator)){
561  parTestWENO = false;
562  }
563  std::vector<int> numThreadsDim;
564  if(simGrid.SetupThreads(numThreadsDim)){
565  std::cout << "Setting up threads for grid object failed." << std::endl;
566  return;
567  }
568  rhsWENO.InitThreadIntervals();
569 
570  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Jump:InitRHS",parTestWENO);
571 
572  if(!parTestWENO)
573  return;
574 
575  if(rhsWENO.SetupRHS(simGrid,simState,rhsState)){
576  parTestWENO = false;
577  }
578 
579  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Jump:SetupRHS",parTestWENO);
580 
581  if(!parTestWENO)
582  return;
583 
584  //std::cout << "RHO BEFORE: " << std::endl;
585  //pcpp::io::DumpContents(std::cout,rho," ");
586  //std::cout << std::endl;
587 
588  // Exchange the state
589  rhsWENO.ExchangeState(0);
590 
591  //std::cout << "RHO AFTER: " << std::endl;
592  //pcpp::io::DumpContents(std::cout,rho," ");
593  //std::cout << std::endl;
594 
595  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
596  // computed correctly
597  size_t kFinal = 1;
598  if (numDim == 3) {
599  kFinal = bufferSizes[2];
600  }
601 
602  for (size_t kIndex = 0; kIndex < kFinal && parTestWENO; kIndex++) {
603  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
604  for (size_t jIndex = 0; jIndex < bufferSizes[1] && parTestWENO; jIndex++) {
605  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
606  for (size_t iIndex = 0; iIndex < bufferSizes[0] && parTestWENO; iIndex++) {
607  size_t bufferIndex = jBufferIndex + iIndex;
608 
609  int c = 0;
610  if (iIndex < iStart || iIndex > iEnd) c++;
611  if (jIndex < jStart || jIndex > jEnd) c++;
612  if (numDim == 3 && (kIndex < kStart || kIndex > kEnd)) c++;
613  if (c != 1) continue; // then we are not in a halo region
614 
615  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
616 
617  if (cloc < loc1 || cloc > loc2) {
618  if (std::abs(rho[bufferIndex] - rhoL) > eps) {
619  parTestWENO = false;
620  std::cout << "Failed state exchange, rho " << rho[bufferIndex] << std::endl;
621  }
622  for (int vDim=0; vDim<numDim; vDim++) {
623  if (std::abs(rhoV[bufferIndex + vDim*numPointsBuffer] - rhoVL) > eps) {
624  parTestWENO = false;
625  std::cout << "Failed state exchange, rhoV" << vDim << " "
626  << rhoV[bufferIndex + vDim*numPointsBuffer] << std::endl;
627  }
628  }
629  if (std::abs(rhoE[bufferIndex] - rhoEL) > eps) {
630  parTestWENO = false;
631  std::cout << "Failed state exchange, rhoE " << rhoE[bufferIndex] << std::endl;
632  }
633  } else {
634  if (std::abs(rho[bufferIndex] - rhoR) > eps) {
635  parTestWENO = false;
636  std::cout << "Failed state exchange, rho " << rho[bufferIndex] << std::endl;
637  }
638  for (int vDim=0; vDim<numDim; vDim++) {
639  if (std::abs(rhoV[bufferIndex + vDim*numPointsBuffer] - rhoVR) > eps) {
640  parTestWENO = false;
641  std::cout << "Failed state exchange, rhoV" << vDim << " "
642  << rhoV[bufferIndex + vDim*numPointsBuffer] << std::endl;
643  }
644  }
645  if (std::abs(rhoE[bufferIndex] - rhoER) > eps) {
646  parTestWENO = false;
647  std::cout << "Failed state exchange, rhoE " << rhoE[bufferIndex] << std::endl;
648  }
649  }
650 
651  if (!parTestWENO) {
652  std::cout << "iIndex = " << iIndex << ", jIndex = " << jIndex
653  << ", kIndex = " << kIndex << std::endl;
654  std::cout << "iStart = " << iStart << ", jStart = " << jStart
655  << ", kStart = " << kStart << std::endl;
656  std::cout << "iEnd = " << iEnd << ", jEnd = " << jEnd
657  << ", kEnd = " << kEnd << std::endl;
658  }
659  }
660  }
661  }
662 
663  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Jump:ExchangeState",parTestWENO);
664 
665  if(!parTestWENO)
666  return;
667 
668  double *inviscidFluxBuffer(rhsWENO.InviscidFluxBuffer());
669  double *dFluxBuffer(rhsWENO.DFluxBuffer());
670 
671  for(int iDim = 0;iDim < numDim;iDim++){
672 
673  // Running WENO in dimension (iDim)
674  std::cout << "Checking WENO in dimension(" << iDim << ")"
675  << std::endl;
676 
677  // initialize dFluxBuffer so test can't pass because data was initialized to zero
678  for (size_t i=0; i<numPointsBuffer; i++) {
679  for (int j=0; j<numDim+2; j++) {
680  size_t bufferIndex = j*numPointsBuffer + i;
681  dFluxBuffer[bufferIndex] = -1000;
682  }
683  }
684 
685  rhsWENO.ComputeDV(0);
686  rhsWENO.InviscidFlux(iDim,0);
687  rhsWENO.ExchangeInviscidFlux(0);
688 
689  //if(numDim == 2){
690  // std::cout << "Velocity = " << std::endl;
691  // for(int iVel = 0;iVel < numPointsBuffer;iVel++){
692  // std::cout << "Vel(" << iVel << ") = ("
693  // << V[iVel] << "," << V[iVel+numPointsBuffer]
694  // << ")" << std::endl;
695  // std::cout << "PRESSURE: " << p[iVel] << std::endl;
696  // std::cout << "FLUX_RHO: " << inviscidFluxBuffer[iVel] << std::endl;
697  // std::cout << "FLUX_RHOV1: " << inviscidFluxBuffer[iVel+numPointsBuffer] << std::endl;
698  // std::cout << "FLUX_RHOV2: " << inviscidFluxBuffer[iVel+2*numPointsBuffer] << std::endl;
699  // std::cout << "FLUX_E: " << inviscidFluxBuffer[iVel+3*numPointsBuffer] << std::endl;
700  // }
701  //}
702 
703  // "left" state (edges of domain in direction jumpDim)
704  double fluxRhoL = 2.0;
705  double fluxRhoVMainL = 5.6;
706  double fluxRhoVAltL = 4.0;
707  double fluxRhoEL = 11.2 + 4.0*numDim;
708 
709  // "right" state (center of domain in direction jumpDim)
710  double fluxRhoR = 6.0;
711  double fluxRhoVMainR = 21.2;
712  double fluxRhoVAltR = 18.0;
713  double fluxRhoER = 33.6 + 27.0*numDim;
714 
715  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
716  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
717  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
718  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
719  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
720  size_t bufferIndex = jBufferIndex + iIndex;
721 
722  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
723 
724  if (cloc < loc1 || cloc > loc2) {
725  if(std::abs(inviscidFluxBuffer[bufferIndex] - fluxRhoL) > eps) {
726  parTestWENO = false;
727  std::cout << "Flux check failed: rho = " << inviscidFluxBuffer[bufferIndex]
728  << ", expected " << fluxRhoL << std::endl;
729  }
730  for(int vDim = 0;vDim < numDim;vDim++){
731  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
732  if(vDim == iDim){
733  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVMainL) > eps) {
734  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
735  << ", expected " << fluxRhoVMainL << std::endl;
736  parTestWENO = false;
737  }
738  } else {
739  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVAltL) > eps) {
740  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
741  << ", expected " << fluxRhoVAltL << std::endl;
742  parTestWENO = false;
743  }
744  }
745  }
746  if(std::abs(inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] -
747  fluxRhoEL) > eps) {
748  std::cout << "Flux check failed: rhoE = "
749  << inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
750  << ", expected " << fluxRhoEL << std::endl;
751  parTestWENO = false;
752  }
753  } else {
754  if(std::abs(inviscidFluxBuffer[bufferIndex] - fluxRhoR) > eps) {
755  parTestWENO = false;
756  std::cout << "Flux check failed: rho = " << inviscidFluxBuffer[bufferIndex]
757  << ", expected " << fluxRhoR << std::endl;
758  }
759  for(int vDim = 0;vDim < numDim;vDim++){
760  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
761  if(vDim == iDim){
762  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVMainR) > eps) {
763  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
764  << ", expected " << fluxRhoVMainR << std::endl;
765  parTestWENO = false;
766  }
767  } else {
768  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVAltR) > eps) {
769  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
770  << ", expected " << fluxRhoVAltR << std::endl;
771  parTestWENO = false;
772  }
773  }
774  }
775  if(std::abs(inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] -
776  fluxRhoER) > eps) {
777  std::cout << "Flux check failed: rhoE = "
778  << inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
779  << ", expected " << fluxRhoER << std::endl;
780  parTestWENO = false;
781  }
782  } // cloc if statement
783 
784  if (!parTestWENO) {
785  std::cout << "i = " << iIndex << ", j = " << jIndex << ", k = " << kIndex
786  << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
787  std::cout << "rho = " << rho[bufferIndex];
788  for (int vDim=0; vDim<numDim; vDim++) {
789  std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
790  }
791  std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
792  }
793  }
794  }
795  }
796 
797  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Jump:FluxCheck",parTestWENO);
798 
799  if(!parTestWENO)
800  return;
801 
802  int statusWENO = rhsWENO.ApplyWENO(iDim,0);
803  if(statusWENO) {
804  parTestWENO = false;
805  std::cout << "ApplyWENO returned error code " << statusWENO << std::endl;
806  }
807 
808  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Jump:ErrorApplyWENO",parTestWENO);
809 
810  if(!parTestWENO)
811  return;
812 
813  // setup expected values
814  double dFluxRho = fluxRhoR - fluxRhoL;
815  double dFluxRhoVMain = fluxRhoVMainR - fluxRhoVMainL;
816  double dFluxRhoVAlt = fluxRhoVAltR - fluxRhoVAltL;
817  double dFluxRhoE = fluxRhoER - fluxRhoEL;
818 
819  // Check dFluxBuffer -- should be zero everywhere except at initial discontinuity
820  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
821  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
822  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
823  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
824  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
825  size_t bufferIndex = jBufferIndex + iIndex;
826 
827  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
828 
829  if (iDim == jumpDim && std::abs(cloc - checkLoc1) < eps) {
830  if(std::abs(dFluxBuffer[bufferIndex] - dFluxRho) > epsLarge) {
831  parTestWENO = false;
832  std::cout << "dFlux check failed: rho = " << dFluxBuffer[bufferIndex]
833  << ", expected " << dFluxRho << std::endl;
834  }
835  for(int vDim = 0;vDim < numDim;vDim++){
836  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
837  if(vDim == iDim){
838  if(std::abs(dFluxBuffer[vIndex] - dFluxRhoVMain) > epsLarge) {
839  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
840  << ", expected " << dFluxRhoVMain << std::endl;
841  parTestWENO = false;
842  }
843  } else {
844  if(std::abs(dFluxBuffer[vIndex] - dFluxRhoVAlt) > epsLarge) {
845  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
846  << ", expected " << dFluxRhoVAlt << std::endl;
847  parTestWENO = false;
848  }
849  }
850  }
851  if(std::abs(dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] -
852  dFluxRhoE) > epsLarge) {
853  std::cout << "dFlux check failed: rhoE = "
854  << dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
855  << ", expected " << dFluxRhoE << std::endl;
856  parTestWENO = false;
857  }
858  } else if (iDim == jumpDim && std::abs(cloc - checkLoc2) < eps) {
859  if(std::abs(dFluxBuffer[bufferIndex] + dFluxRho) > epsLarge) {
860  parTestWENO = false;
861  std::cout << "dFlux check failed: rho = " << dFluxBuffer[bufferIndex]
862  << ", expected " << -dFluxRho << std::endl;
863  }
864  for(int vDim = 0;vDim < numDim;vDim++){
865  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
866  if(vDim == iDim){
867  if(std::abs(dFluxBuffer[vIndex] + dFluxRhoVMain) > epsLarge) {
868  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
869  << ", expected " << -dFluxRhoVMain << std::endl;
870  parTestWENO = false;
871  }
872  } else {
873  if(std::abs(dFluxBuffer[vIndex] + dFluxRhoVAlt) > epsLarge) {
874  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
875  << ", expected " << -dFluxRhoVAlt << std::endl;
876  parTestWENO = false;
877  }
878  }
879  }
880  if(std::abs(dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] +
881  dFluxRhoE) > epsLarge) {
882  std::cout << "dFlux check failed: rhoE = "
883  << dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
884  << ", expected " << -dFluxRhoE << std::endl;
885  parTestWENO = false;
886  }
887  } else if (iDim == jumpDim && (std::abs(cloc - checkLoc1) < stenWidth*dx+eps)
888  || (std::abs(cloc - checkLoc2) < stenWidth*dx+eps)) {
889  // near discontinuities WENO can have larger values of "0"
890  // (i.e. numerical dissipation)
891  for(int eqn = 0; eqn < numDim+2; eqn++){
892  size_t eIndex = bufferIndex + eqn*numPointsBuffer;
893  if(std::abs(dFluxBuffer[eIndex]) > epsLarge)
894  {
895  parTestWENO = false;
896  std::cout << "Failed dFluxBuffer " << dFluxBuffer[eIndex] << std::endl;
897  }
898  }
899  } else {
900  for(int eqn = 0; eqn < numDim+2; eqn++){
901  size_t eIndex = bufferIndex + eqn*numPointsBuffer;
902  if(std::abs(dFluxBuffer[eIndex]) > eps)
903  {
904  parTestWENO = false;
905  std::cout << "Failed dFluxBuffer " << dFluxBuffer[eIndex] << std::endl;
906  }
907  }
908  } // checkLoc if statement
909 
910  if (!parTestWENO) {
911  std::cout << "i = " << iIndex << ", j = " << jIndex << ", k = " << kIndex
912  << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
913  std::cout << "checkLoc1 = " << checkLoc1 << ", checkLoc2 = " << checkLoc2 << std::endl;
914  std::cout << "rho = " << rho[bufferIndex];
915  for (int vDim=0; vDim<numDim; vDim++) {
916  std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
917  }
918  std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
919  }
920  }
921  }
922  }
923 
924  parallelUnitResults.UpdateResult("WENO:ApplyWENO:Jump:dFluxCheck",parTestWENO);
925 
926  if(!parTestWENO)
927  return;
928 
929  } // Loop over dimensions
930  } // jumpDim loop
931  } // Loop over numDim = 2,3
932 
933  // --- TransonicJump ---
934  // Tests that when states have initial discontinuity, dFluxBuffer is
935  // computed correctly when there is a transonic rarefaction and entropy fix
936  // is needed for Roe scheme
937 
938  for(int numDim = 2;numDim < 4;numDim++){
939 
940  std::cout << "TestWENO:ApplyWENO:TransonicJump:Testing WENO in " << numDim << " dimensions." << std::endl;
941 
942  for (int jumpDim=0; jumpDim<numDim; jumpDim++) {
943  grid_t simGrid;
944  operator_t sbpOperator;
945  halo_t &simHalo(simGrid.Halo());
946  state_t simState;
947  state_t paramState;
948  rhs_t rhsWENO;
949  weno_t &coeffsWENO(rhsWENO.WENOCoefficients());
950 
951  bool parTestWENO = true;
952  int numPoints = 21;
953  double length = 20.0;
954  std::vector<size_t> gridSizes(numDim,numPoints);
955  double dx = length/(numPoints - 1);
956 
957  if(testfixtures::CreateSimulationFixtures(simGrid,simHalo,sbpOperator,coeffsWENO,
958  gridSizes,2,testComm,&messageStream)){
959  parTestWENO = false;
960  }
961 
962  if(simGrid.GenerateCoordinates(messageStream))
963  parTestWENO = false;
964 
965  std::vector<double> &gridCoords(simGrid.CoordinateData());
966 
967  fixtures::CommunicatorType *gridCommPtr;
968  std::vector<int> cartNeighbors;
969  gridCommPtr = &simGrid.Communicator();
970  gridCommPtr->CartNeighbors(cartNeighbors);
971 
972  // Fill halo regions with coordinates
973  int xyzMessageId = simHalo.CreateMessageBuffers(numDim);
974  simHalo.PackMessageBuffers(xyzMessageId,0,numDim,&gridCoords[0]);
975 
976  simHalo.SendMessage(0,cartNeighbors,testComm);
977  simHalo.ReceiveMessage(0,cartNeighbors,testComm);
978  simHalo.UnPackMessageBuffers(0,0,numDim,&gridCoords[0]);
979 
980  parallelUnitResults.UpdateResult("WENO:ApplyWENO:TransonicJump:InitFixtures",parTestWENO);
981 
982  if(!parTestWENO)
983  return;
984 
985  pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());
986  pcpp::IndexIntervalType &partitionInterval(simGrid.PartitionInterval());
987  std::vector<size_t> &bufferSizes(simGrid.BufferSizes());
988 
989  if(testfixtures::euler::SetupEulerState(simGrid,simState,paramState,0)){
990  parTestWENO = false;
991  }
992 
993  parallelUnitResults.UpdateResult("WENO:ApplyWENO:TransonicJump:SetupState",parTestWENO);
994 
995  if(!parTestWENO)
996  return;
997 
998  std::cout << messageStream.str() << std::endl;
999  std::cout << "Testing WENO jump in " << jumpDim << " dimension." << std::endl;
1000 
1001  size_t numPointsBuffer = simGrid.BufferSize();
1002 
1003  double t = 0.0;
1004  double gamma = 1.4;
1005  double dt = .001;
1006  double cfl = 1.0;
1007  double Re = 0.0;
1008  std::vector<double> numbers(4,0.0);
1009  int flags = USEWENO;
1010 
1011  // Add some fields that are not usual
1012  paramState.Create(numPointsBuffer,0);
1013 
1014  paramState.SetFieldBuffer("gamma",&gamma);
1015  paramState.SetFieldBuffer("inputDT",&dt);
1016  paramState.SetFieldBuffer("inputCFL",&cfl);
1017  paramState.SetFieldBuffer("refRe",&Re);
1018  paramState.SetFieldBuffer("Flags",&flags);
1019  paramState.SetFieldBuffer("Numbers",&numbers[0]);
1020 
1021  std::vector<double> rho(numPointsBuffer,-1000.0);
1022  std::vector<double> rhoV(numDim*numPointsBuffer,-1000.0);
1023  std::vector<double> rhoE(numPointsBuffer,-1000.0);
1024 
1025  // "left" state (edges of domain in direction jumpDim)
1026  // corresponding primitive variables:
1027  // rho = 1, u = 1, v = w = 0, p = 1.6
1028  double rhoL = 1.0;
1029  double rhoVL = 1.0;
1030  double rhoEL = 4.5;
1031 
1032  // "right" state (center of domain in direction jumpDim)
1033  // corresponding primitive variables:
1034  // rho = 0.1, u = 0, v = w = 0, p = 0.4
1035  double rhoR = 0.1;
1036  double rhoVR = 0.0;
1037  double rhoER = 1.0;
1038 
1039  // locations of jumps, using physical coordinates
1040  double loc1 = 0.25*length;
1041  double loc2 = 0.75*length;
1042 
1043  // locations of jumps, using physical coordinates
1044  double checkLoc1 = loc1;
1045  double checkLoc2 = loc2 + dx;
1046 
1047  size_t iStart = partitionBufferInterval[0].first;
1048  size_t iEnd = partitionBufferInterval[0].second;
1049  size_t jStart = partitionBufferInterval[1].first;
1050  size_t jEnd = partitionBufferInterval[1].second;
1051 
1052  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
1053  // computed correctly
1054  size_t kStart = 0;
1055  size_t kEnd = 0;
1056  if (numDim == 3) {
1057  kStart = partitionBufferInterval[2].first;
1058  kEnd = partitionBufferInterval[2].second;
1059  }
1060 
1061  for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
1062  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
1063  for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
1064  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
1065  for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
1066  size_t bufferIndex = jBufferIndex + iIndex;
1067 
1068  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
1069 
1070  if (cloc < loc1 || cloc > loc2) {
1071  rho[bufferIndex] = rhoL;
1072  for (int vDim=0; vDim<numDim; vDim++) {
1073  if (vDim == jumpDim) rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVL;
1074  else rhoV[bufferIndex + vDim*numPointsBuffer] = 0.0;
1075  }
1076  rhoE[bufferIndex] = rhoEL;
1077  } else {
1078  rho[bufferIndex] = rhoR;
1079  for (int vDim=0; vDim<numDim; vDim++) {
1080  if (vDim == jumpDim) rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVR;
1081  else rhoV[bufferIndex + vDim*numPointsBuffer] = 0.0;
1082  }
1083  rhoE[bufferIndex] = rhoER;
1084  }
1085  }
1086  }
1087  }
1088 
1089  std::vector<double> dRho(numPointsBuffer,-1000.0);
1090  std::vector<double> dRhoV(numDim*numPointsBuffer,-1000.0);
1091  std::vector<double> dRhoE(numPointsBuffer,-1000.0);
1092 
1093  std::vector<double> T(numPointsBuffer,0.0);
1094  std::vector<double> p(numPointsBuffer,0.0);
1095  std::vector<double> V(numDim*numPointsBuffer,0.0);
1096 
1097  std::vector<double> rhom1(numPointsBuffer,0.0);
1098 
1099  simState.SetFieldBuffer("simTime",&t);
1100  simState.SetFieldBuffer("rho",&rho[0]);
1101  simState.SetFieldBuffer("rhoV",&rhoV[0]);
1102  simState.SetFieldBuffer("rhoE",&rhoE[0]);
1103 
1104  simState.SetFieldBuffer("pressure",&p[0]);
1105  simState.SetFieldBuffer("temperature",&T[0]);
1106  simState.SetFieldBuffer("velocity",&V[0]);
1107  simState.SetFieldBuffer("rhom1",&rhom1[0]);
1108 
1109  state_t rhsState(simState);
1110 
1111  rhsState.SetFieldBuffer("rho",&dRho[0]);
1112  rhsState.SetFieldBuffer("rhoV",&dRhoV[0]);
1113  rhsState.SetFieldBuffer("rhoE",&dRhoE[0]);
1114 
1115  // State/Soln initialization
1116  std::vector<double> inParams;
1117  std::vector<int> inFlags;
1118 
1119  if(rhsWENO.Initialize(simGrid,simState,paramState,sbpOperator)){
1120  parTestWENO = false;
1121  }
1122  std::vector<int> numThreadsDim;
1123  if(simGrid.SetupThreads(numThreadsDim)){
1124  std::cout << "Setting up threads for grid object failed." << std::endl;
1125  return;
1126  }
1127  rhsWENO.InitThreadIntervals();
1128 
1129  parallelUnitResults.UpdateResult("WENO:ApplyWENO:TransonicJump:InitRHS",parTestWENO);
1130 
1131  if(!parTestWENO)
1132  return;
1133 
1134  if(rhsWENO.SetupRHS(simGrid,simState,rhsState)){
1135  parTestWENO = false;
1136  }
1137 
1138  parallelUnitResults.UpdateResult("WENO:ApplyWENO:TransonicJump:SetupRHS",parTestWENO);
1139 
1140  if(!parTestWENO)
1141  return;
1142 
1143  //std::cout << "RHO BEFORE: " << std::endl;
1144  //pcpp::io::DumpContents(std::cout,rho," ");
1145  //std::cout << std::endl;
1146 
1147  // Exchange the state
1148  rhsWENO.ExchangeState(0);
1149 
1150  //std::cout << "RHO AFTER: " << std::endl;
1151  //pcpp::io::DumpContents(std::cout,rho," ");
1152  //std::cout << std::endl;
1153 
1154  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
1155  // computed correctly
1156  size_t kFinal = 1;
1157  if (numDim == 3) {
1158  kFinal = bufferSizes[2];
1159  }
1160 
1161  for (size_t kIndex = 0; kIndex < kFinal && parTestWENO; kIndex++) {
1162  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
1163  for (size_t jIndex = 0; jIndex < bufferSizes[1] && parTestWENO; jIndex++) {
1164  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
1165  for (size_t iIndex = 0; iIndex < bufferSizes[0] && parTestWENO; iIndex++) {
1166  size_t bufferIndex = jBufferIndex + iIndex;
1167 
1168  int c = 0;
1169  if (iIndex < iStart || iIndex > iEnd) c++;
1170  if (jIndex < jStart || jIndex > jEnd) c++;
1171  if (numDim == 3 && (kIndex < kStart || kIndex > kEnd)) c++;
1172  if (c != 1) continue; // then we are not in a halo region
1173 
1174  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
1175 
1176  if (cloc < loc1 || cloc > loc2) {
1177  if (std::abs(rho[bufferIndex] - rhoL) > eps) {
1178  parTestWENO = false;
1179  std::cout << "Failed state exchange, rho " << rho[bufferIndex] << std::endl;
1180  }
1181  for (int vDim=0; vDim<numDim; vDim++) {
1182  if (vDim == jumpDim) {
1183  if (std::abs(rhoV[bufferIndex + vDim*numPointsBuffer] - rhoVL) > eps) {
1184  parTestWENO = false;
1185  std::cout << "Failed state exchange, rhoV" << vDim << " "
1186  << rhoV[bufferIndex + vDim*numPointsBuffer] << std::endl;
1187  }
1188  } else {
1189  if (std::abs(rhoV[bufferIndex + vDim*numPointsBuffer]) > eps) {
1190  parTestWENO = false;
1191  std::cout << "Failed state exchange, rhoV" << vDim << " "
1192  << rhoV[bufferIndex + vDim*numPointsBuffer] << std::endl;
1193  }
1194  }
1195  }
1196  if (std::abs(rhoE[bufferIndex] - rhoEL) > eps) {
1197  parTestWENO = false;
1198  std::cout << "Failed state exchange, rhoE " << rhoE[bufferIndex] << std::endl;
1199  }
1200  } else {
1201  if (std::abs(rho[bufferIndex] - rhoR) > eps) {
1202  parTestWENO = false;
1203  std::cout << "Failed state exchange, rho " << rho[bufferIndex] << std::endl;
1204  }
1205  for (int vDim=0; vDim<numDim; vDim++) {
1206  if (vDim == jumpDim) {
1207  if (std::abs(rhoV[bufferIndex + vDim*numPointsBuffer] - rhoVR) > eps) {
1208  parTestWENO = false;
1209  std::cout << "Failed state exchange, rhoV" << vDim << " "
1210  << rhoV[bufferIndex + vDim*numPointsBuffer] << std::endl;
1211  }
1212  } else {
1213  if (std::abs(rhoV[bufferIndex + vDim*numPointsBuffer]) > eps) {
1214  parTestWENO = false;
1215  std::cout << "Failed state exchange, rhoV" << vDim << " "
1216  << rhoV[bufferIndex + vDim*numPointsBuffer] << std::endl;
1217  }
1218  }
1219  }
1220  if (std::abs(rhoE[bufferIndex] - rhoER) > eps) {
1221  parTestWENO = false;
1222  std::cout << "Failed state exchange, rhoE " << rhoE[bufferIndex] << std::endl;
1223  }
1224  }
1225 
1226  if (!parTestWENO) {
1227  std::cout << "iIndex = " << iIndex << ", jIndex = " << jIndex
1228  << ", kIndex = " << kIndex << std::endl;
1229  std::cout << "iStart = " << iStart << ", jStart = " << jStart
1230  << ", kStart = " << kStart << std::endl;
1231  std::cout << "iEnd = " << iEnd << ", jEnd = " << jEnd
1232  << ", kEnd = " << kEnd << std::endl;
1233  }
1234  }
1235  }
1236  }
1237 
1238  parallelUnitResults.UpdateResult("WENO:ApplyWENO:TransonicJump:ExchangeState",parTestWENO);
1239 
1240  if(!parTestWENO)
1241  return;
1242 
1243  double *inviscidFluxBuffer(rhsWENO.InviscidFluxBuffer());
1244  double *dFluxBuffer(rhsWENO.DFluxBuffer());
1245 
1246  for(int iDim = 0;iDim < numDim;iDim++){
1247 
1248  // Running WENO in dimension (iDim)
1249  std::cout << "Checking WENO in dimension(" << iDim << ")"
1250  << std::endl;
1251 
1252  // initialize dFluxBuffer so test can't pass because data was initialized to zero
1253  for (size_t i=0; i<numPointsBuffer; i++) {
1254  for (int j=0; j<numDim+2; j++) {
1255  size_t bufferIndex = j*numPointsBuffer + i;
1256  dFluxBuffer[bufferIndex] = -1000;
1257  }
1258  }
1259 
1260  rhsWENO.ComputeDV(0);
1261  rhsWENO.InviscidFlux(iDim,0);
1262  rhsWENO.ExchangeInviscidFlux(0);
1263 
1264  //if(numDim == 2){
1265  // std::cout << "Velocity = " << std::endl;
1266  // for(int iVel = 0;iVel < numPointsBuffer;iVel++){
1267  // std::cout << "Vel(" << iVel << ") = ("
1268  // << V[iVel] << "," << V[iVel+numPointsBuffer]
1269  // << ")" << std::endl;
1270  // std::cout << "PRESSURE: " << p[iVel] << std::endl;
1271  // std::cout << "FLUX_RHO: " << inviscidFluxBuffer[iVel] << std::endl;
1272  // std::cout << "FLUX_RHOV1: " << inviscidFluxBuffer[iVel+numPointsBuffer] << std::endl;
1273  // std::cout << "FLUX_RHOV2: " << inviscidFluxBuffer[iVel+2*numPointsBuffer] << std::endl;
1274  // std::cout << "FLUX_E: " << inviscidFluxBuffer[iVel+3*numPointsBuffer] << std::endl;
1275  // }
1276  //}
1277 
1278  // "left" state (edges of domain in direction jumpDim)
1279  double fluxRhoL = 1.0;
1280  double fluxRhoVMainL = 2.6;
1281  double fluxRhoVAltL = 0.0;
1282  double fluxRhoEL = 6.1;
1283  double fluxRhoVTransverseL = 1.6;
1284 
1285  // "right" state (center of domain in direction jumpDim)
1286  double fluxRhoR = 0.0;
1287  double fluxRhoVMainR = 0.4;
1288  double fluxRhoVAltR = 0.0;
1289  double fluxRhoER = 0.0;
1290  double fluxRhoVTransverseR = 0.4;
1291 
1292  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
1293  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
1294  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
1295  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
1296  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
1297  size_t bufferIndex = jBufferIndex + iIndex;
1298 
1299  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
1300 
1301  if (iDim == jumpDim) {
1302  if (cloc < loc1 || cloc > loc2) {
1303  if(std::abs(inviscidFluxBuffer[bufferIndex] - fluxRhoL) > eps) {
1304  parTestWENO = false;
1305  std::cout << "Flux check failed: rho = " << inviscidFluxBuffer[bufferIndex]
1306  << ", expected " << fluxRhoL << std::endl;
1307  }
1308  for(int vDim = 0;vDim < numDim;vDim++){
1309  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
1310  if(vDim == iDim){
1311  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVMainL) > eps) {
1312  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
1313  << ", expected " << fluxRhoVMainL << std::endl;
1314  parTestWENO = false;
1315  }
1316  } else {
1317  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVAltL) > eps) {
1318  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
1319  << ", expected " << fluxRhoVAltL << std::endl;
1320  parTestWENO = false;
1321  }
1322  }
1323  }
1324  if(std::abs(inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] -
1325  fluxRhoEL) > eps) {
1326  std::cout << "Flux check failed: rhoE = "
1327  << inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
1328  << ", expected " << fluxRhoEL << std::endl;
1329  parTestWENO = false;
1330  }
1331  } else {
1332  if(std::abs(inviscidFluxBuffer[bufferIndex] - fluxRhoR) > eps) {
1333  parTestWENO = false;
1334  std::cout << "Flux check failed: rho = " << inviscidFluxBuffer[bufferIndex]
1335  << ", expected " << fluxRhoR << std::endl;
1336  }
1337  for(int vDim = 0;vDim < numDim;vDim++){
1338  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
1339  if(vDim == iDim){
1340  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVMainR) > eps) {
1341  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
1342  << ", expected " << fluxRhoVMainR << std::endl;
1343  parTestWENO = false;
1344  }
1345  } else {
1346  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVAltR) > eps) {
1347  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
1348  << ", expected " << fluxRhoVAltR << std::endl;
1349  parTestWENO = false;
1350  }
1351  }
1352  }
1353  if(std::abs(inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] -
1354  fluxRhoER) > eps) {
1355  std::cout << "Flux check failed: rhoE = "
1356  << inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
1357  << ", expected " << fluxRhoER << std::endl;
1358  parTestWENO = false;
1359  }
1360  } // cloc if statement
1361  } else {
1362  for(int eqn = 0; eqn < numDim+2; eqn++){
1363  size_t eIndex = bufferIndex + eqn*numPointsBuffer;
1364  if(eqn == iDim+1) { // momentum flux in transverse direction still includes pressure
1365  if (cloc < loc1 || cloc > loc2) {
1366  if(std::abs(inviscidFluxBuffer[eIndex] - fluxRhoVTransverseL) > eps)
1367  {
1368  parTestWENO = false;
1369  std::cout << "Flux check failed in transverse direction "
1370  << inviscidFluxBuffer[eIndex] << std::endl;
1371  }
1372  } else {
1373  if(std::abs(inviscidFluxBuffer[eIndex] - fluxRhoVTransverseR) > eps)
1374  {
1375  parTestWENO = false;
1376  std::cout << "Flux check failed in transverse direction "
1377  << inviscidFluxBuffer[eIndex] << std::endl;
1378  }
1379  }
1380  } else {
1381  if(std::abs(inviscidFluxBuffer[eIndex]) > eps)
1382  {
1383  parTestWENO = false;
1384  std::cout << "Flux check failed in transverse direction "
1385  << inviscidFluxBuffer[eIndex] << std::endl;
1386  }
1387  }
1388  }
1389  } // iDim == jumpDim if statement
1390 
1391  if (!parTestWENO) {
1392  std::cout << "i = " << iIndex << ", j = " << jIndex << ", k = " << kIndex
1393  << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
1394  std::cout << "rho = " << rho[bufferIndex];
1395  for (int vDim=0; vDim<numDim; vDim++) {
1396  std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
1397  }
1398  std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
1399  }
1400  }
1401  }
1402  }
1403 
1404  parallelUnitResults.UpdateResult("WENO:ApplyWENO:TransonicJump:FluxCheck",parTestWENO);
1405 
1406  if(!parTestWENO)
1407  return;
1408 
1409  int statusWENO = rhsWENO.ApplyWENO(iDim,0);
1410  if(statusWENO) {
1411  parTestWENO = false;
1412  std::cout << "ApplyWENO returned error code " << statusWENO << std::endl;
1413  }
1414 
1415  parallelUnitResults.UpdateResult("WENO:ApplyWENO:TransonicJump:ErrorApplyWENO",parTestWENO);
1416 
1417  if(!parTestWENO)
1418  return;
1419 
1421  //if (jumpDim == 0) {
1422  // size_t kBufferIndex = kStart*bufferSizes[0]*bufferSizes[1];
1423  // size_t jBufferIndex = jStart*bufferSizes[0] + kBufferIndex;
1424  // for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
1425  // size_t bufferIndex = jBufferIndex + iIndex;
1426 
1427  // double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
1428 
1429  // //std::cout << "i = " << iIndex << ", j = " << jStart << ", k = " << kStart
1430  // // << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
1431  // //std::cout << "checkLoc1 = " << checkLoc1 << ", checkLoc2 = " << checkLoc2 << std::endl;
1432  // //std::cout << "rho = " << rho[bufferIndex];
1433  // //for (int vDim=0; vDim<numDim; vDim++) {
1434  // // std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
1435  // //}
1436  // //std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
1437 
1438  // std::cout << "cloc = " << cloc << ", ";
1439 
1440  // std::cout << "dFluxRho = " << dFluxBuffer[bufferIndex];
1441  // for (int vDim=0; vDim<numDim; vDim++) {
1442  // std::cout << ", dFluxRhoV[" << vDim << "] = " << dFluxBuffer[bufferIndex + (vDim+1)*numPointsBuffer];
1443  // }
1444  // std::cout << ", dFluxRhoE = " << dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] << std::endl;
1445  // }
1446  //}
1447 
1448  // setup expected values
1449  // larger values
1450  double dFluxLargeRho = -1.149;
1451  double dFluxLargeRhoVMain = -2.130;
1452  double dFluxLargeRhoVAlt = 0.0;
1453  double dFluxLargeRhoE = -6.807;
1454 
1455  // smaller values
1456  double dFluxSmallRho = 0.149;
1457  double dFluxSmallRhoVMain = -0.0698;
1458  double dFluxSmallRhoVAlt = 0.0;
1459  double dFluxSmallRhoE = 0.707;
1460 
1461  // Check dFluxBuffer -- should be zero everywhere except at initial discontinuity
1462  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
1463  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
1464  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
1465  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
1466  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
1467  size_t bufferIndex = jBufferIndex + iIndex;
1468 
1469  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
1470 
1471  if (iDim == jumpDim && std::abs(cloc - checkLoc1) < eps) {
1472  if(std::abs(dFluxBuffer[bufferIndex] - dFluxLargeRho) > epsLarge) {
1473  parTestWENO = false;
1474  std::cout << "dFlux check failed: rho = " << dFluxBuffer[bufferIndex]
1475  << ", expected " << dFluxLargeRho << std::endl;
1476  }
1477  for(int vDim = 0;vDim < numDim;vDim++){
1478  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
1479  if(vDim == iDim){
1480  if(std::abs(dFluxBuffer[vIndex] - dFluxLargeRhoVMain) > epsLarge) {
1481  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
1482  << ", expected " << dFluxLargeRhoVMain << std::endl;
1483  parTestWENO = false;
1484  }
1485  } else {
1486  if(std::abs(dFluxBuffer[vIndex] - dFluxLargeRhoVAlt) > epsLarge) {
1487  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
1488  << ", expected " << dFluxLargeRhoVAlt << std::endl;
1489  parTestWENO = false;
1490  }
1491  }
1492  }
1493  if(std::abs(dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] -
1494  dFluxLargeRhoE) > epsLarge) {
1495  std::cout << "dFlux check failed: rhoE = "
1496  << dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
1497  << ", expected " << dFluxLargeRhoE << std::endl;
1498  parTestWENO = false;
1499  }
1500  } else if (iDim == jumpDim && std::abs(cloc - (checkLoc1-dx)) < eps) {
1501  if(std::abs(dFluxBuffer[bufferIndex] - dFluxSmallRho) > epsLarge) {
1502  parTestWENO = false;
1503  std::cout << "dFlux check failed: rho = " << dFluxBuffer[bufferIndex]
1504  << ", expected " << dFluxSmallRho << std::endl;
1505  }
1506  for(int vDim = 0;vDim < numDim;vDim++){
1507  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
1508  if(vDim == iDim){
1509  if(std::abs(dFluxBuffer[vIndex] - dFluxSmallRhoVMain) > epsLarge) {
1510  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
1511  << ", expected " << dFluxSmallRhoVMain << std::endl;
1512  parTestWENO = false;
1513  }
1514  } else {
1515  if(std::abs(dFluxBuffer[vIndex] - dFluxSmallRhoVAlt) > epsLarge) {
1516  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
1517  << ", expected " << dFluxSmallRhoVAlt << std::endl;
1518  parTestWENO = false;
1519  }
1520  }
1521  }
1522  if(std::abs(dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] -
1523  dFluxSmallRhoE) > epsLarge) {
1524  std::cout << "dFlux check failed: rhoE = "
1525  << dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
1526  << ", expected " << dFluxSmallRhoE << std::endl;
1527  parTestWENO = false;
1528  }
1529  } else if (iDim == jumpDim && std::abs(cloc - checkLoc2) < eps) {
1530  if(std::abs(dFluxBuffer[bufferIndex] + dFluxLargeRho) > epsLarge) {
1531  parTestWENO = false;
1532  std::cout << "dFlux check failed: rho = " << dFluxBuffer[bufferIndex]
1533  << ", expected " << -dFluxLargeRho << std::endl;
1534  }
1535  for(int vDim = 0;vDim < numDim;vDim++){
1536  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
1537  if(vDim == iDim){
1538  if(std::abs(dFluxBuffer[vIndex] + dFluxLargeRhoVMain) > epsLarge) {
1539  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
1540  << ", expected " << -dFluxLargeRhoVMain << std::endl;
1541  parTestWENO = false;
1542  }
1543  } else {
1544  if(std::abs(dFluxBuffer[vIndex] + dFluxLargeRhoVAlt) > epsLarge) {
1545  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
1546  << ", expected " << -dFluxLargeRhoVAlt << std::endl;
1547  parTestWENO = false;
1548  }
1549  }
1550  }
1551  if(std::abs(dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] +
1552  dFluxLargeRhoE) > epsLarge) {
1553  std::cout << "dFlux check failed: rhoE = "
1554  << dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
1555  << ", expected " << -dFluxLargeRhoE << std::endl;
1556  parTestWENO = false;
1557  }
1558  } else if (iDim == jumpDim && std::abs(cloc - (checkLoc2-dx)) < eps) {
1559  if(std::abs(dFluxBuffer[bufferIndex] + dFluxSmallRho) > epsLarge) {
1560  parTestWENO = false;
1561  std::cout << "dFlux check failed: rho = " << dFluxBuffer[bufferIndex]
1562  << ", expected " << -dFluxSmallRho << std::endl;
1563  }
1564  for(int vDim = 0;vDim < numDim;vDim++){
1565  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
1566  if(vDim == iDim){
1567  if(std::abs(dFluxBuffer[vIndex] + dFluxSmallRhoVMain) > epsLarge) {
1568  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
1569  << ", expected " << -dFluxSmallRhoVMain << std::endl;
1570  parTestWENO = false;
1571  }
1572  } else {
1573  if(std::abs(dFluxBuffer[vIndex] + dFluxSmallRhoVAlt) > epsLarge) {
1574  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
1575  << ", expected " << -dFluxSmallRhoVAlt << std::endl;
1576  parTestWENO = false;
1577  }
1578  }
1579  }
1580  if(std::abs(dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] +
1581  dFluxSmallRhoE) > epsLarge) {
1582  std::cout << "dFlux check failed: rhoE = "
1583  << dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
1584  << ", expected " << -dFluxSmallRhoE << std::endl;
1585  parTestWENO = false;
1586  }
1587  } else if (iDim == jumpDim && (std::abs(cloc - checkLoc1) < stenWidth*dx+eps)
1588  || (std::abs(cloc - checkLoc2) < stenWidth*dx+eps)) {
1589  // near discontinuities WENO can have larger values of "0"
1590  // (i.e. numerical dissipation)
1591  for(int eqn = 0; eqn < numDim+2; eqn++){
1592  size_t eIndex = bufferIndex + eqn*numPointsBuffer;
1593  if(std::abs(dFluxBuffer[eIndex]) > epsLarge)
1594  {
1595  parTestWENO = false;
1596  std::cout << "Failed dFluxBuffer " << dFluxBuffer[eIndex] << std::endl;
1597  }
1598  }
1599  } else {
1600  for(int eqn = 0; eqn < numDim+2; eqn++){
1601  size_t eIndex = bufferIndex + eqn*numPointsBuffer;
1602  if(std::abs(dFluxBuffer[eIndex]) > eps)
1603  {
1604  parTestWENO = false;
1605  std::cout << "Failed dFluxBuffer " << dFluxBuffer[eIndex] << std::endl;
1606  }
1607  }
1608  } // checkLoc if statement
1609 
1610  if (!parTestWENO) {
1611  std::cout << "i = " << iIndex << ", j = " << jIndex << ", k = " << kIndex
1612  << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
1613  std::cout << "checkLoc1 = " << checkLoc1 << ", checkLoc2 = " << checkLoc2 << std::endl;
1614  std::cout << "rho = " << rho[bufferIndex];
1615  for (int vDim=0; vDim<numDim; vDim++) {
1616  std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
1617  }
1618  std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
1619  }
1620  }
1621  }
1622  }
1623 
1624  parallelUnitResults.UpdateResult("WENO:ApplyWENO:TransonicJump:dFluxCheck",parTestWENO);
1625 
1626  if(!parTestWENO)
1627  return;
1628 
1629  } // Loop over dimensions
1630  } // jumpDim loop
1631  } // Loop over numDim = 2,3
1632 
1633  // --- JumpMetric ---
1634  // Tests that when states have initial discontinuity, dFluxBuffer is
1635  // computed correctly
1636  // Grid setup so metric terms will be involved
1637 
1638  for(int numDim = 2;numDim < 4;numDim++){
1639 
1640  std::cout << "TestWENO:ApplyWENO:JumpMetric:Testing WENO in " << numDim << " dimensions." << std::endl;
1641 
1642  for (int jumpDim=0; jumpDim<numDim; jumpDim++) {
1643  grid_t simGrid;
1644  operator_t sbpOperator;
1645  halo_t &simHalo(simGrid.Halo());
1646  state_t simState;
1647  state_t paramState;
1648  rhs_t rhsWENO;
1649  weno_t &coeffsWENO(rhsWENO.WENOCoefficients());
1650 
1651  bool parTestWENO = true;
1652  int numPoints = 41;
1653  double length = 20.0;
1654  std::vector<size_t> gridSizes(numDim,numPoints);
1655  double dx = length/(numPoints - 1);
1656 
1657  if(testfixtures::CreateSimulationFixtures(simGrid,simHalo,sbpOperator,coeffsWENO,
1658  gridSizes,2,testComm,&messageStream)){
1659  parTestWENO = false;
1660  }
1661 
1662  if(simGrid.GenerateCoordinates(messageStream))
1663  parTestWENO = false;
1664 
1665  std::vector<double> &gridCoords(simGrid.CoordinateData());
1666 
1667  fixtures::CommunicatorType *gridCommPtr;
1668  std::vector<int> cartNeighbors;
1669  gridCommPtr = &simGrid.Communicator();
1670  gridCommPtr->CartNeighbors(cartNeighbors);
1671 
1672  // Fill halo regions with coordinates
1673  int xyzMessageId = simHalo.CreateMessageBuffers(numDim);
1674  simHalo.PackMessageBuffers(xyzMessageId,0,numDim,&gridCoords[0]);
1675 
1676  simHalo.SendMessage(0,cartNeighbors,testComm);
1677  simHalo.ReceiveMessage(0,cartNeighbors,testComm);
1678  simHalo.UnPackMessageBuffers(0,0,numDim,&gridCoords[0]);
1679 
1680  parallelUnitResults.UpdateResult("WENO:ApplyWENO:JumpMetric:InitFixtures",parTestWENO);
1681 
1682  if(!parTestWENO)
1683  return;
1684 
1685  pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());
1686  pcpp::IndexIntervalType &partitionInterval(simGrid.PartitionInterval());
1687  std::vector<size_t> &bufferSizes(simGrid.BufferSizes());
1688 
1689  if(testfixtures::euler::SetupEulerState(simGrid,simState,paramState,0)){
1690  parTestWENO = false;
1691  }
1692 
1693  parallelUnitResults.UpdateResult("WENO:ApplyWENO:JumpMetric:SetupState",parTestWENO);
1694 
1695  if(!parTestWENO)
1696  return;
1697 
1698  std::cout << messageStream.str() << std::endl;
1699  std::cout << "Testing WENO jump in " << jumpDim << " dimension." << std::endl;
1700 
1701  size_t numPointsBuffer = simGrid.BufferSize();
1702 
1703  double t = 0.0;
1704  double gamma = 1.4;
1705  double dt = .001;
1706  double cfl = 1.0;
1707  double Re = 0.0;
1708  std::vector<double> numbers(4,0.0);
1709  int flags = USEWENO;
1710 
1711  // Add some fields that are not usual
1712  paramState.Create(numPointsBuffer,0);
1713 
1714  paramState.SetFieldBuffer("gamma",&gamma);
1715  paramState.SetFieldBuffer("inputDT",&dt);
1716  paramState.SetFieldBuffer("inputCFL",&cfl);
1717  paramState.SetFieldBuffer("refRe",&Re);
1718  paramState.SetFieldBuffer("Flags",&flags);
1719  paramState.SetFieldBuffer("Numbers",&numbers[0]);
1720 
1721  std::vector<double> rho(numPointsBuffer,-1000.0);
1722  std::vector<double> rhoV(numDim*numPointsBuffer,-1000.0);
1723  std::vector<double> rhoE(numPointsBuffer,-1000.0);
1724 
1725  // "left" state (edges of domain in direction jumpDim)
1726  double rhoL = 1.0;
1727  double rhoVL = 2.0;
1728  double rhoEL = 4.0 + 2.0*numDim;
1729 
1730  // "right" state (center of domain in direction jumpDim)
1731  double rhoR = 2.0;
1732  double rhoVR = 6.0;
1733  double rhoER = 8.0 + 9.0*numDim;
1734 
1735  // locations of jumps, using physical coordinates
1736  double loc1 = 0.25*length;
1737  double loc2 = 0.75*length;
1738 
1739  // locations of jumps, using physical coordinates
1740  double checkLoc1 = loc1;
1741  double checkLoc2 = loc2 + dx;
1742 
1743  size_t iStart = partitionBufferInterval[0].first;
1744  size_t iEnd = partitionBufferInterval[0].second;
1745  size_t jStart = partitionBufferInterval[1].first;
1746  size_t jEnd = partitionBufferInterval[1].second;
1747 
1748  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
1749  // computed correctly
1750  size_t kStart = 0;
1751  size_t kEnd = 0;
1752  if (numDim == 3) {
1753  kStart = partitionBufferInterval[2].first;
1754  kEnd = partitionBufferInterval[2].second;
1755  }
1756 
1757  for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
1758  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
1759  for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
1760  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
1761  for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
1762  size_t bufferIndex = jBufferIndex + iIndex;
1763 
1764  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
1765 
1766  if (cloc < loc1 || cloc > loc2) {
1767  rho[bufferIndex] = rhoL;
1768  for (int vDim=0; vDim<numDim; vDim++) {
1769  rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVL;
1770  }
1771  rhoE[bufferIndex] = rhoEL;
1772  } else {
1773  rho[bufferIndex] = rhoR;
1774  for (int vDim=0; vDim<numDim; vDim++) {
1775  rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVR;
1776  }
1777  rhoE[bufferIndex] = rhoER;
1778  }
1779  }
1780  }
1781  }
1782 
1783  std::vector<double> dRho(numPointsBuffer,-1000.0);
1784  std::vector<double> dRhoV(numDim*numPointsBuffer,-1000.0);
1785  std::vector<double> dRhoE(numPointsBuffer,-1000.0);
1786 
1787  std::vector<double> T(numPointsBuffer,0.0);
1788  std::vector<double> p(numPointsBuffer,0.0);
1789  std::vector<double> V(numDim*numPointsBuffer,0.0);
1790 
1791  std::vector<double> rhom1(numPointsBuffer,0.0);
1792 
1793  simState.SetFieldBuffer("simTime",&t);
1794  simState.SetFieldBuffer("rho",&rho[0]);
1795  simState.SetFieldBuffer("rhoV",&rhoV[0]);
1796  simState.SetFieldBuffer("rhoE",&rhoE[0]);
1797 
1798  simState.SetFieldBuffer("pressure",&p[0]);
1799  simState.SetFieldBuffer("temperature",&T[0]);
1800  simState.SetFieldBuffer("velocity",&V[0]);
1801  simState.SetFieldBuffer("rhom1",&rhom1[0]);
1802 
1803  state_t rhsState(simState);
1804 
1805  rhsState.SetFieldBuffer("rho",&dRho[0]);
1806  rhsState.SetFieldBuffer("rhoV",&dRhoV[0]);
1807  rhsState.SetFieldBuffer("rhoE",&dRhoE[0]);
1808 
1809  // State/Soln initialization
1810  std::vector<double> inParams;
1811  std::vector<int> inFlags;
1812 
1813  if(rhsWENO.Initialize(simGrid,simState,paramState,sbpOperator)){
1814  parTestWENO = false;
1815  }
1816  std::vector<int> numThreadsDim;
1817  if(simGrid.SetupThreads(numThreadsDim)){
1818  std::cout << "Setting up threads for grid object failed." << std::endl;
1819  return;
1820  }
1821  rhsWENO.InitThreadIntervals();
1822 
1823  parallelUnitResults.UpdateResult("WENO:ApplyWENO:JumpMetric:InitRHS",parTestWENO);
1824 
1825  if(!parTestWENO)
1826  return;
1827 
1828  if(rhsWENO.SetupRHS(simGrid,simState,rhsState)){
1829  parTestWENO = false;
1830  }
1831 
1832  parallelUnitResults.UpdateResult("WENO:ApplyWENO:JumpMetric:SetupRHS",parTestWENO);
1833 
1834  if(!parTestWENO)
1835  return;
1836 
1837  //std::cout << "RHO BEFORE: " << std::endl;
1838  //pcpp::io::DumpContents(std::cout,rho," ");
1839  //std::cout << std::endl;
1840 
1841  // Exchange the state
1842  rhsWENO.ExchangeState(0);
1843 
1844  //std::cout << "RHO AFTER: " << std::endl;
1845  //pcpp::io::DumpContents(std::cout,rho," ");
1846  //std::cout << std::endl;
1847 
1848  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
1849  // computed correctly
1850  size_t kFinal = 1;
1851  if (numDim == 3) {
1852  kFinal = bufferSizes[2];
1853  }
1854 
1855  for (size_t kIndex = 0; kIndex < kFinal && parTestWENO; kIndex++) {
1856  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
1857  for (size_t jIndex = 0; jIndex < bufferSizes[1] && parTestWENO; jIndex++) {
1858  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
1859  for (size_t iIndex = 0; iIndex < bufferSizes[0] && parTestWENO; iIndex++) {
1860  size_t bufferIndex = jBufferIndex + iIndex;
1861 
1862  int c = 0;
1863  if (iIndex < iStart || iIndex > iEnd) c++;
1864  if (jIndex < jStart || jIndex > jEnd) c++;
1865  if (numDim == 3 && (kIndex < kStart || kIndex > kEnd)) c++;
1866  if (c != 1) continue; // then we are not in a halo region
1867 
1868  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
1869 
1870  if (cloc < loc1 || cloc > loc2) {
1871  if (std::abs(rho[bufferIndex] - rhoL) > eps) {
1872  parTestWENO = false;
1873  std::cout << "Failed state exchange, rho " << rho[bufferIndex] << std::endl;
1874  }
1875  for (int vDim=0; vDim<numDim; vDim++) {
1876  if (std::abs(rhoV[bufferIndex + vDim*numPointsBuffer] - rhoVL) > eps) {
1877  parTestWENO = false;
1878  std::cout << "Failed state exchange, rhoV" << vDim << " "
1879  << rhoV[bufferIndex + vDim*numPointsBuffer] << std::endl;
1880  }
1881  }
1882  if (std::abs(rhoE[bufferIndex] - rhoEL) > eps) {
1883  parTestWENO = false;
1884  std::cout << "Failed state exchange, rhoE " << rhoE[bufferIndex] << std::endl;
1885  }
1886  } else {
1887  if (std::abs(rho[bufferIndex] - rhoR) > eps) {
1888  parTestWENO = false;
1889  std::cout << "Failed state exchange, rho " << rho[bufferIndex] << std::endl;
1890  }
1891  for (int vDim=0; vDim<numDim; vDim++) {
1892  if (std::abs(rhoV[bufferIndex + vDim*numPointsBuffer] - rhoVR) > eps) {
1893  parTestWENO = false;
1894  std::cout << "Failed state exchange, rhoV" << vDim << " "
1895  << rhoV[bufferIndex + vDim*numPointsBuffer] << std::endl;
1896  }
1897  }
1898  if (std::abs(rhoE[bufferIndex] - rhoER) > eps) {
1899  parTestWENO = false;
1900  std::cout << "Failed state exchange, rhoE " << rhoE[bufferIndex] << std::endl;
1901  }
1902  }
1903 
1904  if (!parTestWENO) {
1905  std::cout << "iIndex = " << iIndex << ", jIndex = " << jIndex
1906  << ", kIndex = " << kIndex << std::endl;
1907  std::cout << "iStart = " << iStart << ", jStart = " << jStart
1908  << ", kStart = " << kStart << std::endl;
1909  std::cout << "iEnd = " << iEnd << ", jEnd = " << jEnd
1910  << ", kEnd = " << kEnd << std::endl;
1911  }
1912  }
1913  }
1914  }
1915 
1916  parallelUnitResults.UpdateResult("WENO:ApplyWENO:JumpMetric:ExchangeState",parTestWENO);
1917 
1918  if(!parTestWENO)
1919  return;
1920 
1921  double *inviscidFluxBuffer(rhsWENO.InviscidFluxBuffer());
1922  double *dFluxBuffer(rhsWENO.DFluxBuffer());
1923 
1924  for(int iDim = 0;iDim < numDim;iDim++){
1925 
1926  // Running WENO in dimension (iDim)
1927  std::cout << "Checking WENO in dimension(" << iDim << ")"
1928  << std::endl;
1929 
1930  // initialize dFluxBuffer so test can't pass because data was initialized to zero
1931  for (size_t i=0; i<numPointsBuffer; i++) {
1932  for (int j=0; j<numDim+2; j++) {
1933  size_t bufferIndex = j*numPointsBuffer + i;
1934  dFluxBuffer[bufferIndex] = -1000;
1935  }
1936  }
1937 
1938  rhsWENO.ComputeDV(0);
1939  rhsWENO.InviscidFlux(iDim,0);
1940  rhsWENO.ExchangeInviscidFlux(0);
1941 
1942  //if(numDim == 2){
1943  // std::cout << "Velocity = " << std::endl;
1944  // for(int iVel = 0;iVel < numPointsBuffer;iVel++){
1945  // std::cout << "Vel(" << iVel << ") = ("
1946  // << V[iVel] << "," << V[iVel+numPointsBuffer]
1947  // << ")" << std::endl;
1948  // std::cout << "PRESSURE: " << p[iVel] << std::endl;
1949  // std::cout << "FLUX_RHO: " << inviscidFluxBuffer[iVel] << std::endl;
1950  // std::cout << "FLUX_RHOV1: " << inviscidFluxBuffer[iVel+numPointsBuffer] << std::endl;
1951  // std::cout << "FLUX_RHOV2: " << inviscidFluxBuffer[iVel+2*numPointsBuffer] << std::endl;
1952  // std::cout << "FLUX_E: " << inviscidFluxBuffer[iVel+3*numPointsBuffer] << std::endl;
1953  // }
1954  //}
1955 
1956  double metric = 1.0/(2*(numDim-1)); // actually 1/(2**(numDim-1)), but this is equivalent for numDim = 2,3
1957 
1958  // "left" state (edges of domain in direction jumpDim)
1959  double fluxRhoL = 2.0*metric;
1960  double fluxRhoVMainL = 5.6*metric;
1961  double fluxRhoVAltL = 4.0*metric;
1962  double fluxRhoEL = (11.2 + 4.0*numDim)*metric;
1963 
1964  // "right" state (center of domain in direction jumpDim)
1965  double fluxRhoR = 6.0*metric;
1966  double fluxRhoVMainR = 21.2*metric;
1967  double fluxRhoVAltR = 18.0*metric;
1968  double fluxRhoER = (33.6 + 27.0*numDim)*metric;
1969 
1970  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
1971  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
1972  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
1973  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
1974  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
1975  size_t bufferIndex = jBufferIndex + iIndex;
1976 
1977  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
1978 
1979  if (cloc < loc1 || cloc > loc2) {
1980  if(std::abs(inviscidFluxBuffer[bufferIndex] - fluxRhoL) > eps) {
1981  parTestWENO = false;
1982  std::cout << "Flux check failed: rho = " << inviscidFluxBuffer[bufferIndex]
1983  << ", expected " << fluxRhoL << std::endl;
1984  }
1985  for(int vDim = 0;vDim < numDim;vDim++){
1986  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
1987  if(vDim == iDim){
1988  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVMainL) > eps) {
1989  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
1990  << ", expected " << fluxRhoVMainL << std::endl;
1991  parTestWENO = false;
1992  }
1993  } else {
1994  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVAltL) > eps) {
1995  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
1996  << ", expected " << fluxRhoVAltL << std::endl;
1997  parTestWENO = false;
1998  }
1999  }
2000  }
2001  if(std::abs(inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] -
2002  fluxRhoEL) > eps) {
2003  std::cout << "Flux check failed: rhoE = "
2004  << inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
2005  << ", expected " << fluxRhoEL << std::endl;
2006  parTestWENO = false;
2007  }
2008  } else {
2009  if(std::abs(inviscidFluxBuffer[bufferIndex] - fluxRhoR) > eps) {
2010  parTestWENO = false;
2011  std::cout << "Flux check failed: rho = " << inviscidFluxBuffer[bufferIndex]
2012  << ", expected " << fluxRhoR << std::endl;
2013  }
2014  for(int vDim = 0;vDim < numDim;vDim++){
2015  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
2016  if(vDim == iDim){
2017  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVMainR) > eps) {
2018  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
2019  << ", expected " << fluxRhoVMainR << std::endl;
2020  parTestWENO = false;
2021  }
2022  } else {
2023  if(std::abs(inviscidFluxBuffer[vIndex] - fluxRhoVAltR) > eps) {
2024  std::cout << "Flux check failed: rhoV = " << inviscidFluxBuffer[vIndex]
2025  << ", expected " << fluxRhoVAltR << std::endl;
2026  parTestWENO = false;
2027  }
2028  }
2029  }
2030  if(std::abs(inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] -
2031  fluxRhoER) > eps) {
2032  std::cout << "Flux check failed: rhoE = "
2033  << inviscidFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
2034  << ", expected " << fluxRhoER << std::endl;
2035  parTestWENO = false;
2036  }
2037  } // cloc if statement
2038 
2039  if (!parTestWENO) {
2040  std::cout << "i = " << iIndex << ", j = " << jIndex << ", k = " << kIndex
2041  << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
2042  std::cout << "rho = " << rho[bufferIndex];
2043  for (int vDim=0; vDim<numDim; vDim++) {
2044  std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
2045  }
2046  std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
2047  }
2048  }
2049  }
2050  }
2051 
2052  parallelUnitResults.UpdateResult("WENO:ApplyWENO:JumpMetric:FluxCheck",parTestWENO);
2053 
2054  if(!parTestWENO)
2055  return;
2056 
2057  int statusWENO = rhsWENO.ApplyWENO(iDim,0);
2058  if(statusWENO) {
2059  parTestWENO = false;
2060  std::cout << "ApplyWENO returned error code " << statusWENO << std::endl;
2061  }
2062 
2063  parallelUnitResults.UpdateResult("WENO:ApplyWENO:JumpMetric:ErrorApplyWENO",parTestWENO);
2064 
2065  if(!parTestWENO)
2066  return;
2067 
2068  // setup expected values
2069  double dFluxRho = fluxRhoR - fluxRhoL;
2070  double dFluxRhoVMain = fluxRhoVMainR - fluxRhoVMainL;
2071  double dFluxRhoVAlt = fluxRhoVAltR - fluxRhoVAltL;
2072  double dFluxRhoE = fluxRhoER - fluxRhoEL;
2073 
2074  // Check dFluxBuffer -- should be zero everywhere except at initial discontinuity
2075  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
2076  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
2077  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
2078  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
2079  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
2080  size_t bufferIndex = jBufferIndex + iIndex;
2081 
2082  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
2083 
2084  if (iDim == jumpDim && std::abs(cloc - checkLoc1) < eps) {
2085  if(std::abs(dFluxBuffer[bufferIndex] - dFluxRho) > epsLarge) {
2086  parTestWENO = false;
2087  std::cout << "dFlux check failed: rho = " << dFluxBuffer[bufferIndex]
2088  << ", expected " << dFluxRho << std::endl;
2089  }
2090  for(int vDim = 0;vDim < numDim;vDim++){
2091  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
2092  if(vDim == iDim){
2093  if(std::abs(dFluxBuffer[vIndex] - dFluxRhoVMain) > epsLarge) {
2094  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
2095  << ", expected " << dFluxRhoVMain << std::endl;
2096  parTestWENO = false;
2097  }
2098  } else {
2099  if(std::abs(dFluxBuffer[vIndex] - dFluxRhoVAlt) > epsLarge) {
2100  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
2101  << ", expected " << dFluxRhoVAlt << std::endl;
2102  parTestWENO = false;
2103  }
2104  }
2105  }
2106  if(std::abs(dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] -
2107  dFluxRhoE) > epsLarge) {
2108  std::cout << "dFlux check failed: rhoE = "
2109  << dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
2110  << ", expected " << dFluxRhoE << std::endl;
2111  parTestWENO = false;
2112  }
2113  } else if (iDim == jumpDim && std::abs(cloc - checkLoc2) < eps) {
2114  if(std::abs(dFluxBuffer[bufferIndex] + dFluxRho) > epsLarge) {
2115  parTestWENO = false;
2116  std::cout << "dFlux check failed: rho = " << dFluxBuffer[bufferIndex]
2117  << ", expected " << -dFluxRho << std::endl;
2118  }
2119  for(int vDim = 0;vDim < numDim;vDim++){
2120  size_t vIndex = bufferIndex + (vDim+1)*numPointsBuffer;
2121  if(vDim == iDim){
2122  if(std::abs(dFluxBuffer[vIndex] + dFluxRhoVMain) > epsLarge) {
2123  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
2124  << ", expected " << -dFluxRhoVMain << std::endl;
2125  parTestWENO = false;
2126  }
2127  } else {
2128  if(std::abs(dFluxBuffer[vIndex] + dFluxRhoVAlt) > epsLarge) {
2129  std::cout << "dFlux check failed: rhoV = " << dFluxBuffer[vIndex]
2130  << ", expected " << -dFluxRhoVAlt << std::endl;
2131  parTestWENO = false;
2132  }
2133  }
2134  }
2135  if(std::abs(dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer] +
2136  dFluxRhoE) > epsLarge) {
2137  std::cout << "dFlux check failed: rhoE = "
2138  << dFluxBuffer[bufferIndex + (numDim+1)*numPointsBuffer]
2139  << ", expected " << -dFluxRhoE << std::endl;
2140  parTestWENO = false;
2141  }
2142  } else if (iDim == jumpDim && (std::abs(cloc - checkLoc1) < stenWidth*dx+eps)
2143  || (std::abs(cloc - checkLoc2) < stenWidth*dx+eps)) {
2144  // near discontinuities WENO can have larger values of "0"
2145  // (i.e. numerical dissipation)
2146  for(int eqn = 0; eqn < numDim+2; eqn++){
2147  size_t eIndex = bufferIndex + eqn*numPointsBuffer;
2148  if(std::abs(dFluxBuffer[eIndex]) > epsLarge)
2149  {
2150  parTestWENO = false;
2151  std::cout << "Failed dFluxBuffer " << dFluxBuffer[eIndex] << std::endl;
2152  }
2153  }
2154  } else {
2155  for(int eqn = 0; eqn < numDim+2; eqn++){
2156  size_t eIndex = bufferIndex + eqn*numPointsBuffer;
2157  if(std::abs(dFluxBuffer[eIndex]) > eps)
2158  {
2159  parTestWENO = false;
2160  std::cout << "Failed dFluxBuffer " << dFluxBuffer[eIndex] << std::endl;
2161  }
2162  }
2163  } // checkLoc if statement
2164 
2165  if (!parTestWENO) {
2166  std::cout << "i = " << iIndex << ", j = " << jIndex << ", k = " << kIndex
2167  << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
2168  std::cout << "checkLoc1 = " << checkLoc1 << ", checkLoc2 = " << checkLoc2 << std::endl;
2169  std::cout << "rho = " << rho[bufferIndex];
2170  for (int vDim=0; vDim<numDim; vDim++) {
2171  std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
2172  }
2173  std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
2174  }
2175  }
2176  }
2177  }
2178 
2179  parallelUnitResults.UpdateResult("WENO:ApplyWENO:JumpMetric:dFluxCheck",parTestWENO);
2180 
2181  if(!parTestWENO)
2182  return;
2183 
2184  } // Loop over dimensions
2185  } // jumpDim loop
2186  } // Loop over numDim = 2,3
2187 };
2188 
2189 void TestWENO_RHS(ix::test::results &parallelUnitResults,pcpp::CommunicatorType &testComm)
2190 {
2191  // Test whole RHS with WENO
2192 
2193  double eps = 1e-15;
2194  double epsLarge = 1e-2; // right near jumps, WENO reports values that are not quite zero
2195  int stenWidth = 3; // max length of one WENO stencil
2196 
2197  typedef plascom2::operators::sbp::base operator_t;
2198  typedef simulation::grid::halo halo_t;
2199  typedef simulation::state::base state_t;
2200  typedef pcpp::ParallelGlobalType global_t;
2201  typedef pcpp::CommunicatorType comm_t;
2202  typedef pcpp::IndexIntervalType interval_t;
2203 
2204  typedef simulation::grid::parallel_blockstructured grid_t;
2205  typedef simulation::domain::base<grid_t,state_t> domain_t;
2206  typedef euler::rhs<grid_t,state_t,operator_t> rhs_t;
2207  typedef WENO::CoeffsWENO weno_t;
2208 
2209  std::ostringstream messageStream;
2210 
2211  // --- Equal ---
2212  // Tests that when states are initially constant, RHS is 0
2213  for(int numDim = 2;numDim < 4;numDim++){
2214 
2215  grid_t simGrid;
2216  operator_t sbpOperator;
2217  halo_t &simHalo(simGrid.Halo());
2218  state_t simState;
2219  state_t paramState;
2220  rhs_t rhsWENO;
2221  weno_t &coeffsWENO(rhsWENO.WENOCoefficients());
2222 
2223  bool parTestWENO = true;
2224  std::vector<size_t> gridSizes(numDim,21);
2225 
2226  if(testfixtures::CreateSimulationFixtures(simGrid,simHalo,sbpOperator,coeffsWENO,
2227  gridSizes,2,testComm,&messageStream)){
2228  parTestWENO = false;
2229  }
2230 
2231  parallelUnitResults.UpdateResult("WENO:RHS:Equal:InitFixtures",parTestWENO);
2232 
2233  if(!parTestWENO)
2234  return;
2235 
2236  pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());
2237  pcpp::IndexIntervalType &partitionInterval(simGrid.PartitionInterval());
2238  std::vector<size_t> &bufferSizes(simGrid.BufferSizes());
2239 
2240  if(testfixtures::euler::SetupEulerState(simGrid,simState,paramState,0)){
2241  parTestWENO = false;
2242  }
2243 
2244  parallelUnitResults.UpdateResult("WENO:RHS:Equal:SetupState",parTestWENO);
2245 
2246  if(!parTestWENO)
2247  return;
2248 
2249  std::cout << messageStream.str() << std::endl;
2250  std::cout << "TestWENO_RHS:Equal:Testing WENO in " << numDim << " dimensions." << std::endl;
2251  size_t numPointsBuffer = simGrid.BufferSize();
2252 
2253  double t = 0.0;
2254  double gamma = 1.4;
2255  double dt = .001;
2256  double cfl = 1.0;
2257  double Re = 0.0;
2258  std::vector<double> numbers(4,0.0);
2259  int flags = USEWENO;
2260 
2261  // Add some fields that are not usual
2262  paramState.Create(numPointsBuffer,0);
2263 
2264  paramState.SetFieldBuffer("gamma",&gamma);
2265  paramState.SetFieldBuffer("inputDT",&dt);
2266  paramState.SetFieldBuffer("inputCFL",&cfl);
2267  paramState.SetFieldBuffer("refRe",&Re);
2268  paramState.SetFieldBuffer("Flags",&flags);
2269  paramState.SetFieldBuffer("Numbers",&numbers[0]);
2270 
2271  std::vector<double> rho(numPointsBuffer,-1000.0);
2272  std::vector<double> rhoV(numDim*numPointsBuffer,-1000.0);
2273  std::vector<double> rhoE(numPointsBuffer,-1000.0);
2274 
2275  size_t iStart = partitionBufferInterval[0].first;
2276  size_t iEnd = partitionBufferInterval[0].second;
2277  size_t jStart = partitionBufferInterval[1].first;
2278  size_t jEnd = partitionBufferInterval[1].second;
2279 
2280  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
2281  // computed correctly
2282  size_t kStart = 0;
2283  size_t kEnd = 0;
2284  if (numDim == 3) {
2285  kStart = partitionBufferInterval[2].first;
2286  kEnd = partitionBufferInterval[2].second;
2287  }
2288 
2289  for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
2290  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
2291  for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
2292  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
2293  for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
2294  size_t bufferIndex = jBufferIndex + iIndex;
2295  rho[bufferIndex] = 1.0;
2296  for (int vDim=0; vDim<numDim; vDim++) {
2297  rhoV[bufferIndex + vDim*numPointsBuffer] = 1.0;
2298  }
2299  rhoE[bufferIndex] = 4.0 + 0.5*numDim;
2300  }
2301  }
2302  }
2303 
2304  std::vector<double> dRho(numPointsBuffer,-1000.0);
2305  std::vector<double> dRhoV(numDim*numPointsBuffer,-1000.0);
2306  std::vector<double> dRhoE(numPointsBuffer,-1000.0);
2307 
2308  std::vector<double> T(numPointsBuffer,0.0);
2309  std::vector<double> p(numPointsBuffer,0.0);
2310  std::vector<double> V(numDim*numPointsBuffer,0.0);
2311 
2312  std::vector<double> rhom1(numPointsBuffer,0.0);
2313 
2314 
2315  simState.SetFieldBuffer("simTime",&t);
2316  simState.SetFieldBuffer("rho",&rho[0]);
2317  simState.SetFieldBuffer("rhoV",&rhoV[0]);
2318  simState.SetFieldBuffer("rhoE",&rhoE[0]);
2319 
2320  simState.SetFieldBuffer("pressure",&p[0]);
2321  simState.SetFieldBuffer("temperature",&T[0]);
2322  simState.SetFieldBuffer("velocity",&V[0]);
2323  simState.SetFieldBuffer("rhom1",&rhom1[0]);
2324 
2325  state_t rhsState(simState);
2326 
2327  rhsState.SetFieldBuffer("rho",&dRho[0]);
2328  rhsState.SetFieldBuffer("rhoV",&dRhoV[0]);
2329  rhsState.SetFieldBuffer("rhoE",&dRhoE[0]);
2330 
2331  // State/Soln initialization
2332  std::vector<double> inParams;
2333  std::vector<int> inFlags;
2334 
2335  // if(InitializeAcousticPulse(simGrid,simState,
2336  // paramState,inParams,
2337  // inFlags,0)){
2338  // // then it failed.
2339  // }
2340 
2341 
2342  if(rhsWENO.Initialize(simGrid,simState,paramState,sbpOperator)){
2343  parTestWENO = false;
2344  }
2345  std::vector<int> numThreadsDim;
2346  if(simGrid.SetupThreads(numThreadsDim)){
2347  std::cout << "Setting up threads for grid object failed." << std::endl;
2348  return;
2349  }
2350  rhsWENO.InitThreadIntervals();
2351 
2352  parallelUnitResults.UpdateResult("WENO:RHS:Equal:InitRHS",parTestWENO);
2353 
2354  if(!parTestWENO)
2355  return;
2356 
2357  if(rhsWENO.SetupRHS(simGrid,simState,rhsState)){
2358  parTestWENO = false;
2359  }
2360 
2361  parallelUnitResults.UpdateResult("WENO:RHS:Equal:SetupRHS",parTestWENO);
2362 
2363  if(!parTestWENO)
2364  return;
2365 
2366  int status = rhsWENO.RHS(0);
2367  if(status){
2368  parTestWENO = false;
2369  std::cout << "RHS returned error code " << status << std::endl;
2370  }
2371 
2372  parallelUnitResults.UpdateResult("WENO:RHS:Equal:ErrorRHS",parTestWENO);
2373 
2374  if(!parTestWENO)
2375  return;
2376 
2377  std::vector<double* > &rhsBuffers(rhsWENO.RHSBuffers());
2378 
2379  for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
2380  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
2381  for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
2382  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
2383  for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
2384  size_t bufferIndex = jBufferIndex + iIndex;
2385  for (int eqn=0; eqn<numDim+2; eqn++) {
2386  if (rhsBuffers[eqn][bufferIndex] > eps) {
2387  parTestWENO = false;
2388  std::cout << "Failed rhsBuffers[" << eqn << "] "
2389  << rhsBuffers[eqn][bufferIndex]
2390  << ", i = " << iIndex << ", j = " << jIndex
2391  << ", k = " << kIndex
2392  << std::endl;
2393  }
2394  }
2395  }
2396  }
2397  }
2398 
2399  parallelUnitResults.UpdateResult("WENO:RHS:Equal:CheckRHS",parTestWENO);
2400 
2401  if(!parTestWENO)
2402  return;
2403 
2404  } // Loop over numDim = 2,3
2405 
2406  // --- Jump ---
2407  // Tests that when states have initial discontinuity, RHS is
2408  // computed correctly
2409 
2410  for(int numDim = 2;numDim < 4;numDim++){
2411 
2412  std::cout << "TestWENO:RHS:Jump:Testing WENO in " << numDim << " dimensions." << std::endl;
2413 
2414  for (int jumpDim=0; jumpDim<numDim; jumpDim++) {
2415  grid_t simGrid;
2416  operator_t sbpOperator;
2417  halo_t &simHalo(simGrid.Halo());
2418  state_t simState;
2419  state_t paramState;
2420  rhs_t rhsWENO;
2421  weno_t &coeffsWENO(rhsWENO.WENOCoefficients());
2422 
2423  bool parTestWENO = true;
2424  int numPoints = 21;
2425  double length = 20.0;
2426  std::vector<size_t> gridSizes(numDim,numPoints);
2427  double dx = length/(numPoints - 1);
2428 
2429  if(testfixtures::CreateSimulationFixtures(simGrid,simHalo,sbpOperator,coeffsWENO,
2430  gridSizes,2,testComm,&messageStream)){
2431  parTestWENO = false;
2432  }
2433 
2434  if(simGrid.GenerateCoordinates(messageStream))
2435  parTestWENO = false;
2436 
2437  std::vector<double> &gridCoords(simGrid.CoordinateData());
2438 
2439  fixtures::CommunicatorType *gridCommPtr;
2440  std::vector<int> cartNeighbors;
2441  gridCommPtr = &simGrid.Communicator();
2442  gridCommPtr->CartNeighbors(cartNeighbors);
2443 
2444  // Fill halo regions with coordinates
2445  int xyzMessageId = simHalo.CreateMessageBuffers(numDim);
2446  simHalo.PackMessageBuffers(xyzMessageId,0,numDim,&gridCoords[0]);
2447 
2448  simHalo.SendMessage(0,cartNeighbors,testComm);
2449  simHalo.ReceiveMessage(0,cartNeighbors,testComm);
2450  simHalo.UnPackMessageBuffers(0,0,numDim,&gridCoords[0]);
2451 
2452  parallelUnitResults.UpdateResult("WENO:RHS:Jump:InitFixtures",parTestWENO);
2453 
2454  if(!parTestWENO)
2455  return;
2456 
2457  pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());
2458  pcpp::IndexIntervalType &partitionInterval(simGrid.PartitionInterval());
2459  std::vector<size_t> &bufferSizes(simGrid.BufferSizes());
2460 
2461  if(testfixtures::euler::SetupEulerState(simGrid,simState,paramState,0)){
2462  parTestWENO = false;
2463  }
2464 
2465  parallelUnitResults.UpdateResult("WENO:RHS:Jump:SetupState",parTestWENO);
2466 
2467  if(!parTestWENO)
2468  return;
2469 
2470  std::cout << messageStream.str() << std::endl;
2471  std::cout << "Testing WENO jump in " << jumpDim << " dimension." << std::endl;
2472 
2473  size_t numPointsBuffer = simGrid.BufferSize();
2474 
2475  double t = 0.0;
2476  double gamma = 1.4;
2477  double dt = .001;
2478  double cfl = 1.0;
2479  double Re = 0.0;
2480  std::vector<double> numbers(4,0.0);
2481  int flags = USEWENO;
2482 
2483  // Add some fields that are not usual
2484  paramState.Create(numPointsBuffer,0);
2485 
2486  paramState.SetFieldBuffer("gamma",&gamma);
2487  paramState.SetFieldBuffer("inputDT",&dt);
2488  paramState.SetFieldBuffer("inputCFL",&cfl);
2489  paramState.SetFieldBuffer("refRe",&Re);
2490  paramState.SetFieldBuffer("Flags",&flags);
2491  paramState.SetFieldBuffer("Numbers",&numbers[0]);
2492 
2493  std::vector<double> rho(numPointsBuffer,-1000.0);
2494  std::vector<double> rhoV(numDim*numPointsBuffer,-1000.0);
2495  std::vector<double> rhoE(numPointsBuffer,-1000.0);
2496 
2497  // "left" state (edges of domain in direction jumpDim)
2498  double rhoL = 1.0;
2499  double rhoVL = 2.0;
2500  double rhoEL = 4.0 + 2.0*numDim;
2501 
2502  // "right" state (center of domain in direction jumpDim)
2503  double rhoR = 2.0;
2504  double rhoVR = 6.0;
2505  double rhoER = 8.0 + 9.0*numDim;
2506 
2507  // locations of jumps, using physical coordinates
2508  double loc1 = 0.25*length;
2509  double loc2 = 0.75*length;
2510 
2511  // locations of jumps, using physical coordinates
2512  double checkLoc1 = loc1;
2513  double checkLoc2 = loc2 + dx;
2514 
2515  size_t iStart = partitionBufferInterval[0].first;
2516  size_t iEnd = partitionBufferInterval[0].second;
2517  size_t jStart = partitionBufferInterval[1].first;
2518  size_t jEnd = partitionBufferInterval[1].second;
2519 
2520  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
2521  // computed correctly
2522  size_t kStart = 0;
2523  size_t kEnd = 0;
2524  if (numDim == 3) {
2525  kStart = partitionBufferInterval[2].first;
2526  kEnd = partitionBufferInterval[2].second;
2527  }
2528 
2529  for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
2530  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
2531  for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
2532  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
2533  for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
2534  size_t bufferIndex = jBufferIndex + iIndex;
2535 
2536  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
2537 
2538  if (cloc < loc1 || cloc > loc2) {
2539  rho[bufferIndex] = rhoL;
2540  for (int vDim=0; vDim<numDim; vDim++) {
2541  rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVL;
2542  }
2543  rhoE[bufferIndex] = rhoEL;
2544  } else {
2545  rho[bufferIndex] = rhoR;
2546  for (int vDim=0; vDim<numDim; vDim++) {
2547  rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVR;
2548  }
2549  rhoE[bufferIndex] = rhoER;
2550  }
2551  }
2552  }
2553  }
2554 
2555  std::vector<double> dRho(numPointsBuffer,-1000.0);
2556  std::vector<double> dRhoV(numDim*numPointsBuffer,-1000.0);
2557  std::vector<double> dRhoE(numPointsBuffer,-1000.0);
2558 
2559  std::vector<double> T(numPointsBuffer,0.0);
2560  std::vector<double> p(numPointsBuffer,0.0);
2561  std::vector<double> V(numDim*numPointsBuffer,0.0);
2562 
2563  std::vector<double> rhom1(numPointsBuffer,0.0);
2564 
2565  simState.SetFieldBuffer("simTime",&t);
2566  simState.SetFieldBuffer("rho",&rho[0]);
2567  simState.SetFieldBuffer("rhoV",&rhoV[0]);
2568  simState.SetFieldBuffer("rhoE",&rhoE[0]);
2569 
2570  simState.SetFieldBuffer("pressure",&p[0]);
2571  simState.SetFieldBuffer("temperature",&T[0]);
2572  simState.SetFieldBuffer("velocity",&V[0]);
2573  simState.SetFieldBuffer("rhom1",&rhom1[0]);
2574 
2575  state_t rhsState(simState);
2576 
2577  rhsState.SetFieldBuffer("rho",&dRho[0]);
2578  rhsState.SetFieldBuffer("rhoV",&dRhoV[0]);
2579  rhsState.SetFieldBuffer("rhoE",&dRhoE[0]);
2580 
2581  // State/Soln initialization
2582  std::vector<double> inParams;
2583  std::vector<int> inFlags;
2584 
2585  if(rhsWENO.Initialize(simGrid,simState,paramState,sbpOperator)){
2586  parTestWENO = false;
2587  }
2588  std::vector<int> numThreadsDim;
2589  if(simGrid.SetupThreads(numThreadsDim)){
2590  std::cout << "Setting up threads for grid object failed." << std::endl;
2591  return;
2592  }
2593  rhsWENO.InitThreadIntervals();
2594 
2595  parallelUnitResults.UpdateResult("WENO:RHS:Jump:InitRHS",parTestWENO);
2596 
2597  if(!parTestWENO)
2598  return;
2599 
2600  if(rhsWENO.SetupRHS(simGrid,simState,rhsState)){
2601  parTestWENO = false;
2602  }
2603 
2604  parallelUnitResults.UpdateResult("WENO:RHS:Jump:SetupRHS",parTestWENO);
2605 
2606  if(!parTestWENO)
2607  return;
2608 
2609  int statusWENO = rhsWENO.RHS(0);
2610  if(statusWENO) {
2611  parTestWENO = false;
2612  std::cout << "RHS returned error code " << statusWENO << std::endl;
2613  }
2614 
2615  parallelUnitResults.UpdateResult("WENO:RHS:Jump:ErrorRHS",parTestWENO);
2616 
2617  if(!parTestWENO)
2618  return;
2619 
2620  std::vector<double* > &rhsBuffers(rhsWENO.RHSBuffers());
2621 
2622  // -- setup expected values
2623 
2624  // "left" state (edges of domain in direction jumpDim)
2625  double fluxRhoL = 2.0;
2626  double fluxRhoVMainL = 5.6;
2627  double fluxRhoVAltL = 4.0;
2628  double fluxRhoEL = 11.2 + 4.0*numDim;
2629 
2630  // "right" state (center of domain in direction jumpDim)
2631  double fluxRhoR = 6.0;
2632  double fluxRhoVMainR = 21.2;
2633  double fluxRhoVAltR = 18.0;
2634  double fluxRhoER = 33.6 + 27.0*numDim;
2635 
2636  double dFluxRho = (fluxRhoR - fluxRhoL)/dx;
2637  double dFluxRhoVMain = (fluxRhoVMainR - fluxRhoVMainL)/dx;
2638  double dFluxRhoVAlt = (fluxRhoVAltR - fluxRhoVAltL)/dx;
2639  double dFluxRhoE = (fluxRhoER - fluxRhoEL)/dx;
2640 
2641 // // Print data along a pencil for debugging purposes
2642 // if (jumpDim == 0) {
2643 // size_t kBufferIndex = kStart*bufferSizes[0]*bufferSizes[1];
2644 // size_t jBufferIndex = jStart*bufferSizes[0] + kBufferIndex;
2645 // for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
2646 // size_t bufferIndex = jBufferIndex + iIndex;
2647 //
2648 // double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
2649 //
2650 // //std::cout << "i = " << iIndex << ", j = " << jStart << ", k = " << kStart
2651 // // << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
2652 // //std::cout << "checkLoc1 = " << checkLoc1 << ", checkLoc2 = " << checkLoc2 << std::endl;
2653 // //std::cout << "rho = " << rho[bufferIndex];
2654 // //for (int vDim=0; vDim<numDim; vDim++) {
2655 // // std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
2656 // //}
2657 // //std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
2658 //
2659 // std::cout << "cloc = " << cloc << ", ";
2660 //
2661 // std::cout << "rhoRHS = " << rhsBuffers[0][bufferIndex];
2662 // for (int vDim=0; vDim<numDim; vDim++) {
2663 // std::cout << ", rhoVRHS[" << vDim << "] = " << rhsBuffers[vDim+1][bufferIndex];
2664 // }
2665 // std::cout << ", rhoERHS = " << rhsBuffers[numDim+1][bufferIndex] << std::endl;
2666 // }
2667 // }
2668 
2669  // Check RHS -- should be zero everywhere except at initial discontinuity
2670  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
2671  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
2672  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
2673  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
2674  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
2675  size_t bufferIndex = jBufferIndex + iIndex;
2676 
2677  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
2678 
2679  if (std::abs(cloc - checkLoc1) < eps) {
2680  if(std::abs(rhsBuffers[0][bufferIndex] + dFluxRho) > epsLarge) {
2681  parTestWENO = false;
2682  std::cout << "RHS check failed: rho = " << rhsBuffers[0][bufferIndex]
2683  << ", expected " << -dFluxRho << std::endl;
2684  }
2685  for(int vDim = 0;vDim < numDim;vDim++){
2686  if(vDim == jumpDim){
2687  if(std::abs(rhsBuffers[vDim+1][bufferIndex] + dFluxRhoVMain) > epsLarge) {
2688  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
2689  << ", expected " << -dFluxRhoVMain << std::endl;
2690  parTestWENO = false;
2691  }
2692  } else {
2693  if(std::abs(rhsBuffers[vDim+1][bufferIndex] + dFluxRhoVAlt) > epsLarge) {
2694  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
2695  << ", expected " << -dFluxRhoVAlt << std::endl;
2696  parTestWENO = false;
2697  }
2698  }
2699  }
2700  if(std::abs(rhsBuffers[numDim+1][bufferIndex] + dFluxRhoE) > epsLarge) {
2701  std::cout << "RHS check failed: rhoE = "
2702  << rhsBuffers[numDim+1][bufferIndex]
2703  << ", expected " << -dFluxRhoE << std::endl;
2704  parTestWENO = false;
2705  }
2706  } else if (std::abs(cloc - checkLoc2) < eps) {
2707  if(std::abs(rhsBuffers[0][bufferIndex] - dFluxRho) > epsLarge) {
2708  parTestWENO = false;
2709  std::cout << "RHS check failed: rho = " << rhsBuffers[0][bufferIndex]
2710  << ", expected " << dFluxRho << std::endl;
2711  }
2712  for(int vDim = 0;vDim < numDim;vDim++){
2713  if(vDim == jumpDim){
2714  if(std::abs(rhsBuffers[vDim+1][bufferIndex] - dFluxRhoVMain) > epsLarge) {
2715  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
2716  << ", expected " << dFluxRhoVMain << std::endl;
2717  parTestWENO = false;
2718  }
2719  } else {
2720  if(std::abs(rhsBuffers[vDim+1][bufferIndex] - dFluxRhoVAlt) > epsLarge) {
2721  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
2722  << ", expected " << dFluxRhoVAlt << std::endl;
2723  parTestWENO = false;
2724  }
2725  }
2726  }
2727  if(std::abs(rhsBuffers[numDim+1][bufferIndex] - dFluxRhoE) > epsLarge) {
2728  std::cout << "RHS check failed: rhoE = "
2729  << rhsBuffers[numDim+1][bufferIndex]
2730  << ", expected " << dFluxRhoE << std::endl;
2731  parTestWENO = false;
2732  }
2733  } else if ((std::abs(cloc - checkLoc1) < stenWidth*dx+eps)
2734  || (std::abs(cloc - checkLoc2) < stenWidth*dx+eps)) {
2735  // near discontinuities WENO can have larger values of "0"
2736  // (i.e. numerical dissipation)
2737  for(int eqn = 0; eqn < numDim+2; eqn++){
2738  if(std::abs(rhsBuffers[eqn][bufferIndex]) > epsLarge)
2739  {
2740  parTestWENO = false;
2741  std::cout << "Failed RHS " << rhsBuffers[eqn][bufferIndex] << std::endl;
2742  }
2743  }
2744  } else {
2745  for(int eqn = 0; eqn < numDim+2; eqn++){
2746  if(std::abs(rhsBuffers[eqn][bufferIndex]) > eps)
2747  {
2748  parTestWENO = false;
2749  std::cout << "Failed RHS " << rhsBuffers[eqn][bufferIndex] << std::endl;
2750  }
2751  }
2752  } // checkLoc if statement
2753 
2754  if (!parTestWENO) {
2755  std::cout << "i = " << iIndex << ", j = " << jIndex << ", k = " << kIndex
2756  << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
2757  std::cout << "checkLoc1 = " << checkLoc1 << ", checkLoc2 = " << checkLoc2 << std::endl;
2758  std::cout << "rho = " << rho[bufferIndex];
2759  for (int vDim=0; vDim<numDim; vDim++) {
2760  std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
2761  }
2762  std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
2763  }
2764  }
2765  }
2766 
2767  parallelUnitResults.UpdateResult("WENO:RHS:Jump:CheckRHS",parTestWENO);
2768 
2769  if(!parTestWENO)
2770  return;
2771 
2772  } // Loop over dimensions
2773  } // jumpDim loop
2774  } // Loop over numDim = 2,3
2775 
2776  // --- TransonicJump ---
2777  // Tests that when states have initial discontinuity, RHS is
2778  // computed correctly when there is a transonic rarefaction and entropy fix
2779  // is needed for Roe scheme
2780 
2781  for(int numDim = 2;numDim < 4;numDim++){
2782 
2783  std::cout << "TestWENO:RHS:TransonicJump:Testing WENO in " << numDim << " dimensions." << std::endl;
2784 
2785  for (int jumpDim=0; jumpDim<numDim; jumpDim++) {
2786  grid_t simGrid;
2787  operator_t sbpOperator;
2788  halo_t &simHalo(simGrid.Halo());
2789  state_t simState;
2790  state_t paramState;
2791  rhs_t rhsWENO;
2792  weno_t &coeffsWENO(rhsWENO.WENOCoefficients());
2793 
2794  bool parTestWENO = true;
2795  int numPoints = 21;
2796  double length = 20.0;
2797  std::vector<size_t> gridSizes(numDim,numPoints);
2798  double dx = length/(numPoints - 1);
2799 
2800  if(testfixtures::CreateSimulationFixtures(simGrid,simHalo,sbpOperator,coeffsWENO,
2801  gridSizes,2,testComm,&messageStream)){
2802  parTestWENO = false;
2803  }
2804 
2805  if(simGrid.GenerateCoordinates(messageStream))
2806  parTestWENO = false;
2807 
2808  std::vector<double> &gridCoords(simGrid.CoordinateData());
2809 
2810  fixtures::CommunicatorType *gridCommPtr;
2811  std::vector<int> cartNeighbors;
2812  gridCommPtr = &simGrid.Communicator();
2813  gridCommPtr->CartNeighbors(cartNeighbors);
2814 
2815  // Fill halo regions with coordinates
2816  int xyzMessageId = simHalo.CreateMessageBuffers(numDim);
2817  simHalo.PackMessageBuffers(xyzMessageId,0,numDim,&gridCoords[0]);
2818 
2819  simHalo.SendMessage(0,cartNeighbors,testComm);
2820  simHalo.ReceiveMessage(0,cartNeighbors,testComm);
2821  simHalo.UnPackMessageBuffers(0,0,numDim,&gridCoords[0]);
2822 
2823  parallelUnitResults.UpdateResult("WENO:RHS:TransonicJump:InitFixtures",parTestWENO);
2824 
2825  if(!parTestWENO)
2826  return;
2827 
2828  pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());
2829  pcpp::IndexIntervalType &partitionInterval(simGrid.PartitionInterval());
2830  std::vector<size_t> &bufferSizes(simGrid.BufferSizes());
2831 
2832  if(testfixtures::euler::SetupEulerState(simGrid,simState,paramState,0)){
2833  parTestWENO = false;
2834  }
2835 
2836  parallelUnitResults.UpdateResult("WENO:RHS:TransonicJump:SetupState",parTestWENO);
2837 
2838  if(!parTestWENO)
2839  return;
2840 
2841  std::cout << messageStream.str() << std::endl;
2842  std::cout << "Testing WENO jump in " << jumpDim << " dimension." << std::endl;
2843 
2844  size_t numPointsBuffer = simGrid.BufferSize();
2845 
2846  double t = 0.0;
2847  double gamma = 1.4;
2848  double dt = .001;
2849  double cfl = 1.0;
2850  double Re = 0.0;
2851  std::vector<double> numbers(4,0.0);
2852  int flags = USEWENO;
2853 
2854  // Add some fields that are not usual
2855  paramState.Create(numPointsBuffer,0);
2856 
2857  paramState.SetFieldBuffer("gamma",&gamma);
2858  paramState.SetFieldBuffer("inputDT",&dt);
2859  paramState.SetFieldBuffer("inputCFL",&cfl);
2860  paramState.SetFieldBuffer("refRe",&Re);
2861  paramState.SetFieldBuffer("Flags",&flags);
2862  paramState.SetFieldBuffer("Numbers",&numbers[0]);
2863 
2864  std::vector<double> rho(numPointsBuffer,-1000.0);
2865  std::vector<double> rhoV(numDim*numPointsBuffer,-1000.0);
2866  std::vector<double> rhoE(numPointsBuffer,-1000.0);
2867 
2868  // "left" state (edges of domain in direction jumpDim)
2869  // corresponding primitive variables:
2870  // rho = 1, u = 1, v = w = 0, p = 1.6
2871  double rhoL = 1.0;
2872  double rhoVL = 1.0;
2873  double rhoEL = 4.5;
2874 
2875  // "right" state (center of domain in direction jumpDim)
2876  // corresponding primitive variables:
2877  // rho = 0.1, u = 0, v = w = 0, p = 0.4
2878  double rhoR = 0.1;
2879  double rhoVR = 0.0;
2880  double rhoER = 1.0;
2881 
2882  // locations of jumps, using physical coordinates
2883  double loc1 = 0.25*length;
2884  double loc2 = 0.75*length;
2885 
2886  // locations of jumps, using physical coordinates
2887  double checkLoc1 = loc1;
2888  double checkLoc2 = loc2 + dx;
2889 
2890  size_t iStart = partitionBufferInterval[0].first;
2891  size_t iEnd = partitionBufferInterval[0].second;
2892  size_t jStart = partitionBufferInterval[1].first;
2893  size_t jEnd = partitionBufferInterval[1].second;
2894 
2895  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
2896  // computed correctly
2897  size_t kStart = 0;
2898  size_t kEnd = 0;
2899  if (numDim == 3) {
2900  kStart = partitionBufferInterval[2].first;
2901  kEnd = partitionBufferInterval[2].second;
2902  }
2903 
2904  for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
2905  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
2906  for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
2907  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
2908  for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
2909  size_t bufferIndex = jBufferIndex + iIndex;
2910 
2911  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
2912 
2913  if (cloc < loc1 || cloc > loc2) {
2914  rho[bufferIndex] = rhoL;
2915  for (int vDim=0; vDim<numDim; vDim++) {
2916  if (vDim == jumpDim) rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVL;
2917  else rhoV[bufferIndex + vDim*numPointsBuffer] = 0.0;
2918  }
2919  rhoE[bufferIndex] = rhoEL;
2920  } else {
2921  rho[bufferIndex] = rhoR;
2922  for (int vDim=0; vDim<numDim; vDim++) {
2923  if (vDim == jumpDim) rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVR;
2924  else rhoV[bufferIndex + vDim*numPointsBuffer] = 0.0;
2925  }
2926  rhoE[bufferIndex] = rhoER;
2927  }
2928  }
2929  }
2930  }
2931 
2932  std::vector<double> dRho(numPointsBuffer,-1000.0);
2933  std::vector<double> dRhoV(numDim*numPointsBuffer,-1000.0);
2934  std::vector<double> dRhoE(numPointsBuffer,-1000.0);
2935 
2936  std::vector<double> T(numPointsBuffer,0.0);
2937  std::vector<double> p(numPointsBuffer,0.0);
2938  std::vector<double> V(numDim*numPointsBuffer,0.0);
2939 
2940  std::vector<double> rhom1(numPointsBuffer,0.0);
2941 
2942  simState.SetFieldBuffer("simTime",&t);
2943  simState.SetFieldBuffer("rho",&rho[0]);
2944  simState.SetFieldBuffer("rhoV",&rhoV[0]);
2945  simState.SetFieldBuffer("rhoE",&rhoE[0]);
2946 
2947  simState.SetFieldBuffer("pressure",&p[0]);
2948  simState.SetFieldBuffer("temperature",&T[0]);
2949  simState.SetFieldBuffer("velocity",&V[0]);
2950  simState.SetFieldBuffer("rhom1",&rhom1[0]);
2951 
2952  state_t rhsState(simState);
2953 
2954  rhsState.SetFieldBuffer("rho",&dRho[0]);
2955  rhsState.SetFieldBuffer("rhoV",&dRhoV[0]);
2956  rhsState.SetFieldBuffer("rhoE",&dRhoE[0]);
2957 
2958  // State/Soln initialization
2959  std::vector<double> inParams;
2960  std::vector<int> inFlags;
2961 
2962  if(rhsWENO.Initialize(simGrid,simState,paramState,sbpOperator)){
2963  parTestWENO = false;
2964  }
2965  std::vector<int> numThreadsDim;
2966  if(simGrid.SetupThreads(numThreadsDim)){
2967  std::cout << "Setting up threads for grid object failed." << std::endl;
2968  return;
2969  }
2970  rhsWENO.InitThreadIntervals();
2971 
2972  parallelUnitResults.UpdateResult("WENO:RHS:TransonicJump:InitRHS",parTestWENO);
2973 
2974  if(!parTestWENO)
2975  return;
2976 
2977  if(rhsWENO.SetupRHS(simGrid,simState,rhsState)){
2978  parTestWENO = false;
2979  }
2980 
2981  parallelUnitResults.UpdateResult("WENO:RHS:TransonicJump:SetupRHS",parTestWENO);
2982 
2983  if(!parTestWENO)
2984  return;
2985 
2986  int statusWENO = rhsWENO.RHS(0);
2987  if(statusWENO) {
2988  parTestWENO = false;
2989  std::cout << "RHS returned error code " << statusWENO << std::endl;
2990  }
2991 
2992  parallelUnitResults.UpdateResult("WENO:RHS:TransonicJump:ErrorRHS",parTestWENO);
2993 
2994  if(!parTestWENO)
2995  return;
2996 
2997  std::vector<double* > &rhsBuffers(rhsWENO.RHSBuffers());
2998 
2999  // // Print data along a pencil for debugging purposes
3000  // if (jumpDim == 0) {
3001  // size_t kBufferIndex = kStart*bufferSizes[0]*bufferSizes[1];
3002  // size_t jBufferIndex = jStart*bufferSizes[0] + kBufferIndex;
3003  // for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
3004  // size_t bufferIndex = jBufferIndex + iIndex;
3005  //
3006  // double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
3007  //
3008  // //std::cout << "i = " << iIndex << ", j = " << jStart << ", k = " << kStart
3009  // // << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
3010  // //std::cout << "checkLoc1 = " << checkLoc1 << ", checkLoc2 = " << checkLoc2 << std::endl;
3011  // //std::cout << "rho = " << rho[bufferIndex];
3012  // //for (int vDim=0; vDim<numDim; vDim++) {
3013  // // std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
3014  // //}
3015  // //std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
3016  //
3017  // std::cout << "cloc = " << cloc << ", ";
3018  //
3019  // std::cout << "rhoRHS = " << rhsBuffers[0][bufferIndex];
3020  // for (int vDim=0; vDim<numDim; vDim++) {
3021  // std::cout << ", rhoVRHS[" << vDim << "] = " << rhsBuffers[vDim+1][bufferIndex];
3022  // }
3023  // std::cout << ", rhoERHS = " << rhsBuffers[numDim+1][bufferIndex] << std::endl;
3024  // }
3025  // }
3026 
3027  // -- setup expected values
3028 
3029  // larger values
3030  double dFluxLargeRho = 1.149/dx;
3031  double dFluxLargeRhoVMain = 2.130/dx;
3032  double dFluxLargeRhoVAlt = 0.0/dx;
3033  double dFluxLargeRhoE = 6.807/dx;
3034 
3035  // smaller values
3036  double dFluxSmallRho = -0.149/dx;
3037  double dFluxSmallRhoVMain = 0.0698/dx;
3038  double dFluxSmallRhoVAlt = 0.0/dx;
3039  double dFluxSmallRhoE = -0.707/dx;
3040 
3041  // Check RHS -- should be zero everywhere except at initial discontinuity
3042  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
3043  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
3044  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
3045  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
3046  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
3047  size_t bufferIndex = jBufferIndex + iIndex;
3048 
3049  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
3050 
3051  if (std::abs(cloc - checkLoc1) < eps) {
3052  if(std::abs(rhsBuffers[0][bufferIndex] - dFluxLargeRho) > epsLarge) {
3053  parTestWENO = false;
3054  std::cout << "RHS check failed: rho = " << rhsBuffers[0][bufferIndex]
3055  << ", expected " << dFluxLargeRho << std::endl;
3056  }
3057  for(int vDim = 0;vDim < numDim;vDim++){
3058  if(vDim == jumpDim){
3059  if(std::abs(rhsBuffers[vDim+1][bufferIndex] - dFluxLargeRhoVMain) > epsLarge) {
3060  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3061  << ", expected " << dFluxLargeRhoVMain << std::endl;
3062  parTestWENO = false;
3063  }
3064  } else {
3065  if(std::abs(rhsBuffers[vDim+1][bufferIndex] - dFluxLargeRhoVAlt) > epsLarge) {
3066  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3067  << ", expected " << dFluxLargeRhoVAlt << std::endl;
3068  parTestWENO = false;
3069  }
3070  }
3071  }
3072  if(std::abs(rhsBuffers[numDim+1][bufferIndex] -
3073  dFluxLargeRhoE) > epsLarge) {
3074  std::cout << "RHS check failed: rhoE = "
3075  << rhsBuffers[numDim+1][bufferIndex]
3076  << ", expected " << dFluxLargeRhoE << std::endl;
3077  parTestWENO = false;
3078  }
3079  } else if (std::abs(cloc - (checkLoc1-dx)) < eps) {
3080  if(std::abs(rhsBuffers[0][bufferIndex] - dFluxSmallRho) > epsLarge) {
3081  parTestWENO = false;
3082  std::cout << "RHS check failed: rho = " << rhsBuffers[0][bufferIndex]
3083  << ", expected " << dFluxSmallRho << std::endl;
3084  }
3085  for(int vDim = 0;vDim < numDim;vDim++){
3086  if(vDim == jumpDim){
3087  if(std::abs(rhsBuffers[vDim+1][bufferIndex] - dFluxSmallRhoVMain) > epsLarge) {
3088  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3089  << ", expected " << dFluxSmallRhoVMain << std::endl;
3090  parTestWENO = false;
3091  }
3092  } else {
3093  if(std::abs(rhsBuffers[vDim+1][bufferIndex] - dFluxSmallRhoVAlt) > epsLarge) {
3094  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3095  << ", expected " << dFluxSmallRhoVAlt << std::endl;
3096  parTestWENO = false;
3097  }
3098  }
3099  }
3100  if(std::abs(rhsBuffers[numDim+1][bufferIndex] -
3101  dFluxSmallRhoE) > epsLarge) {
3102  std::cout << "RHS check failed: rhoE = "
3103  << rhsBuffers[numDim+1][bufferIndex]
3104  << ", expected " << dFluxSmallRhoE << std::endl;
3105  parTestWENO = false;
3106  }
3107  } else if (std::abs(cloc - checkLoc2) < eps) {
3108  if(std::abs(rhsBuffers[0][bufferIndex] + dFluxLargeRho) > epsLarge) {
3109  parTestWENO = false;
3110  std::cout << "RHS check failed: rho = " << rhsBuffers[0][bufferIndex]
3111  << ", expected " << -dFluxLargeRho << std::endl;
3112  }
3113  for(int vDim = 0;vDim < numDim;vDim++){
3114  if(vDim == jumpDim){
3115  if(std::abs(rhsBuffers[vDim+1][bufferIndex] + dFluxLargeRhoVMain) > epsLarge) {
3116  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3117  << ", expected " << -dFluxLargeRhoVMain << std::endl;
3118  parTestWENO = false;
3119  }
3120  } else {
3121  if(std::abs(rhsBuffers[vDim+1][bufferIndex] + dFluxLargeRhoVAlt) > epsLarge) {
3122  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3123  << ", expected " << -dFluxLargeRhoVAlt << std::endl;
3124  parTestWENO = false;
3125  }
3126  }
3127  }
3128  if(std::abs(rhsBuffers[numDim+1][bufferIndex] +
3129  dFluxLargeRhoE) > epsLarge) {
3130  std::cout << "RHS check failed: rhoE = "
3131  << rhsBuffers[numDim+1][bufferIndex]
3132  << ", expected " << -dFluxLargeRhoE << std::endl;
3133  parTestWENO = false;
3134  }
3135  } else if (std::abs(cloc - (checkLoc2-dx)) < eps) {
3136  if(std::abs(rhsBuffers[0][bufferIndex] + dFluxSmallRho) > epsLarge) {
3137  parTestWENO = false;
3138  std::cout << "RHS check failed: rho = " << rhsBuffers[0][bufferIndex]
3139  << ", expected " << -dFluxSmallRho << std::endl;
3140  }
3141  for(int vDim = 0;vDim < numDim;vDim++){
3142  if(vDim == jumpDim){
3143  if(std::abs(rhsBuffers[vDim+1][bufferIndex] + dFluxSmallRhoVMain) > epsLarge) {
3144  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3145  << ", expected " << -dFluxSmallRhoVMain << std::endl;
3146  parTestWENO = false;
3147  }
3148  } else {
3149  if(std::abs(rhsBuffers[vDim+1][bufferIndex] + dFluxSmallRhoVAlt) > epsLarge) {
3150  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3151  << ", expected " << -dFluxSmallRhoVAlt << std::endl;
3152  parTestWENO = false;
3153  }
3154  }
3155  }
3156  if(std::abs(rhsBuffers[numDim+1][bufferIndex] +
3157  dFluxSmallRhoE) > epsLarge) {
3158  std::cout << "RHS check failed: rhoE = "
3159  << rhsBuffers[numDim+1][bufferIndex]
3160  << ", expected " << -dFluxSmallRhoE << std::endl;
3161  parTestWENO = false;
3162  }
3163  } else if ((std::abs(cloc - checkLoc1) < stenWidth*dx+eps)
3164  || (std::abs(cloc - checkLoc2) < stenWidth*dx+eps)) {
3165  // near discontinuities WENO can have larger values of "0"
3166  // (i.e. numerical dissipation)
3167  for(int eqn = 0; eqn < numDim+2; eqn++){
3168  if(std::abs(rhsBuffers[eqn][bufferIndex]) > epsLarge)
3169  {
3170  parTestWENO = false;
3171  std::cout << "RHS check failed " << rhsBuffers[eqn][bufferIndex] << std::endl;
3172  }
3173  }
3174  } else {
3175  for(int eqn = 0; eqn < numDim+2; eqn++){
3176  if(std::abs(rhsBuffers[eqn][bufferIndex]) > eps)
3177  {
3178  parTestWENO = false;
3179  std::cout << "RHS check failed " << rhsBuffers[eqn][bufferIndex] << std::endl;
3180  }
3181  }
3182  } // checkLoc if statement
3183 
3184  if (!parTestWENO) {
3185  std::cout << "i = " << iIndex << ", j = " << jIndex << ", k = " << kIndex
3186  << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
3187  std::cout << "checkLoc1 = " << checkLoc1 << ", checkLoc2 = " << checkLoc2 << std::endl;
3188  std::cout << "rho = " << rho[bufferIndex];
3189  for (int vDim=0; vDim<numDim; vDim++) {
3190  std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
3191  }
3192  std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
3193  }
3194  }
3195  }
3196 
3197  parallelUnitResults.UpdateResult("WENO:RHS:TransonicJump:CheckRHS",parTestWENO);
3198 
3199  if(!parTestWENO)
3200  return;
3201 
3202  } // Loop over dimensions
3203  } // jumpDim loop
3204  } // Loop over numDim = 2,3
3205 
3206  // --- JumpMetric ---
3207  // Tests that when states have initial discontinuity, RHS is
3208  // computed correctly
3209  // Grid setup so that metric terms are involved
3210 
3211  for(int numDim = 2;numDim < 4;numDim++){
3212 
3213  std::cout << "TestWENO:RHS:JumpMetric:Testing WENO in " << numDim << " dimensions." << std::endl;
3214 
3215  for (int jumpDim=0; jumpDim<numDim; jumpDim++) {
3216  grid_t simGrid;
3217  operator_t sbpOperator;
3218  halo_t &simHalo(simGrid.Halo());
3219  state_t simState;
3220  state_t paramState;
3221  rhs_t rhsWENO;
3222  weno_t &coeffsWENO(rhsWENO.WENOCoefficients());
3223 
3224  bool parTestWENO = true;
3225  int numPoints = 41;
3226  double length = 20.0;
3227  std::vector<size_t> gridSizes(numDim,numPoints);
3228  double dx = length/(numPoints - 1);
3229 
3230  if(testfixtures::CreateSimulationFixtures(simGrid,simHalo,sbpOperator,coeffsWENO,
3231  gridSizes,2,testComm,&messageStream)){
3232  parTestWENO = false;
3233  }
3234 
3235  if(simGrid.GenerateCoordinates(messageStream))
3236  parTestWENO = false;
3237 
3238  std::vector<double> &gridCoords(simGrid.CoordinateData());
3239 
3240  fixtures::CommunicatorType *gridCommPtr;
3241  std::vector<int> cartNeighbors;
3242  gridCommPtr = &simGrid.Communicator();
3243  gridCommPtr->CartNeighbors(cartNeighbors);
3244 
3245  // Fill halo regions with coordinates
3246  int xyzMessageId = simHalo.CreateMessageBuffers(numDim);
3247  simHalo.PackMessageBuffers(xyzMessageId,0,numDim,&gridCoords[0]);
3248 
3249  simHalo.SendMessage(0,cartNeighbors,testComm);
3250  simHalo.ReceiveMessage(0,cartNeighbors,testComm);
3251  simHalo.UnPackMessageBuffers(0,0,numDim,&gridCoords[0]);
3252 
3253  parallelUnitResults.UpdateResult("WENO:RHS:JumpMetric:InitFixtures",parTestWENO);
3254 
3255  if(!parTestWENO)
3256  return;
3257 
3258  pcpp::IndexIntervalType &partitionBufferInterval(simGrid.PartitionBufferInterval());
3259  pcpp::IndexIntervalType &partitionInterval(simGrid.PartitionInterval());
3260  std::vector<size_t> &bufferSizes(simGrid.BufferSizes());
3261 
3262  if(testfixtures::euler::SetupEulerState(simGrid,simState,paramState,0)){
3263  parTestWENO = false;
3264  }
3265 
3266  parallelUnitResults.UpdateResult("WENO:RHS:JumpMetric:SetupState",parTestWENO);
3267 
3268  if(!parTestWENO)
3269  return;
3270 
3271  std::cout << messageStream.str() << std::endl;
3272  std::cout << "Testing WENO jump in " << jumpDim << " dimension." << std::endl;
3273 
3274  size_t numPointsBuffer = simGrid.BufferSize();
3275 
3276  double t = 0.0;
3277  double gamma = 1.4;
3278  double dt = .001;
3279  double cfl = 1.0;
3280  double Re = 0.0;
3281  std::vector<double> numbers(4,0.0);
3282  int flags = USEWENO;
3283 
3284  // Add some fields that are not usual
3285  paramState.Create(numPointsBuffer,0);
3286 
3287  paramState.SetFieldBuffer("gamma",&gamma);
3288  paramState.SetFieldBuffer("inputDT",&dt);
3289  paramState.SetFieldBuffer("inputCFL",&cfl);
3290  paramState.SetFieldBuffer("refRe",&Re);
3291  paramState.SetFieldBuffer("Flags",&flags);
3292  paramState.SetFieldBuffer("Numbers",&numbers[0]);
3293 
3294  std::vector<double> rho(numPointsBuffer,-1000.0);
3295  std::vector<double> rhoV(numDim*numPointsBuffer,-1000.0);
3296  std::vector<double> rhoE(numPointsBuffer,-1000.0);
3297 
3298  // "left" state (edges of domain in direction jumpDim)
3299  double rhoL = 1.0;
3300  double rhoVL = 2.0;
3301  double rhoEL = 4.0 + 2.0*numDim;
3302 
3303  // "right" state (center of domain in direction jumpDim)
3304  double rhoR = 2.0;
3305  double rhoVR = 6.0;
3306  double rhoER = 8.0 + 9.0*numDim;
3307 
3308  // locations of jumps, using physical coordinates
3309  double loc1 = 0.25*length;
3310  double loc2 = 0.75*length;
3311 
3312  // locations of jumps, using physical coordinates
3313  double checkLoc1 = loc1;
3314  double checkLoc2 = loc2 + dx;
3315 
3316  size_t iStart = partitionBufferInterval[0].first;
3317  size_t iEnd = partitionBufferInterval[0].second;
3318  size_t jStart = partitionBufferInterval[1].first;
3319  size_t jEnd = partitionBufferInterval[1].second;
3320 
3321  // if numDim == 2 this will set kIndex = 0, so bufferIndex will still be
3322  // computed correctly
3323  size_t kStart = 0;
3324  size_t kEnd = 0;
3325  if (numDim == 3) {
3326  kStart = partitionBufferInterval[2].first;
3327  kEnd = partitionBufferInterval[2].second;
3328  }
3329 
3330  for (size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
3331  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
3332  for (size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
3333  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
3334  for (size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
3335  size_t bufferIndex = jBufferIndex + iIndex;
3336 
3337  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
3338 
3339  if (cloc < loc1 || cloc > loc2) {
3340  rho[bufferIndex] = rhoL;
3341  for (int vDim=0; vDim<numDim; vDim++) {
3342  rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVL;
3343  }
3344  rhoE[bufferIndex] = rhoEL;
3345  } else {
3346  rho[bufferIndex] = rhoR;
3347  for (int vDim=0; vDim<numDim; vDim++) {
3348  rhoV[bufferIndex + vDim*numPointsBuffer] = rhoVR;
3349  }
3350  rhoE[bufferIndex] = rhoER;
3351  }
3352  }
3353  }
3354  }
3355 
3356  std::vector<double> dRho(numPointsBuffer,-1000.0);
3357  std::vector<double> dRhoV(numDim*numPointsBuffer,-1000.0);
3358  std::vector<double> dRhoE(numPointsBuffer,-1000.0);
3359 
3360  std::vector<double> T(numPointsBuffer,0.0);
3361  std::vector<double> p(numPointsBuffer,0.0);
3362  std::vector<double> V(numDim*numPointsBuffer,0.0);
3363 
3364  std::vector<double> rhom1(numPointsBuffer,0.0);
3365 
3366  simState.SetFieldBuffer("simTime",&t);
3367  simState.SetFieldBuffer("rho",&rho[0]);
3368  simState.SetFieldBuffer("rhoV",&rhoV[0]);
3369  simState.SetFieldBuffer("rhoE",&rhoE[0]);
3370 
3371  simState.SetFieldBuffer("pressure",&p[0]);
3372  simState.SetFieldBuffer("temperature",&T[0]);
3373  simState.SetFieldBuffer("velocity",&V[0]);
3374  simState.SetFieldBuffer("rhom1",&rhom1[0]);
3375 
3376  state_t rhsState(simState);
3377 
3378  rhsState.SetFieldBuffer("rho",&dRho[0]);
3379  rhsState.SetFieldBuffer("rhoV",&dRhoV[0]);
3380  rhsState.SetFieldBuffer("rhoE",&dRhoE[0]);
3381 
3382  // State/Soln initialization
3383  std::vector<double> inParams;
3384  std::vector<int> inFlags;
3385 
3386  if(rhsWENO.Initialize(simGrid,simState,paramState,sbpOperator)){
3387  parTestWENO = false;
3388  }
3389  std::vector<int> numThreadsDim;
3390  if(simGrid.SetupThreads(numThreadsDim)){
3391  std::cout << "Setting up threads for grid object failed." << std::endl;
3392  return;
3393  }
3394  rhsWENO.InitThreadIntervals();
3395 
3396  parallelUnitResults.UpdateResult("WENO:RHS:JumpMetric:InitRHS",parTestWENO);
3397 
3398  if(!parTestWENO)
3399  return;
3400 
3401  if(rhsWENO.SetupRHS(simGrid,simState,rhsState)){
3402  parTestWENO = false;
3403  }
3404 
3405  parallelUnitResults.UpdateResult("WENO:RHS:JumpMetric:SetupRHS",parTestWENO);
3406 
3407  if(!parTestWENO)
3408  return;
3409 
3410  int statusWENO = rhsWENO.RHS(0);
3411  if(statusWENO) {
3412  parTestWENO = false;
3413  std::cout << "RHS returned error code " << statusWENO << std::endl;
3414  }
3415 
3416  parallelUnitResults.UpdateResult("WENO:RHS:JumpMetric:ErrorRHS",parTestWENO);
3417 
3418  if(!parTestWENO)
3419  return;
3420 
3421  std::vector<double* > &rhsBuffers(rhsWENO.RHSBuffers());
3422 
3423  // -- setup expected values
3424 
3425  // "left" state (edges of domain in direction jumpDim)
3426  double fluxRhoL = 2.0;
3427  double fluxRhoVMainL = 5.6;
3428  double fluxRhoVAltL = 4.0;
3429  double fluxRhoEL = 11.2 + 4.0*numDim;
3430 
3431  // "right" state (center of domain in direction jumpDim)
3432  double fluxRhoR = 6.0;
3433  double fluxRhoVMainR = 21.2;
3434  double fluxRhoVAltR = 18.0;
3435  double fluxRhoER = 33.6 + 27.0*numDim;
3436 
3437  double dFluxRho = (fluxRhoR - fluxRhoL)/dx;
3438  double dFluxRhoVMain = (fluxRhoVMainR - fluxRhoVMainL)/dx;
3439  double dFluxRhoVAlt = (fluxRhoVAltR - fluxRhoVAltL)/dx;
3440  double dFluxRhoE = (fluxRhoER - fluxRhoEL)/dx;
3441 
3442 // // Print data along a pencil for debugging purposes
3443 // if (jumpDim == 0) {
3444 // size_t kBufferIndex = kStart*bufferSizes[0]*bufferSizes[1];
3445 // size_t jBufferIndex = jStart*bufferSizes[0] + kBufferIndex;
3446 // for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
3447 // size_t bufferIndex = jBufferIndex + iIndex;
3448 //
3449 // double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
3450 //
3451 // //std::cout << "i = " << iIndex << ", j = " << jStart << ", k = " << kStart
3452 // // << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
3453 // //std::cout << "checkLoc1 = " << checkLoc1 << ", checkLoc2 = " << checkLoc2 << std::endl;
3454 // //std::cout << "rho = " << rho[bufferIndex];
3455 // //for (int vDim=0; vDim<numDim; vDim++) {
3456 // // std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
3457 // //}
3458 // //std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
3459 //
3460 // std::cout << "cloc = " << cloc << ", ";
3461 //
3462 // std::cout << "rhoRHS = " << rhsBuffers[0][bufferIndex];
3463 // for (int vDim=0; vDim<numDim; vDim++) {
3464 // std::cout << ", rhoVRHS[" << vDim << "] = " << rhsBuffers[vDim+1][bufferIndex];
3465 // }
3466 // std::cout << ", rhoERHS = " << rhsBuffers[numDim+1][bufferIndex] << std::endl;
3467 // }
3468 // }
3469 
3470  // Check RHS -- should be zero everywhere except at initial discontinuity
3471  for(size_t kIndex = kStart;((kIndex <= kEnd)&&parTestWENO);kIndex++){
3472  size_t kBufferIndex = kIndex*bufferSizes[0]*bufferSizes[1];
3473  for(size_t jIndex = jStart;((jIndex <= jEnd)&&parTestWENO);jIndex++){
3474  size_t jBufferIndex = jIndex*bufferSizes[0] + kBufferIndex;
3475  for(size_t iIndex = iStart;((iIndex <= iEnd) && parTestWENO);iIndex++){
3476  size_t bufferIndex = jBufferIndex + iIndex;
3477 
3478  double cloc = gridCoords[bufferIndex + jumpDim*numPointsBuffer];
3479 
3480  if (std::abs(cloc - checkLoc1) < eps) {
3481  if(std::abs(rhsBuffers[0][bufferIndex] + dFluxRho) > epsLarge) {
3482  parTestWENO = false;
3483  std::cout << "RHS check failed: rho = " << rhsBuffers[0][bufferIndex]
3484  << ", expected " << -dFluxRho << std::endl;
3485  }
3486  for(int vDim = 0;vDim < numDim;vDim++){
3487  if(vDim == jumpDim){
3488  if(std::abs(rhsBuffers[vDim+1][bufferIndex] + dFluxRhoVMain) > epsLarge) {
3489  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3490  << ", expected " << -dFluxRhoVMain << std::endl;
3491  parTestWENO = false;
3492  }
3493  } else {
3494  if(std::abs(rhsBuffers[vDim+1][bufferIndex] + dFluxRhoVAlt) > epsLarge) {
3495  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3496  << ", expected " << -dFluxRhoVAlt << std::endl;
3497  parTestWENO = false;
3498  }
3499  }
3500  }
3501  if(std::abs(rhsBuffers[numDim+1][bufferIndex] + dFluxRhoE) > epsLarge) {
3502  std::cout << "RHS check failed: rhoE = "
3503  << rhsBuffers[numDim+1][bufferIndex]
3504  << ", expected " << -dFluxRhoE << std::endl;
3505  parTestWENO = false;
3506  }
3507  } else if (std::abs(cloc - checkLoc2) < eps) {
3508  if(std::abs(rhsBuffers[0][bufferIndex] - dFluxRho) > epsLarge) {
3509  parTestWENO = false;
3510  std::cout << "RHS check failed: rho = " << rhsBuffers[0][bufferIndex]
3511  << ", expected " << dFluxRho << std::endl;
3512  }
3513  for(int vDim = 0;vDim < numDim;vDim++){
3514  if(vDim == jumpDim){
3515  if(std::abs(rhsBuffers[vDim+1][bufferIndex] - dFluxRhoVMain) > epsLarge) {
3516  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3517  << ", expected " << dFluxRhoVMain << std::endl;
3518  parTestWENO = false;
3519  }
3520  } else {
3521  if(std::abs(rhsBuffers[vDim+1][bufferIndex] - dFluxRhoVAlt) > epsLarge) {
3522  std::cout << "RHS check failed: rhoV = " << rhsBuffers[vDim+1][bufferIndex]
3523  << ", expected " << dFluxRhoVAlt << std::endl;
3524  parTestWENO = false;
3525  }
3526  }
3527  }
3528  if(std::abs(rhsBuffers[numDim+1][bufferIndex] - dFluxRhoE) > epsLarge) {
3529  std::cout << "RHS check failed: rhoE = "
3530  << rhsBuffers[numDim+1][bufferIndex]
3531  << ", expected " << dFluxRhoE << std::endl;
3532  parTestWENO = false;
3533  }
3534  } else if ((std::abs(cloc - checkLoc1) < stenWidth*dx+eps)
3535  || (std::abs(cloc - checkLoc2) < stenWidth*dx+eps)) {
3536  // near discontinuities WENO can have larger values of "0"
3537  // (i.e. numerical dissipation)
3538  for(int eqn = 0; eqn < numDim+2; eqn++){
3539  if(std::abs(rhsBuffers[eqn][bufferIndex]) > epsLarge)
3540  {
3541  parTestWENO = false;
3542  std::cout << "Failed RHS " << rhsBuffers[eqn][bufferIndex] << std::endl;
3543  }
3544  }
3545  } else {
3546  for(int eqn = 0; eqn < numDim+2; eqn++){
3547  if(std::abs(rhsBuffers[eqn][bufferIndex]) > eps)
3548  {
3549  parTestWENO = false;
3550  std::cout << "Failed RHS " << rhsBuffers[eqn][bufferIndex] << std::endl;
3551  }
3552  }
3553  } // checkLoc if statement
3554 
3555  if (!parTestWENO) {
3556  std::cout << "i = " << iIndex << ", j = " << jIndex << ", k = " << kIndex
3557  << ", jumpDim = " << jumpDim << ", cloc = " << cloc << std::endl;
3558  std::cout << "checkLoc1 = " << checkLoc1 << ", checkLoc2 = " << checkLoc2 << std::endl;
3559  std::cout << "rho = " << rho[bufferIndex];
3560  for (int vDim=0; vDim<numDim; vDim++) {
3561  std::cout << ", rhoV[" << vDim << "] = " << rhoV[bufferIndex + vDim*numPointsBuffer];
3562  }
3563  std::cout << ", rhoE = " << rhoE[bufferIndex] << std::endl;
3564  }
3565  }
3566  }
3567 
3568  parallelUnitResults.UpdateResult("WENO:RHS:JumpMetric:CheckRHS",parTestWENO);
3569 
3570  if(!parTestWENO)
3571  return;
3572 
3573  } // Loop over dimensions
3574  } // jumpDim loop
3575  } // Loop over numDim = 2,3
3576 };
plascom2::state_t
simulation::state::base state_t
Definition: PlasCom2.H:17
numPoints
void const size_t * numPoints
Definition: EulerKernels.H:10
operator_t
plascom2::operators::sbp::base operator_t
Definition: TestMaxwellSolver.C:1568
simulation::grid::parallel_blockstructured::PartitionInterval
pcpp::IndexIntervalType & PartitionInterval()
Definition: Grid.H:500
simulation::grid::parallel_blockstructured::GenerateCoordinates
int GenerateCoordinates(std::ostream &)
Definition: Grid.C:1628
gridSizes
void const size_t const size_t * gridSizes
Definition: EulerKernels.H:10
comm_t
pcpp::CommunicatorType comm_t
Definition: TestRK4Advancer.C:17
pcpp::field::dataset::Create
virtual size_t Create(size_t number_of_nodes=0, size_t number_of_cells=0)
Definition: PCPPFieldData.H:578
interval_t
pcpp::IndexIntervalType interval_t
Definition: TestMaxwellSolver.C:1567
testfixtures::euler::SetupEulerState
int SetupEulerState(const GridType &inGrid, StateType &inState, StateType &inParams, int numScalars, bool withFluxes=false)
Definition: EulerTestFixtures.H:26
global_t
pcpp::ParallelGlobalType global_t
Definition: TestHDF5.C:16
simulation::grid::parallel_blockstructured::BufferSizes
const std::vector< size_t > & BufferSizes() const
Definition: Grid.H:459
TestWENO_RHS
void TestWENO_RHS(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
Definition: TestWENOParallel.C:2189
USEWENO
static const int USEWENO
Definition: EulerRHS.H:26
ix::test::results
Encapsulating class for collections of test results.
Definition: Testing.H:18
plascom2::domain_t
simulation::domain::base< grid_t, state_t > domain_t
Definition: PlasCom2.H:19
simulation::grid::parallel_blockstructured::PartitionBufferInterval
pcpp::IndexIntervalType & PartitionBufferInterval()
Definition: Grid.H:519
ix::comm::CommunicatorObject
Main encapsulation of MPI.
Definition: COMM.H:62
TestWENO_ApplyWENO
void TestWENO_ApplyWENO(ix::test::results &parallelUnitResults, pcpp::CommunicatorType &testComm)
Definition: TestWENOParallel.C:17
Testing.H
Testing constructs for unit testing.
simulation::grid::parallel_blockstructured::SetupThreads
int SetupThreads(std::vector< int > &numThreadPartitions)
Definition: Grid.C:2174
plascom2::operators::stencilset
Encapsulation for a collection of operator stencils.
Definition: Stencil.H:26
testfixtures::CreateSimulationFixtures
int CreateSimulationFixtures(GridType &inGrid, HaloType &inHalo, plascom2::operators::sbp::base &inOperator, const std::vector< size_t > &gridSizes, int interiorOrder, pcpp::CommunicatorType &inComm, std::ostream *messageStreamPtr=NULL)
Definition: TestFixtures.H:31
simulation::grid::parallel_blockstructured::Halo
HaloType & Halo()
Definition: Grid.H:544
bufferSizes
void const size_t const size_t * bufferSizes
Definition: MetricKernels.H:19
V
void const size_t const size_t const size_t const int const double * V
Definition: OperatorKernels.H:76
ix::test::results::UpdateResult
void UpdateResult(const std::string &name, const ValueType &result)
Updates an existing test result.
Definition: Testing.H:55
simulation::grid::parallel_blockstructured::CoordinateData
std::vector< double > & CoordinateData()
Definition: Grid.H:503
plascom2::halo_t
simulation::grid::halo halo_t
Definition: PlasCom2.H:18
ix::util::sizeextent
Simple Block Structured Mesh object.
Definition: IndexUtil.H:21
simulation::grid::parallel_base::Communicator
fixtures::CommunicatorType & Communicator() const
Definition: Grid.H:350
plascom2::grid_t
simulation::grid::parallel_blockstructured grid_t
Definition: PlasCom2.H:16
pcpp::field::dataset::SetFieldBuffer
void SetFieldBuffer(const std::string &name, void *buf)
Definition: PCPPFieldData.H:716
simulation::grid::parallel_blockstructured::BufferSize
size_t BufferSize() const
Definition: Grid.H:430
numPointsBuffer
void const size_t * numPointsBuffer
Definition: MetricKernels.H:19
simulation::grid::parallel_blockstructured
Definition: Grid.H:387