Execute post-processing step on a sublattice.
85 {
86
87 int newX0, newX1, newY0, newY1;
89 x0, x1, y0, y1,
90 x0_, x1_, y0_, y1_,
91 newX0, newX1, newY0, newY1 ) ) {
92
93 for (int iX=newX0; iX<=newX1; ++iX) {
94 for (int iY=newY0; iY<=newY1; ++iY) {
95 T materialNumber = blockLattice.get(iX,iY).template getField<descriptors::CELL_TYPE>();
96
97
98
99
100
101 T physDiffusion_c;
102 T tmp_source_c;
103 T physDiffusion_phi;
104 T tmp_source_phi;
105
106 std::vector<T> conc;
107 std::vector<T> conc_plusX;
108 std::vector<T> conc_minusX;
109 std::vector<T> conc_plusY;
110 std::vector<T> conc_minusY;
111
112 conc.emplace_back(blockLattice.get(iX,iY).computeRho());
113 for (int iter_component = 0; iter_component<component_number-1; ++ iter_component){
114 conc.emplace_back(tpartners[iter_component]->
get(iX,iY).computeRho());
115 }
116
117 T tmp_Diffusion_E = 1.2e-21*
util::pow(conc[0],4)- 6.5e-18*
util::pow(conc[0],3)+ 1.14e-14*
util::pow(conc[0],2)-8.06e-12*conc[0]+2.24e-9;
119
120 if ((int) materialNumber == 6 || (int) materialNumber == 7 || (int) materialNumber == 1 || (int) materialNumber == 3 || (int) materialNumber == 4 || (int) materialNumber == 5){
121
122 if (iX== newX0 ){
123 conc_minusX.emplace_back(0.0);
124 conc_minusX.emplace_back(0.0);
125 }
126 else {
127 conc_minusX.emplace_back(blockLattice.get(iX-1,iY).computeRho());
128 conc_minusX.emplace_back(tpartners[0]->
get(iX-1,iY).computeRho());
129 }
130
131 if (iX== newX1){
132 conc_plusX.emplace_back(0.0);
133 conc_plusX.emplace_back(0.0);
134 }
135 else {
136 conc_plusX.emplace_back(blockLattice.get(iX+1,iY).computeRho());
137 conc_plusX.emplace_back(tpartners[0]->
get(iX+1,iY).computeRho());
138 }
139
140 if (iY== newY0){
141 conc_minusY.emplace_back(blockLattice.get(iX, newY1).computeRho());
142 conc_minusY.emplace_back(tpartners[0]->
get(iX, newY1).computeRho());
143 }
144 else {
145 conc_minusY.emplace_back(blockLattice.get(iX,iY-1).computeRho());
146 conc_minusY.emplace_back(tpartners[0]->
get(iX,iY-1).computeRho());
147 }
148
149 if (iY == newY1){
150 conc_plusY.emplace_back(blockLattice.get(iX,newY0).computeRho());
151 conc_plusY.emplace_back(tpartners[0]->
get(iX,newY0).computeRho());
152 }
153 else {
154 conc_plusY.emplace_back(blockLattice.get(iX,iY+1).computeRho());
155 conc_plusY.emplace_back(tpartners[0]->
get(iX,iY+1).computeRho());
156 }
157
158 T kappa_E_minusX = -2.39e-11*
util::pow(conc_minusX[0],4)+1.21e-7*
util::pow(conc_minusX[0],3)-2.89e-4*
util::pow(conc_minusX[0],2)+0.32*conc_minusX[0]-2.789;
159 T kappa_E_plusX = -2.39e-11*
util::pow(conc_plusX[0],4)+1.21e-7*
util::pow(conc_plusX[0],3)-2.89e-4*
util::pow(conc_plusX[0],2)+0.32*conc_plusX[0]-2.789;
160 T kappa_E_minusY = -2.39e-11*
util::pow(conc_minusY[0],4)+1.21e-7*
util::pow(conc_minusY[0],3)-2.89e-4*
util::pow(conc_minusY[0],2)+0.32*conc_minusY[0]-2.789;
161 T kappa_E_plusY = -2.39e-11*
util::pow(conc_plusY[0],4)+1.21e-7*
util::pow(conc_plusY[0],3)-2.89e-4*
util::pow(conc_plusY[0],2)+0.32*conc_plusY[0]-2.789;
162
163 if ((int) materialNumber == 6 || (int)materialNumber == 5){
164 physDiffusion_c = diffusion_S;
165 physDiffusion_phi = kappa_S;
166 tmp_source_c = 0.0;
167 tmp_source_phi= 0.0;
168 }
169 else if ((int) materialNumber == 7 || (int) materialNumber == 3){
170 physDiffusion_c = fac_sep*(tmp_Diffusion_E- factor1*factor2*(kappa_E/conc[0]));
171 physDiffusion_phi = kappa_E;
172 tmp_source_c = fac_sep*(factor1* (1./4.)*((kappa_E_plusX-kappa_E_minusX)*(conc_plusX[1]-conc_minusX[1])+
173 + (kappa_E_plusY-kappa_E_minusY)*(conc_plusY[1]-conc_minusY[1]))
174 -kappa_E*(conc_minusX[1]+conc_plusX[1]+ conc_minusY[1]+conc_plusY[1]-4*conc[1]));
175 tmp_source_phi = factor2*((1./4.)*((kappa_E_plusX/conc_plusX[0]-kappa_E_minusX/conc_minusX[0])*(conc_plusX[0]-conc_minusX[0])
176 + (kappa_E_plusY/conc_plusY[0]-kappa_E_minusY/conc_minusY[0])*(conc_plusY[0]-conc_minusY[0]))
177 -(kappa_E/conc[0])*(conc_minusX[0]+conc_plusX[0]+conc_minusY[0]+conc_plusY[0]-4*conc[0]));
178 }
179 else if((int) materialNumber == 1 || (int) materialNumber == 4) {
180 physDiffusion_c = tmp_Diffusion_E- factor1*factor2*(kappa_E/conc[0]);
181 physDiffusion_phi = kappa_E;
182 tmp_source_c = (factor1* (1./4.)*((kappa_E_plusX-kappa_E_minusX)*(conc_plusX[1]-conc_minusX[1])+
183 + (kappa_E_plusY-kappa_E_minusY)*(conc_plusY[1]-conc_minusY[1]))
184 -kappa_E*(conc_minusX[1]+conc_plusX[1]+ conc_minusY[1]+conc_plusY[1]-4*conc[1]));
185 tmp_source_phi = factor2*((1./4.)*((kappa_E_plusX/conc_plusX[0]-kappa_E_minusX/conc_minusX[0])*(conc_plusX[0]-conc_minusX[0])
186 + (kappa_E_plusY/conc_plusY[0]-kappa_E_minusY/conc_minusY[0])*(conc_plusY[0]-conc_minusY[0]))
187 -(kappa_E/conc[0])*(conc_minusX[0]+conc_plusX[0]+conc_minusY[0]+conc_plusY[0]-4*conc[0]));
188
189 }
190
191
192 T tau_c = (physDiffusion_c/ conversionDiffusivity * descriptors::invCs2<T,DESCRIPTOR>())+0.5;
193 T tmp_omega_c = 1.0/tau_c;
194
195 T tau_phi = (physDiffusion_phi/ conversionDiffusivity * descriptors::invCs2<T,DESCRIPTOR>())+0.5;
196 T tmp_omega_phi = 1.0/tau_phi;
197
198
199 blockLattice.get(iX,iY).template setField<descriptors::SOURCE>(tmp_source_c*2.5e-13);
200 blockLattice.get(iX,iY).template setField<descriptors::OMEGA>(tmp_omega_c);
201 blockLattice.get(iX,iY).template setField<collision::TRT::MAGIC>(0.1*physDiffusion_c);
202
203 tpartners[0]->get(iX,iY).template setField<descriptors::SOURCE>(tmp_source_phi*2.5e-13);
204 tpartners[0]->get(iX,iY).template setField<descriptors::OMEGA>(tmp_omega_phi);
205 tpartners[0]->get(iX,iY).template setField<collision::TRT::MAGIC>(0.1*physDiffusion_phi);
206
207 }
208
209
210 else{
211 T boundaryInflow = factor1/tmp_Diffusion_E*(kappa_E*(tpartners[0]->get(iX+1,iY).computeRho()-conc[1])
212 +kappa_S/conc[0]*(blockLattice.get(iX+1,iY).computeRho()-conc[0]));
213 blockLattice.get(iX,iY).template setField<descriptors::BOUNDARY> (boundaryInflow);
214 }
215
216
217 }
218 }
219 }
220 }
int get()
Get current device.
cpu::simd::Pack< T > pow(cpu::simd::Pack< T > base, cpu::simd::Pack< T > exp)
bool intersect(int x0, int x1, int y0, int y1, int x0_, int x1_, int y0_, int y1_, int &newX0, int &newX1, int &newY0, int &newY1)