-
Notifications
You must be signed in to change notification settings - Fork 0
/
World.cpp
100 lines (84 loc) · 2.44 KB
/
World.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
/*defines the simulation domain*/
#include <random>
#include <math.h>
#include "World.h"
#include "Field.h"
#include "Species.h"
//make an instance of the Rnd class
Rnd rnd;
using namespace std;
/*constructor*/
World::World(int ni, int nj, int nk):
ni{ni}, nj{nj}, nk{nk}, nn{ni,nj,nk},
phi(ni,nj,nk),rho(ni,nj,nk),node_vol(ni,nk,nk),
ef(ni,nj,nk), object_id(ni,nj,nk) {
time_start = chrono::high_resolution_clock::now(); //save starting time point
}
/*sets domain bounding box and computes mesh spacing*/
void World::setExtents(double3 _x0, double3 _xm) {
/*set origin and the opposite corner*/
x0 = _x0;
xm = _xm;
/*compute spacing by dividing length by the number of cells*/
for (int i=0;i<3;i++)
dh[i] = (xm(i)-x0(i))/(nn(i)-1);
//compute centroid
xc = 0.5*(x0+xm);
/*recompute node volumes*/
computeNodeVolumes();
}
/*marks box boundaries as 0V Dirichlet boundary*/
void World::setBoundaries() {
// Label the boundaries
for (int i=0;i<ni;i++)
for (int j=0;j<nj;j++)
for (int k=0;k<nk;k++) {
if ((i == 0) || (i == ni-1) || (j == 0) || (j == nj-1) || (k == 0) || (k == nk-1)) {
object_id[i][j][k] = 1; // DIRICHLET boundary on box
phi[i][j][k] = 0;
}
}
}
/*returns elapsed wall time in seconds*/
double World::getWallTime() {
auto time_now = chrono::high_resolution_clock::now();
chrono::duration<double> time_delta = time_now-time_start;
return time_delta.count();
}
/*computes charge density from rho = sum(charge*den)*/
void World::computeChargeDensity(vector<Species> &species)
{
rho = 0;
for (Species &sp:species)
{
if (sp.charge==0) continue; //don't bother with neutrals
rho += sp.charge*sp.den;
}
}
/*computes node volumes, dx*dy*dz on internal nodes and fractional
* values on domain boundary faces*/
void World::computeNodeVolumes() {
for (int i=0;i<ni;i++)
for (int j=0;j<nj;j++)
for (int k=0;k<nk;k++)
{
double V = dh[0]*dh[1]*dh[2]; //default volume
if (i==0 || i==ni-1) V*=0.5; //reduce by two for each boundary index
if (j==0 || j==nj-1) V*=0.5;
if (k==0 || k==nk-1) V*=0.5;
node_vol[i][j][k] = V;
}
}
/* computes total potential energy from 0.5*eps0*sum(E^2)*/
double World::getPE() {
double pe = 0;
for (int i=0;i<ni;i++)
for (int j=0;j<nj;j++)
for (int k=0;k<nk;k++)
{
double3 efn = ef[i][j][k]; //ef at this node
double ef2 = efn[0]*efn[0]+efn[1]*efn[1]+efn[2]*efn[2];
pe += ef2*node_vol[i][j][k];
}
return 0.5*Const::EPS_0*pe;
}