PlasCom2  1.0
XPACC Multi-physics simluation application
TestRK4Advancer.C
Go to the documentation of this file.
1 #include "Testing.H"
2 #include "Simulation.H"
3 #include "PCPPCommUtil.H"
4 #include "PCPPReport.H"
5 #include "PCPPIntervalUtils.H"
6 #include "Stencil.H"
7 #include "RKTestFixtures.H"
8 #include "TestFixtures.H"
9 #include "RK4Advancer.H"
10 #include "EulerRHS.H"
11 
12 typedef simulation::grid::parallel_blockstructured grid_t;
13 typedef simulation::state::base state_t;
14 typedef plascom2::operators::sbp::base operator_t;
15 //typedef simulation::grid::halo halo_t;
16 typedef pcpp::ParallelGlobalType global_t;
17 typedef pcpp::CommunicatorType comm_t;
18 typedef pcpp::IndexIntervalType interval_t;
19 
20 void TestRK4Advancer(ix::test::results &serialUnitResults){
21 
22  std::cout << "Entering TestRK4Advancer" << std::endl;
23  // basic types for the advection problem
24  typedef simulation::domain::base<grid_t,state_t> domain_t;
25  typedef testrhs<domain_t> rhs_t;
26 
27  // -- Set up the RHS --
28  rhs_t testRHS;
29 
30  // -- Set up the grid --
31  grid_t testGrid;
32  // operator_t sbpOperator;
33  // halo_t simHalo;
34  state_t testState;
35 
36  std::ostringstream messageStream;
37  std::vector<size_t> gridSizes(3,1);
38  std::vector<size_t> numNodes(gridSizes);
39 
40  if(testfixtures::CreateSimpleGrid(testGrid,gridSizes)){
41  std::cerr << "ERROR! CreateSimpleGrid failed." << std::endl;
42  return;
43  }
44 
45  pcpp::IndexIntervalType globalExtent;
46  globalExtent.InitSimple(numNodes);
47 
48  size_t numNodesDomain = numNodes[0]*numNodes[1]*numNodes[2];
49  size_t numCellsDomain = (numNodes[0]-1)*(numNodes[1]-1)*(numNodes[2]-1);
50 
51  testState.AddField("A",'s',1,8,"");
52  testState.AddField("simTime",'s',1,8,"s");
53  testState.AddField("simDT",'s',1,8,"s");
54  testState.AddField("F0",'n',1,8,"");
55  testState.AddField("F1",'n',1,8,"");
56  testState.AddField("F2",'n',1,8,"");
57  testState.AddField("F3",'n',1,8,"");
58  testState.AddField("F4",'n',1,8,"");
59  testState.AddField("F5",'n',1,8,"");
60  testState.AddField("F6",'n',1,8,"");
61  testState.Create(numNodesDomain,numCellsDomain);
62  testState.InitializeFieldHandles();
63  testState.SetStateFields("F0 F1 F2 F3 F4 F5 F6");
64 
65  double a = 1.0;
66  double t = 1.0;
67  double f0 = a*testRHS.TestF(t,0);
68  double f1 = a*testRHS.TestF(t,1);
69  double f2 = a*testRHS.TestF(t,2);
70  double f3 = a*testRHS.TestF(t,3);
71  double f4 = a*testRHS.TestF(t,4);
72  double f5 = a*testRHS.TestF(t,5);
73  double f6 = a*testRHS.TestF(t,6);
74  double dt = 1.0;
75 
76  // set up the domain
77  domain_t rkDomain;
78  rkDomain.SetNumGrids(1);
79  // rkDomain.RHSs().resize(1,NULL);
80  rkDomain.GridNames().resize(1,"testgrid");
81  rkDomain.SetGrid(0,testGrid);
82  rkDomain.SetGridState(0,testState);
83 // rkDomain.Grids().resize(1,&testGrid);
84 // rkDomain.States().resize(1,&testState);
85  rkDomain.PartitionDomain(0);
86  state_t &domainState(rkDomain.State(0));
87 
88  std::cout << "Initing" << std::endl;
89  rkDomain.SetRHS(testRHS);
90  testRHS.SetTimestep(dt);
91 
92  domainState.SetFieldBuffer("A",&a);
93  domainState.SetFieldBuffer("simTime",&t);
94  domainState.SetFieldBuffer("simDT",&dt);
95  domainState.SetFieldBuffer("F0",&f0);
96  domainState.SetFieldBuffer("F1",&f1);
97  domainState.SetFieldBuffer("F2",&f2);
98  domainState.SetFieldBuffer("F3",&f3);
99  domainState.SetFieldBuffer("F4",&f4);
100  domainState.SetFieldBuffer("F5",&f5);
101  domainState.SetFieldBuffer("F6",&f6);
102 
103  // set up the advancer
104  rk4advancer<domain_t> testAdvancer;
105 
106 // testAdvancer.SetMessageStream(messageStream);
107 // testAdvancer.SetRHS(testRHS);
108  std::cout << "InitializeAdvancer" << std::endl;
109  global_t myGlobal;
110  testAdvancer.InitializeAdvancer(rkDomain,myGlobal,messageStream);
111 
112 
113 // messageStream << "Domain State report: " << std::endl;
114 // domainState.Report(messageStream);
115 
116  std::cout << messageStream.str() << std::endl;
117  pcpp::io::RenewStream(messageStream);
118 
119  int nTrial = 10;
120  bool isAccurate = true;
121  bool isOrder4 = true;
122  double errorTolerance = 1e-15;
123 
124  std::vector<std::vector<double> > errorValues(nTrial);
125  std::vector<double> timeSteps;
126 
127  std::cout << "Beginning trials" << std::endl;
128  for(int iTrial = 0;iTrial < 10;iTrial++){
129 
130  // Reset all the input data for the current trial
131  t = 0.5;
132  double tbefore = t;
133  dt = 1.0/(std::pow(2.0,static_cast<double>(iTrial)));
134  testRHS.SetTimestep(dt);
135 
136  f0 = a*testRHS.TestF(t,0);
137  f1 = a*testRHS.TestF(t,1);
138  f2 = a*testRHS.TestF(t,2);
139  f3 = a*testRHS.TestF(t,3);
140  f4 = a*testRHS.TestF(t,4);
141  f5 = a*testRHS.TestF(t,5);
142  f6 = a*testRHS.TestF(t,6);
143 
144  // Report to stdout some status/diagnostic
145  std::cout << "Driver:Time(before): " << t << std::endl;
146  std::cout << "Driver:TimeStep: " << dt << std::endl;
147 
148  // record the time step values
149  timeSteps.push_back(dt);
150  // domainState.Report(std::cout);
151  std::cout << "Driver:F(before): " << std::endl
152  << "Driver:f0 = " << f0 << std::endl
153  << "Driver:f1 = " << f1 << std::endl
154  << "Driver:f2 = " << f2 << std::endl
155  << "Driver:f3 = " << f3 << std::endl
156  << "Driver:f4 = " << f4 << std::endl
157  << "Driver:f5 = " << f5 << std::endl
158  << "Driver:f6 = " << f6 << std::endl;
159 
160  // Advance one timestep
161  if(testAdvancer.AdvanceDomain()){
162  break;
163  }
164 
165  if(std::abs(t - tbefore - dt) >= 1e-24){
166  std::cout << "Time did not advance, aborting test." << std::endl;
167  return;
168  // Fail test and abort
169  }
170  // sanity check that time was updated
171  // domainState.Report(std::cout);
172  std::cout << "Time(after): " << t << std::endl;
173  std::cout << "Driver:F(after): " << std::endl
174  << "Driver:f0 = " << f0 << std::endl
175  << "Driver:f1 = " << f1 << std::endl
176  << "Driver:f2 = " << f2 << std::endl
177  << "Driver:f3 = " << f3 << std::endl
178  << "Driver:f4 = " << f4 << std::endl
179  << "Driver:f5 = " << f5 << std::endl
180  << "Driver:f6 = " << f6 << std::endl;
181 
182  // create storage for the error values and record them
183  errorValues[iTrial].resize(7,0.0);
184  errorValues[iTrial][0] = std::abs(f0 - a*testRHS.TestF(t,0));
185  errorValues[iTrial][1] = std::abs(f1 - a*testRHS.TestF(t,1));
186  errorValues[iTrial][2] = std::abs(f2 - a*testRHS.TestF(t,2));
187  errorValues[iTrial][3] = std::abs(f3 - a*testRHS.TestF(t,3));
188  errorValues[iTrial][4] = std::abs(f4 - a*testRHS.TestF(t,4));
189  errorValues[iTrial][5] = std::abs(f5 - a*testRHS.TestF(t,5));
190  errorValues[iTrial][6] = std::abs(f6 - a*testRHS.TestF(t,6));
191 
192  // Make sure the fields that should be integrated exactly
193  // have essentially zero error
194  for(int iVal = 0;iVal < 5;iVal++)
195  if(errorValues[iTrial][iVal] > errorTolerance)
196  isAccurate = false;
197 
198  // sanity check to stdout
199  for(int iVal = 0;iVal < 7;iVal++){
200  std::cout << "Driver:Error(" << iVal << ") = " << errorValues[iTrial][iVal] << std::endl;
201  }
202 
203  }
204 
205  // Check the truncation error is correctly O(dt^5)
206  for(int iTrial = 0;iTrial < 8;iTrial++){
207  double logErr = std::log(errorValues[iTrial+1][5]/errorValues[iTrial][5]);
208  double logDT = std::log(timeSteps[iTrial+1]/timeSteps[iTrial]);
209  double methodOrder = logErr/logDT;
210  std::cout << "iTrial(" << iTrial << ") Order: " << methodOrder << std::endl;
211  if(std::abs(methodOrder - 5) > 1e-3)
212  isOrder4 = false;
213  }
214 
215  serialUnitResults.UpdateResult("Advancer:RK4:DidNotCrashAndBurn",true);
216  serialUnitResults.UpdateResult("Advancer:RK4:Accurate",isAccurate);
217  serialUnitResults.UpdateResult("Advancer:RK4:Order4",isOrder4);
218 
219 };
220 
221 //#include "RK4Advancer2.H"
222 
223 void TestRK4Advancer2(ix::test::results &serialUnitResults){
224 
225  std::cout << "TestRK4Advancer2" << std::endl;
226 
227  // basic types for the advection problem
228  typedef simulation::domain::base<grid_t,state_t> domain_t;
229  typedef testrhs<domain_t> rhs_t;
230 
231  // -- Set up the RHS --
232  rhs_t testRHS;
233 
234  // -- Set up the grid --
235  grid_t testGrid;
236  operator_t sbpOperator;
237  state_t testState;
238 
239  std::ostringstream messageStream;
240  std::vector<size_t> gridSizes(3,1);
241  std::vector<size_t> numNodes(gridSizes);
242 
243  if(testfixtures::CreateSimpleGrid(testGrid,gridSizes)){
244  return;
245  }
246 
247  pcpp::IndexIntervalType globalExtent;
248  globalExtent.InitSimple(numNodes);
249 
250  size_t numNodesDomain = numNodes[0]*numNodes[1]*numNodes[2];
251  size_t numCellsDomain = (numNodes[0]-1)*(numNodes[1]-1)*(numNodes[2]-1);
252 
253  // set up the state
254  testState.AddField("A",'s',1,8,"");
255  testState.AddField("simTime",'s',1,8,"s");
256  testState.AddField("simDT",'s',1,8,"s");
257  testState.AddField("inputCFL",'s',1,8,"");
258  testState.AddField("F0",'n',1,8,"");
259  testState.AddField("F1",'n',1,8,"");
260  testState.AddField("F2",'n',1,8,"");
261  testState.AddField("F3",'n',1,8,"");
262  testState.AddField("F4",'n',1,8,"");
263  testState.AddField("F5",'n',1,8,"");
264  testState.AddField("F6",'n',1,8,"");
265  testState.Create(numNodesDomain,numCellsDomain);
266  testState.InitializeFieldHandles();
267  testState.SetStateFields("F0 F1 F2 F3 F4 F5 F6");
268 
269 
270  double a = 1.0;
271  double t = 1.0;
272  double f0 = a*testRHS.TestF(t,0);
273  double f1 = a*testRHS.TestF(t,1);
274  double f2 = a*testRHS.TestF(t,2);
275  double f3 = a*testRHS.TestF(t,3);
276  double f4 = a*testRHS.TestF(t,4);
277  double f5 = a*testRHS.TestF(t,5);
278  double f6 = a*testRHS.TestF(t,6);
279  double dt = 1.0;
280  double cfl = 1.0;
281 
282  // set up the domain
283  domain_t rkDomain;
284  rkDomain.SetNumGrids(1);
285  rkDomain.SetGrid(0,testGrid);
286  rkDomain.SetGridState(0,testState);
287  state_t &domainState(rkDomain.State(0));
288  rkDomain.SetRHS(testRHS);
289  rkDomain.PartitionDomain(0);
290  testRHS.SetTimestep(dt);
291 
292  domainState.SetFieldBuffer("A",&a);
293  domainState.SetFieldBuffer("simTime",&t);
294  domainState.SetFieldBuffer("simDT",&dt);
295  domainState.SetFieldBuffer("inputCFL",&cfl);
296  domainState.SetFieldBuffer("F0",&f0);
297  domainState.SetFieldBuffer("F1",&f1);
298  domainState.SetFieldBuffer("F2",&f2);
299  domainState.SetFieldBuffer("F3",&f3);
300  domainState.SetFieldBuffer("F4",&f4);
301  domainState.SetFieldBuffer("F5",&f5);
302  domainState.SetFieldBuffer("F6",&f6);
303 
304  // set up the advancer
305  rk4advancer<domain_t> testAdvancer;
306 
307  global_t myGlobal;
308 // testAdvancer.SetMessageStream(messageStream);
309 // testAdvancer.SetRHS(testRHS);
310  testAdvancer.InitializeAdvancer(rkDomain,myGlobal,messageStream);
311  // pcpp::ParallelGlobalType myGlobal;
312  // testAdvancer.SetGlobal(myGlobal);
313  // std::vector<pcpp::IndexIntervalType> &threadPartitionIntervals(testGrid.ThreadPartitionBufferIntervals());
314  // pcpp::IndexIntervalType &localThreadInterval(threadPartitionBufferIntervals[0]);
315  // pcpp::IndexIntervalType &localExtent(testGrid.LocalExtent());
316  // localThreadInterval = localExtent;
317  testAdvancer.SyncIntervals();
318 
319 // messageStream << "Domain State report: " << std::endl;
320 // domainState.Report(messageStream);
321 
322  std::cout << messageStream.str() << std::endl;
323  pcpp::io::RenewStream(messageStream);
324 
325  int nTrial = 10;
326  bool isAccurate = true;
327  bool isOrder4 = true;
328  double errorTolerance = 1e-15;
329 
330  std::vector<std::vector<double> > errorValues(nTrial);
331  std::vector<double> timeSteps;
332 
333  for(int iTrial = 0;iTrial < 10;iTrial++){
334 
335  // Reset all the input data for the current trial
336  t = 0.5;
337  double tbefore = t;
338  dt = 1.0/(std::pow(2.0,static_cast<double>(iTrial)));
339  testRHS.SetTimestep(dt);
340  f0 = a*testRHS.TestF(t,0);
341  f1 = a*testRHS.TestF(t,1);
342  f2 = a*testRHS.TestF(t,2);
343  f3 = a*testRHS.TestF(t,3);
344  f4 = a*testRHS.TestF(t,4);
345  f5 = a*testRHS.TestF(t,5);
346  f6 = a*testRHS.TestF(t,6);
347 
348  // Report to stdout some status/diagnostic
349  std::cout << "Driver:Time(before): " << t << std::endl;
350  std::cout << "Driver:TimeStep: " << dt << std::endl;
351 
352  // record the time step values
353  timeSteps.push_back(dt);
354  // domainState.Report(std::cout);
355  std::cout << "Driver:F(before): " << std::endl
356  << "Driver:f0 = " << f0 << std::endl
357  << "Driver:f1 = " << f1 << std::endl
358  << "Driver:f2 = " << f2 << std::endl
359  << "Driver:f3 = " << f3 << std::endl
360  << "Driver:f4 = " << f4 << std::endl
361  << "Driver:f5 = " << f5 << std::endl
362  << "Driver:f6 = " << f6 << std::endl;
363  // Advance one timestep
364  if(testAdvancer.AdvanceDomain()){
365  break;
366  }
367 
368  if(std::abs(t - tbefore - dt) >= 1e-24){
369  std::cout << "Time did not advance, aborting test." << std::endl;
370  return;
371  // Fail test and abort
372  }
373 
374  // sanity check that time was updated
375  // domainState.Report(std::cout);
376  std::cout << "Time(after): " << t << std::endl;
377  std::cout << "Driver:F(after): " << std::endl
378  << "Driver:f0 = " << f0 << std::endl
379  << "Driver:f1 = " << f1 << std::endl
380  << "Driver:f2 = " << f2 << std::endl
381  << "Driver:f3 = " << f3 << std::endl
382  << "Driver:f4 = " << f4 << std::endl
383  << "Driver:f5 = " << f5 << std::endl
384  << "Driver:f6 = " << f6 << std::endl;
385 
386 
387  // create storage for the error values and record them
388  errorValues[iTrial].resize(7,0.0);
389  errorValues[iTrial][0] = std::abs(f0 - a*testRHS.TestF(t,0));
390  errorValues[iTrial][1] = std::abs(f1 - a*testRHS.TestF(t,1));
391  errorValues[iTrial][2] = std::abs(f2 - a*testRHS.TestF(t,2));
392  errorValues[iTrial][3] = std::abs(f3 - a*testRHS.TestF(t,3));
393  errorValues[iTrial][4] = std::abs(f4 - a*testRHS.TestF(t,4));
394  errorValues[iTrial][5] = std::abs(f5 - a*testRHS.TestF(t,5));
395  errorValues[iTrial][6] = std::abs(f6 - a*testRHS.TestF(t,6));
396 
397  // Make sure the fields that should be integrated exactly
398  // have essentially zero error
399  for(int iVal = 0;iVal < 5;iVal++)
400  if(errorValues[iTrial][iVal] > errorTolerance)
401  isAccurate = false;
402 
403  // sanity check to stdout
404  for(int iVal = 0;iVal < 7;iVal++){
405  std::cout << "Driver:Error(" << iVal << ") = " << errorValues[iTrial][iVal] << std::endl;
406  }
407 
408  }
409 
410  // Check the truncation error is correctly O(dt^5)
411  for(int iTrial = 0;iTrial < 8;iTrial++){
412  double logErr = std::log(errorValues[iTrial+1][5]/errorValues[iTrial][5]);
413  double logDT = std::log(timeSteps[iTrial+1]/timeSteps[iTrial]);
414  double methodOrder = logErr/logDT;
415  std::cout << "iTrial(" << iTrial << ") Order: " << methodOrder << std::endl;
416  if(std::abs(methodOrder - 5) > 1e-3)
417  isOrder4 = false;
418  }
419 
420  serialUnitResults.UpdateResult("Advancer:RK4_2:DidNotCrashAndBurn",true);
421  serialUnitResults.UpdateResult("Advancer:RK4_2:Accurate",isAccurate);
422  serialUnitResults.UpdateResult("Advancer:RK4_2:Order4",isOrder4);
423 
424 };
425 
testfixtures::CreateSimpleGrid
int CreateSimpleGrid(GridType &inGrid, const std::vector< size_t > &gridSizes)
Definition: TestFixtures.H:10
global_t
pcpp::ParallelGlobalType global_t
Definition: TestRK4Advancer.C:16
simulation::state::base::SetStateFields
int SetStateFields(const std::vector< int > &inFieldIndices)
Definition: State.H:310
gridSizes
void const size_t const size_t * gridSizes
Definition: EulerKernels.H:10
ix::util::RenewStream
void RenewStream(std::ostringstream &outStream)
Definition: primitive_utilities.H:638
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
operator_t
plascom2::operators::sbp::base operator_t
Definition: TestRK4Advancer.C:14
rk4advancer::AdvanceDomain
virtual int AdvanceDomain()
Advances domain, grid, and state by one step.
Definition: RK4Advancer.H:455
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
ix::comm::CommunicatorObject
Main encapsulation of MPI.
Definition: COMM.H:62
a
void size_t int size_t int size_t int int int int double int int double double *void size_t int size_t int int int int int double int size_t size_t size_t double double *void size_t int size_t int size_t size_t int double int double double *void size_t size_t size_t double * a
Definition: OperatorKernels.H:29
Testing.H
Testing constructs for unit testing.
plascom2::operators::stencilset
Encapsulation for a collection of operator stencils.
Definition: Stencil.H:26
TestRK4Advancer2
void TestRK4Advancer2(ix::test::results &serialUnitResults)
Definition: TestRK4Advancer.C:223
grid_t
simulation::grid::parallel_blockstructured grid_t
Definition: TestRK4Advancer.C:12
state_t
simulation::state::base state_t
Definition: TestRK4Advancer.C:13
interval_t
pcpp::IndexIntervalType interval_t
Definition: TestRK4Advancer.C:18
ix::test::results::UpdateResult
void UpdateResult(const std::string &name, const ValueType &result)
Updates an existing test result.
Definition: Testing.H:55
rk4advancer::InitializeAdvancer
virtual int InitializeAdvancer(DomainType &inDomain, pcpp::ParallelGlobalType &inGlobal, std::ostream &inStream)
Initializes the advancer object.
Definition: RK4Advancer.H:168
TestRK4Advancer
void TestRK4Advancer(ix::test::results &serialUnitResults)
Definition: TestRK4Advancer.C:20
ix::util::sizeextent
Simple Block Structured Mesh object.
Definition: IndexUtil.H:21
pcpp::field::dataset::AddField
void AddField(const std::string &name, char loc, unsigned int ncomp, unsigned int dsize, const std::string &unit)
Definition: PCPPFieldData.H:869
pcpp::field::dataset::SetFieldBuffer
void SetFieldBuffer(const std::string &name, void *buf)
Definition: PCPPFieldData.H:716
ix::util::sizeextent::InitSimple
void InitSimple(const ContainerType &inSize)
Definition: IndexUtil.H:169
rk4advancer::SyncIntervals
void SyncIntervals()
Definition: RK4Advancer.H:896
simulation::grid::parallel_blockstructured
Definition: Grid.H:387