PlasCom2  1.0
XPACC Multi-physics simluation application
MaxwellUtil.H
Go to the documentation of this file.
1 #ifndef __MAXWELL_UTIL_H__
2 #define __MAXWELL_UTIL_H__
3 
4 
5 
6 namespace Maxwell {
7  namespace util {
8 
9  template<typename GridType, typename StateType>
10  int SetupMaxwellStateAndParameters(const GridType &inGrid, StateType &inState, StateType &inParams)
11  {
12 
13  size_t numPointsBuffer = inGrid.BufferSize();
14 
15  const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
16  int numDim = bufferSizes.size();
17 
18  if (numDim != 3) return(1);
19 
20  /**** State variables ****/
21  inState.AddField("simTime",'s',1,8,"s"); // time
22 
23  inState.AddField("Efield",'n',3,8,""); // E-field
24  inState.AddField("Bfield",'n',3,8,""); // B-field
25  inState.AddField("Dfield",'n',3,8,""); // D-field
26  inState.AddField("Hfield",'n',3,8,""); // H-field
27  inState.AddField("J",'n',3,8,""); // current density
28  inState.AddField("epsilonMu",'n',2,8,""); // permittivity and permeability
29 
30  inState.Create(numPointsBuffer, 0);
31  inState.SetStateFields("Efield Hfield");
32 
33  /**** Parameters ****/
34  inParams.AddField("epsilon0",'s',1,8,""); // permittivity of free space
35  inParams.AddField("mu0",'s',1,8,""); // permeability of free space
36 
37  inParams.AddField("inputCFL",'s',1,8,""); // CFL number
38  inParams.AddField("inputDT",'s',1,8,"s"); // timestep size
39  inParams.AddField("refC",'s',1,8,"speed"); // reference speed
40 
41  inParams.Create(numPointsBuffer, 0);
42 
43  return(0);
44 
45  };
46 
47 
48  template<typename GridType, typename StateType>
49  int InitializeMaxwellStateConstFields(const GridType &inGrid, StateType &inState,
50  double *constE, double *constB, double *constJ,
51  double &constEps, double &constMu,
52  bool everyWhere = false)
53  {
54 
55  // Get state buffers
56  int timeHandle = inState.GetDataIndex("simTime");
57  if(timeHandle < 0)
58  return(1);
59  pcpp::field::dataset::DataBufferType &timeData(inState.Field(timeHandle));
60  double *timePtr = timeData.Data<double>();
61  if(!timePtr)
62  return(1);
63 
64  int EfieldHandle = inState.GetDataIndex("Efield");
65  if(EfieldHandle < 0)
66  return(2);
67  pcpp::field::dataset::DataBufferType &EfieldData(inState.Field(EfieldHandle));
68  double *EfieldPtr = EfieldData.Data<double>();
69  if(!EfieldPtr)
70  return(2);
71 
72  int BfieldHandle = inState.GetDataIndex("Bfield");
73  if(BfieldHandle < 0)
74  return(3);
75  pcpp::field::dataset::DataBufferType &BfieldData(inState.Field(BfieldHandle));
76  double *BfieldPtr = BfieldData.Data<double>();
77  if(!BfieldPtr)
78  return(3);
79 
80  int DfieldHandle = inState.GetDataIndex("Dfield");
81  if(DfieldHandle < 0)
82  return(6);
83  pcpp::field::dataset::DataBufferType &DfieldData(inState.Field(DfieldHandle));
84  double *DfieldPtr = DfieldData.Data<double>();
85  if(!DfieldPtr)
86  return(6);
87 
88  int HfieldHandle = inState.GetDataIndex("Hfield");
89  if(HfieldHandle < 0)
90  return(7);
91  pcpp::field::dataset::DataBufferType &HfieldData(inState.Field(HfieldHandle));
92  double *HfieldPtr = HfieldData.Data<double>();
93  if(!HfieldPtr)
94  return(7);
95 
96  int JHandle = inState.GetDataIndex("J");
97  if(JHandle < 0)
98  return(4);
99  pcpp::field::dataset::DataBufferType &JData(inState.Field(JHandle));
100  double *JPtr = JData.Data<double>();
101  if(!JPtr)
102  return(4);
103 
104  int epsilonMuHandle = inState.GetDataIndex("epsilonMu");
105  if(epsilonMuHandle < 0)
106  return(5);
107  pcpp::field::dataset::DataBufferType &epsilonMuData(inState.Field(epsilonMuHandle));
108  double *epsilonMuPtr = epsilonMuData.Data<double>();
109  if(!epsilonMuPtr)
110  return(5);
111 
112  // Get grid, partition and buffer info.
113  size_t numPointsBuffer = inGrid.BufferSize();
114  const std::vector<size_t> &gridSizes(inGrid.GridSizes());
115  const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
116  const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
117  const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
118  int numDim = bufferSizes.size();
119 
120  if (numDim != 3) return(-1); // no. of dimensions must be equal to 3.
121 
122  // Initialize state variables
123  *timePtr = 0.0; // time
124 
125  std::vector<double> constD(3, 0.0);
126  constD[0]=constEps*constE[0];
127  constD[1]=constEps*constE[1];
128  constD[2]=constEps*constE[2];
129 
130  std::vector<double> constH(3, 0.0);
131  constH[0]=constB[0]/constMu;
132  constH[1]=constB[1]/constMu;
133  constH[2]=constB[2]/constMu;
134 
135  if(everyWhere) {
136 
137  // Do it everywhere
138  size_t twoxnumPointsBuffer = 2*numPointsBuffer;
139 
140  for(size_t iPoint=0; iPoint<numPointsBuffer; iPoint++) {
141  epsilonMuPtr[iPoint] = constEps; // epsilon
142  epsilonMuPtr[iPoint+numPointsBuffer] = constMu; // mu
143 
144  EfieldPtr[iPoint] = constE[0]; // E-field
145  EfieldPtr[iPoint+numPointsBuffer] = constE[1];
146  EfieldPtr[iPoint+twoxnumPointsBuffer] = constE[2];
147 
148  BfieldPtr[iPoint] = constB[0]; // B-field
149  BfieldPtr[iPoint+numPointsBuffer] = constB[1];
150  BfieldPtr[iPoint+twoxnumPointsBuffer] = constB[2];
151 
152  DfieldPtr[iPoint] = constD[0]; // D-field
153  DfieldPtr[iPoint+numPointsBuffer] = constD[1];
154  DfieldPtr[iPoint+twoxnumPointsBuffer] = constD[2];
155 
156  HfieldPtr[iPoint] = constH[0]; // H-field
157  HfieldPtr[iPoint+numPointsBuffer] = constH[1];
158  HfieldPtr[iPoint+twoxnumPointsBuffer] = constH[2];
159 
160  JPtr[iPoint] = constJ[0]; // current density
161  JPtr[iPoint+numPointsBuffer] = constJ[1];
162  JPtr[iPoint+twoxnumPointsBuffer] = constJ[2];
163  }
164 
165  } else {
166 
167  // Do it only in the local partition
168  size_t iStart = partitionBufferInterval[0].first;
169  size_t iEnd = partitionBufferInterval[0].second;
170  size_t jStart = partitionBufferInterval[1].first;
171  size_t jEnd = partitionBufferInterval[1].second;
172  size_t kStart = partitionBufferInterval[2].first;
173  size_t kEnd = partitionBufferInterval[2].second;
174 
175  size_t nPlane = bufferSizes[0]*bufferSizes[1];
176  size_t twoxnumPointsBuffer = 2*numPointsBuffer;
177 
178  for(size_t kIndex=kStart; kIndex<=kEnd; kIndex++) {
179  size_t kBufferIndex = kIndex*nPlane;
180  for(size_t jIndex=jStart; jIndex<=jEnd; jIndex++) {
181  size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
182  for(size_t iIndex=iStart; iIndex<=iEnd; iIndex++) {
183  size_t bufferIndex = jkBufferIndex + iIndex;
184 
185  epsilonMuPtr[bufferIndex] = constEps; // epsilon
186  epsilonMuPtr[bufferIndex+numPointsBuffer] = constMu; // mu
187 
188  EfieldPtr[bufferIndex] = constE[0]; // E-field
189  EfieldPtr[bufferIndex+numPointsBuffer] = constE[1];
190  EfieldPtr[bufferIndex+twoxnumPointsBuffer] = constE[2];
191 
192  BfieldPtr[bufferIndex] = constB[0]; // B-field
193  BfieldPtr[bufferIndex+numPointsBuffer] = constB[1];
194  BfieldPtr[bufferIndex+twoxnumPointsBuffer] = constB[2];
195 
196  DfieldPtr[bufferIndex] = constD[0]; // D-field
197  DfieldPtr[bufferIndex+numPointsBuffer] = constD[1];
198  DfieldPtr[bufferIndex+twoxnumPointsBuffer] = constD[2];
199 
200  HfieldPtr[bufferIndex] = constH[0]; // H-field
201  HfieldPtr[bufferIndex+numPointsBuffer] = constH[1];
202  HfieldPtr[bufferIndex+twoxnumPointsBuffer] = constH[2];
203 
204  JPtr[bufferIndex] = constJ[0]; // current density
205  JPtr[bufferIndex+numPointsBuffer] = constJ[1];
206  JPtr[bufferIndex+twoxnumPointsBuffer] = constJ[2];
207  }
208  }
209  }
210 
211  }
212 
213  return(0);
214 
215  };
216 
217 
218  template<typename GridType, typename StateType>
219  int InitializeMaxwellStateGridIndices(const GridType &inGrid, StateType &inState)
220  {
221 
222  // Get state buffers
223  int timeHandle = inState.GetDataIndex("simTime");
224  if(timeHandle < 0)
225  return(1);
226  pcpp::field::dataset::DataBufferType &timeData(inState.Field(timeHandle));
227  double *timePtr = timeData.Data<double>();
228  if(!timePtr)
229  return(1);
230 
231  int EfieldHandle = inState.GetDataIndex("Efield");
232  if(EfieldHandle < 0)
233  return(2);
234  pcpp::field::dataset::DataBufferType &EfieldData(inState.Field(EfieldHandle));
235  double *EfieldPtr = EfieldData.Data<double>();
236  if(!EfieldPtr)
237  return(2);
238 
239  int BfieldHandle = inState.GetDataIndex("Bfield");
240  if(BfieldHandle < 0)
241  return(3);
242  pcpp::field::dataset::DataBufferType &BfieldData(inState.Field(BfieldHandle));
243  double *BfieldPtr = BfieldData.Data<double>();
244  if(!BfieldPtr)
245  return(3);
246 
247  int DfieldHandle = inState.GetDataIndex("Dfield");
248  if(DfieldHandle < 0)
249  return(6);
250  pcpp::field::dataset::DataBufferType &DfieldData(inState.Field(DfieldHandle));
251  double *DfieldPtr = DfieldData.Data<double>();
252  if(!DfieldPtr)
253  return(6);
254 
255  int HfieldHandle = inState.GetDataIndex("Hfield");
256  if(HfieldHandle < 0)
257  return(7);
258  pcpp::field::dataset::DataBufferType &HfieldData(inState.Field(HfieldHandle));
259  double *HfieldPtr = HfieldData.Data<double>();
260  if(!HfieldPtr)
261  return(7);
262 
263  int JHandle = inState.GetDataIndex("J");
264  if(JHandle < 0)
265  return(4);
266  pcpp::field::dataset::DataBufferType &JData(inState.Field(JHandle));
267  double *JPtr = JData.Data<double>();
268  if(!JPtr)
269  return(4);
270 
271  int epsilonMuHandle = inState.GetDataIndex("epsilonMu");
272  if(epsilonMuHandle < 0)
273  return(5);
274  pcpp::field::dataset::DataBufferType &epsilonMuData(inState.Field(epsilonMuHandle));
275  double *epsilonMuPtr = epsilonMuData.Data<double>();
276  if(!epsilonMuPtr)
277  return(5);
278 
279  // Get grid, partition and buffer info.
280  size_t numPointsBuffer = inGrid.BufferSize();
281  const std::vector<size_t> &gridSizes(inGrid.GridSizes());
282  const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
283  const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
284  const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
285  int numDim = bufferSizes.size();
286 
287  if (numDim != 3) return(-1); // no. of dimensions must be equal to 3.
288 
289  // Initialize state variables
290  *timePtr = 0.0; // time
291 
292  // Do it only in the local partition
293  size_t iStart = partitionBufferInterval[0].first;
294  size_t iEnd = partitionBufferInterval[0].second;
295  size_t jStart = partitionBufferInterval[1].first;
296  size_t jEnd = partitionBufferInterval[1].second;
297  size_t kStart = partitionBufferInterval[2].first;
298  size_t kEnd = partitionBufferInterval[2].second;
299 
300  size_t nPlane = bufferSizes[0]*bufferSizes[1];
301  size_t twoxnumPointsBuffer = 2*numPointsBuffer;
302 
303  size_t iStartG = partitionInterval[0].first;
304  size_t iEndG = partitionInterval[0].second;
305  size_t jStartG = partitionInterval[1].first;
306  size_t jEndG = partitionInterval[1].second;
307  size_t kStartG = partitionInterval[2].first;
308  size_t kEndG = partitionInterval[2].second;
309 
310  size_t nPlaneG = gridSizes[0]*gridSizes[1];
311 
312  for(size_t kIndex=kStart; kIndex<=kEnd; kIndex++) {
313  size_t kBufferIndex = kIndex*nPlane;
314  size_t kGridIndex = (kStartG+kIndex-kStart)*nPlaneG;
315  for(size_t jIndex=jStart; jIndex<=jEnd; jIndex++) {
316  size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
317  size_t jkGridIndex = kGridIndex + (jStartG+jIndex-jStart)*gridSizes[0];
318  for(size_t iIndex=iStart; iIndex<=iEnd; iIndex++) {
319  size_t bufferIndex = jkBufferIndex + iIndex;
320  size_t gridIndex = jkGridIndex + (iStartG+iIndex-iStart+1);
321 
322  epsilonMuPtr[bufferIndex] = gridIndex; // epsilon
323  epsilonMuPtr[bufferIndex+numPointsBuffer] = gridIndex; // mu
324 
325  EfieldPtr[bufferIndex] = gridIndex; // E-field
326  EfieldPtr[bufferIndex+numPointsBuffer] = gridIndex;
327  EfieldPtr[bufferIndex+twoxnumPointsBuffer] = gridIndex;
328 
329  BfieldPtr[bufferIndex] = gridIndex; // B-field
330  BfieldPtr[bufferIndex+numPointsBuffer] = gridIndex;
331  BfieldPtr[bufferIndex+twoxnumPointsBuffer] = gridIndex;
332 
333  DfieldPtr[bufferIndex] = gridIndex*gridIndex; // D-field
334  DfieldPtr[bufferIndex+numPointsBuffer] = gridIndex*gridIndex;
335  DfieldPtr[bufferIndex+twoxnumPointsBuffer] = gridIndex*gridIndex;
336 
337  HfieldPtr[bufferIndex] = 1.0; // H-field
338  HfieldPtr[bufferIndex+numPointsBuffer] = 1.0;
339  HfieldPtr[bufferIndex+twoxnumPointsBuffer] = 1.0;
340 
341  JPtr[bufferIndex] = gridIndex; // current density
342  JPtr[bufferIndex+numPointsBuffer] = gridIndex;
343  JPtr[bufferIndex+twoxnumPointsBuffer] = gridIndex;
344  }
345  }
346  }
347 
348  return(0);
349 
350  };
351 
352 
353  template<typename GridType, typename StateType>
354  int InitializeMaxwellParameters(const GridType &inGrid, StateType &inParams, double CFLValue)
355  {
356 
357  // Get parameters buffers
358  int eps0Handle = inParams.GetDataIndex("epsilon0");
359  if(eps0Handle < 0)
360  return(1);
361  pcpp::field::dataset::DataBufferType &eps0Data(inParams.Field(eps0Handle));
362  double *eps0Ptr = eps0Data.Data<double>();
363  if(!eps0Ptr)
364  return(1);
365 
366  int mu0Handle = inParams.GetDataIndex("mu0");
367  if(mu0Handle < 0)
368  return(2);
369  pcpp::field::dataset::DataBufferType &mu0Data(inParams.Field(mu0Handle));
370  double *mu0Ptr = mu0Data.Data<double>();
371  if(!mu0Ptr)
372  return(2);
373 
374  int CFLHandle = inParams.GetDataIndex("inputCFL");
375  if(CFLHandle < 0)
376  return(3);
377  pcpp::field::dataset::DataBufferType &CFLData(inParams.Field(CFLHandle));
378  double *CFLPtr = CFLData.Data<double>();
379  if(!CFLPtr)
380  return(3);
381 
382  int deltaTHandle = inParams.GetDataIndex("inputDT");
383  if(deltaTHandle < 0)
384  return(4);
385  pcpp::field::dataset::DataBufferType &deltaTData(inParams.Field(deltaTHandle));
386  double *deltaTPtr = deltaTData.Data<double>();
387  if(!deltaTPtr)
388  return(4);
389 
390  int cRefHandle = inParams.GetDataIndex("refC");
391  if(cRefHandle < 0)
392  return(5);
393  pcpp::field::dataset::DataBufferType &cRefData(inParams.Field(cRefHandle));
394  double *cRefPtr = cRefData.Data<double>();
395  if(!cRefPtr)
396  return(5);
397 
398  // Initialize parameters
399  const std::vector<double> &dX(inGrid.DX()); // get grid size
400  const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
401  int numDim = bufferSizes.size();
402 
403  if (numDim != 3) return(-1);
404 
405  std::vector<double>::const_iterator dxIt = dX.begin();
406  double minSpacing = *dxIt;
407  while(dxIt != dX.end()) {
408  if(*dxIt < minSpacing) minSpacing = *dxIt;
409  dxIt++;
410  }
411 
412  *eps0Ptr = 8.85418782e-12; // permittivity of free space [F/m]
413  *mu0Ptr = 1.25663706e-6; // permeability of free space [H/m]
414 
415  *CFLPtr = CFLValue; // CFL number
416  *cRefPtr = 299792458.131; // speed of light in free space [m/s]
417  *deltaTPtr = (*CFLPtr)*minSpacing/(*cRefPtr); // timestep size [s]
418 
419  return(0);
420 
421  };
422 
423 
424 // template<typename StateType, typename HaloType>
425 // int CommunicateMaxwellHalo(StateType &inState, HaloType &inHalo)
426 // {
427 //
428 // // Get state buffers
429 // double *EfieldPtr = inState.template GetFieldBuffer<double>("Efield");
430 // double *BfieldPtr = inState.template GetFieldBuffer<double>("Bfield");
431 // double *JPtr = inState.template GetFieldBuffer<double>("J");
432 // double *epsilonMuPtr = inState.template GetFieldBuffer<double>("epsilonMu");
433 //
434 // // Allocate memory for halo buffers
435 // int bufferSize = 11;
436 // int stateBufferMessageID = inHalo.CreateMessageBuffers(bufferSize);
437 //
438 // // Pack message buffers
439 // inHalo.PackMessageBuffers(stateBufferMessageID,0,3,&EfieldPtr[0]);
440 // inHalo.PackMessageBuffers(stateBufferMessageID,3,3,&BfieldPtr[0]);
441 // inHalo.PackMessageBuffers(stateBufferMessageID,6,3,&JPtr[0]);
442 // inHalo.PackMessageBuffers(stateBufferMessageID,9,2,&epsilonMuPtr[0]);
443 //
444 //
445 //
446 //
447 // // Unpack message buffers
448 // inHalo.UnPackMessageBuffers(stateBufferMessageID,0,3,&EfieldPtr[0]);
449 // inHalo.UnPackMessageBuffers(stateBufferMessageID,3,3,&BfieldPtr[0]);
450 // inHalo.UnPackMessageBuffers(stateBufferMessageID,6,3,&JPtr[0]);
451 // inHalo.UnPackMessageBuffers(stateBufferMessageID,9,2,&epsilonMuPtr[0]);
452 //
453 // return(0);
454 //
455 // };
456 
457 
458  template<typename GridType, typename StateType>
459  int InitializeGaussianPulse1DXDir(const GridType &inGrid, StateType &inState, StateType &inParams,
460  const std::vector<double> &inputParams,
461  const std::vector<int> &inputFlags,
462  int threadID)
463  {
464 
465  /* ==== Set pulse parameters ==== */
466  double pulseAmp = 1.0;
467  double pulseWidth = 0.01;
468  double pulseCutOffWidth = 1.0;
469  double pulseCenterX = 0.0;
470 
471  int numInputParams = inputParams.size();
472 
473  if(numInputParams > 0)
474  pulseAmp = inputParams[0];
475  if(numInputParams > 1)
476  pulseWidth = inputParams[1];
477  if(numInputParams > 2)
478  pulseCutOffWidth = inputParams[2];
479  if(numInputParams > 3)
480  pulseCenterX = inputParams[3];
481 
482  /* ==== Get grid info. ==== */
483  size_t numPointsBuffer = inGrid.BufferSize();
484  const std::vector<size_t> &bufferSizes(inGrid.BufferSizes());
485  const std::vector<size_t> &gridSizes(inGrid.GridSizes());
486  const pcpp::IndexIntervalType &partitionInterval(inGrid.PartitionInterval());
487  const pcpp::IndexIntervalType &partitionBufferInterval(inGrid.PartitionBufferInterval());
488  const std::vector<simulation::grid::subregion> &gridSubRegions(inGrid.SubRegions());
489  const std::vector<double> &dX(inGrid.DX());
490  const std::vector<double> &coordinateData(inGrid.CoordinateData());
491  const std::vector<double> &physicalExtent(inGrid.PhysicalExtent());
492  int numDim = gridSizes.size();
493 
494  pcpp::IndexIntervalType gridInterval;
495  gridInterval.InitSimple(gridSizes);
496 
497  std::cout << "GaussianPulse1DXDir: Grid Sizes: (";
498  pcpp::io::DumpContents(std::cout,gridSizes,",");
499  std::cout << ")" << std::endl
500  << "GaussianPulse1DXDir: Partition Interval: ";
501  partitionInterval.PrettyPrint(std::cout);
502  std::cout << std::endl << "GaussianPulse1DXDir: Partition Buffer Interval: ";
503  partitionBufferInterval.PrettyPrint(std::cout);
504  std::cout << std::endl;
505 
506  std::cout << "GaussianPulse1DXDir: Amplitude: "
507  << pulseAmp << std::endl
508  << "GaussianPulse1DXDir: Width: "
509  << pulseWidth << std::endl
510  << "GaussianPulse1DXDir: CutOffWidth: "
511  << pulseCutOffWidth << std::endl
512  << "GaussianPulse1DXDir: CenterX: "
513  << pulseCenterX << std::endl;
514 
515  /* ==== Get state buffers ==== */
516  int timeHandle = inState.GetDataIndex("simTime");
517  if(timeHandle < 0)
518  return(1);
519  pcpp::field::dataset::DataBufferType &timeData(inState.Field(timeHandle));
520  double *timePtr = timeData.Data<double>();
521  if(!timePtr)
522  return(1);
523 
524  int EfieldHandle = inState.GetDataIndex("Efield");
525  if(EfieldHandle < 0)
526  return(2);
527  pcpp::field::dataset::DataBufferType &EfieldData(inState.Field(EfieldHandle));
528  double *EfieldPtr = EfieldData.Data<double>();
529  if(!EfieldPtr)
530  return(2);
531 
532  int BfieldHandle = inState.GetDataIndex("Bfield");
533  if(BfieldHandle < 0)
534  return(3);
535  pcpp::field::dataset::DataBufferType &BfieldData(inState.Field(BfieldHandle));
536  double *BfieldPtr = BfieldData.Data<double>();
537  if(!BfieldPtr)
538  return(3);
539 
540  int DfieldHandle = inState.GetDataIndex("Dfield");
541  if(DfieldHandle < 0)
542  return(6);
543  pcpp::field::dataset::DataBufferType &DfieldData(inState.Field(DfieldHandle));
544  double *DfieldPtr = DfieldData.Data<double>();
545  if(!DfieldPtr)
546  return(6);
547 
548  int HfieldHandle = inState.GetDataIndex("Hfield");
549  if(HfieldHandle < 0)
550  return(7);
551  pcpp::field::dataset::DataBufferType &HfieldData(inState.Field(HfieldHandle));
552  double *HfieldPtr = HfieldData.Data<double>();
553  if(!HfieldPtr)
554  return(7);
555 
556  int JHandle = inState.GetDataIndex("J");
557  if(JHandle < 0)
558  return(4);
559  pcpp::field::dataset::DataBufferType &JData(inState.Field(JHandle));
560  double *JPtr = JData.Data<double>();
561  if(!JPtr)
562  return(4);
563 
564  int epsilonMuHandle = inState.GetDataIndex("epsilonMu");
565  if(epsilonMuHandle < 0)
566  return(5);
567  pcpp::field::dataset::DataBufferType &epsilonMuData(inState.Field(epsilonMuHandle));
568  double *epsilonMuPtr = epsilonMuData.Data<double>();
569  if(!epsilonMuPtr)
570  return(5);
571 
572  /* ==== Get parameters buffers ==== */
573  int eps0Handle = inParams.GetDataIndex("epsilon0");
574  if(eps0Handle < 0)
575  return(8);
576  pcpp::field::dataset::DataBufferType &eps0Data(inParams.Field(eps0Handle));
577  double *eps0Ptr = eps0Data.Data<double>();
578  if(!eps0Ptr)
579  return(8);
580 
581  int mu0Handle = inParams.GetDataIndex("mu0");
582  if(mu0Handle < 0)
583  return(9);
584  pcpp::field::dataset::DataBufferType &mu0Data(inParams.Field(mu0Handle));
585  double *mu0Ptr = mu0Data.Data<double>();
586  if(!mu0Ptr)
587  return(9);
588 
589  int CFLHandle = inParams.GetDataIndex("inputCFL");
590  if(CFLHandle < 0)
591  return(10);
592  pcpp::field::dataset::DataBufferType &CFLData(inParams.Field(CFLHandle));
593  double *CFLPtr = CFLData.Data<double>();
594  if(!CFLPtr)
595  return(10);
596 
597  int deltaTHandle = inParams.GetDataIndex("inputDT");
598  if(deltaTHandle < 0)
599  return(11);
600  pcpp::field::dataset::DataBufferType &deltaTData(inParams.Field(deltaTHandle));
601  double *deltaTPtr = deltaTData.Data<double>();
602  if(!deltaTPtr)
603  return(11);
604 
605  int cRefHandle = inParams.GetDataIndex("refC");
606  if(cRefHandle < 0)
607  return(12);
608  pcpp::field::dataset::DataBufferType &cRefData(inParams.Field(cRefHandle));
609  double *cRefPtr = cRefData.Data<double>();
610  if(!cRefPtr)
611  return(12);
612 
613  /* ==== Initialize fields ==== */
614 
615  // Check no. of dimensions
616  if(numDim != 3) return(-1);
617 
618  double recipCRef = 1.0/(*cRefPtr);
619 
620  size_t iStart = partitionBufferInterval[0].first;
621  size_t iEnd = partitionBufferInterval[0].second;
622  size_t jStart = partitionBufferInterval[1].first;
623  size_t jEnd = partitionBufferInterval[1].second;
624  size_t kStart = partitionBufferInterval[2].first;
625  size_t kEnd = partitionBufferInterval[2].second;
626  size_t nPlane = bufferSizes[0]*bufferSizes[1];
627  for(size_t kIndex = kStart; kIndex <= kEnd; kIndex++) {
628  size_t kBufferIndex = kIndex*nPlane;
629  for(size_t jIndex = jStart; jIndex <= jEnd; jIndex++) {
630  size_t jkBufferIndex = kBufferIndex + jIndex*bufferSizes[0];
631  for(size_t iIndex = iStart; iIndex <= iEnd; iIndex++) {
632  size_t bufferIndex = jkBufferIndex + iIndex;
633 
634  size_t iPartIndex = (iIndex - iStart) + partitionInterval[0].first;
635  double x = iPartIndex*dX[0];
636  double xRad = x - pulseCenterX;
637 
638  if(std::fabs(xRad) <= pulseCutOffWidth) {
639  // Y-component of E-field
640  EfieldPtr[bufferIndex+numPointsBuffer] += euler::util::Pulse(pulseAmp, pulseWidth, pulseCenterX, 1.0, 1.0, x, 1.0, 1.0);
641 
642  // Z-component of B-field
643  BfieldPtr[bufferIndex+2*numPointsBuffer] += recipCRef*euler::util::Pulse(pulseAmp, pulseWidth, pulseCenterX, 1.0, 1.0, x, 1.0, 1.0);
644 
645  // Y-component of D-field
646  DfieldPtr[bufferIndex+numPointsBuffer] += epsilonMuPtr[bufferIndex]*EfieldPtr[bufferIndex+numPointsBuffer];
647 
648  // Z-component of H-field
649  HfieldPtr[bufferIndex+2*numPointsBuffer] += BfieldPtr[bufferIndex+2*numPointsBuffer]/epsilonMuPtr[bufferIndex+numPointsBuffer];
650  }
651  }
652  }
653  }
654 
655  return(0);
656 
657  };
658 
659  }
660 }
661 #endif
Maxwell::util::InitializeMaxwellParameters
int InitializeMaxwellParameters(const GridType &inGrid, StateType &inParams, double CFLValue)
Definition: MaxwellUtil.H:354
Maxwell::util::InitializeMaxwellStateConstFields
int InitializeMaxwellStateConstFields(const GridType &inGrid, StateType &inState, double *constE, double *constB, double *constJ, double &constEps, double &constMu, bool everyWhere=false)
Definition: MaxwellUtil.H:49
Maxwell::util::InitializeGaussianPulse1DXDir
int InitializeGaussianPulse1DXDir(const GridType &inGrid, StateType &inState, StateType &inParams, const std::vector< double > &inputParams, const std::vector< int > &inputFlags, int threadID)
Definition: MaxwellUtil.H:459
gridSizes
void const size_t const size_t * gridSizes
Definition: EulerKernels.H:10
x
void const size_t const size_t const size_t const double const double * x
Definition: OperatorKernels.H:33
Maxwell::util::InitializeMaxwellStateGridIndices
int InitializeMaxwellStateGridIndices(const GridType &inGrid, StateType &inState)
Definition: MaxwellUtil.H:219
Maxwell::util::SetupMaxwellStateAndParameters
int SetupMaxwellStateAndParameters(const GridType &inGrid, StateType &inState, StateType &inParams)
Definition: MaxwellUtil.H:10
apptools::DumpContents
void DumpContents(std::ostream &Ostr, const ContainerType &c, std::string del="\)
Dump container contents.
Definition: AppTools.H:117
bufferSizes
void const size_t const size_t * bufferSizes
Definition: MetricKernels.H:19
ix::util::sizeextent
Simple Block Structured Mesh object.
Definition: IndexUtil.H:21
ix::util::sizeextent::InitSimple
void InitSimple(const ContainerType &inSize)
Definition: IndexUtil.H:169
numPointsBuffer
void const size_t * numPointsBuffer
Definition: MetricKernels.H:19
euler::util::Pulse
double Pulse(double amp, double width, double centerX, double centerY, double centerZ, double x, double y, double z)
Definition: EulerUtil.H:56