85 int x0_,
int x1_,
int y0_,
int y1_)
override {
87 int newX0, newX1, newY0, newY1;
91 newX0, newX1, newY0, newY1 ) ) {
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>();
107 std::vector<T> conc_plusX;
108 std::vector<T> conc_minusX;
109 std::vector<T> conc_plusY;
110 std::vector<T> conc_minusY;
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());
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;
120 if ((
int) materialNumber == 6 || (
int) materialNumber == 7 || (
int) materialNumber == 1 || (
int) materialNumber == 3 || (
int) materialNumber == 4 || (
int) materialNumber == 5){
123 conc_minusX.emplace_back(0.0);
124 conc_minusX.emplace_back(0.0);
127 conc_minusX.emplace_back(blockLattice.
get(iX-1,iY).computeRho());
128 conc_minusX.emplace_back(tpartners[0]->get(iX-1,iY).computeRho());
132 conc_plusX.emplace_back(0.0);
133 conc_plusX.emplace_back(0.0);
136 conc_plusX.emplace_back(blockLattice.
get(iX+1,iY).computeRho());
137 conc_plusX.emplace_back(tpartners[0]->get(iX+1,iY).computeRho());
141 conc_minusY.emplace_back(blockLattice.
get(iX, newY1).computeRho());
142 conc_minusY.emplace_back(tpartners[0]->get(iX, newY1).computeRho());
145 conc_minusY.emplace_back(blockLattice.
get(iX,iY-1).computeRho());
146 conc_minusY.emplace_back(tpartners[0]->get(iX,iY-1).computeRho());
150 conc_plusY.emplace_back(blockLattice.
get(iX,newY0).computeRho());
151 conc_plusY.emplace_back(tpartners[0]->get(iX,newY0).computeRho());
154 conc_plusY.emplace_back(blockLattice.
get(iX,iY+1).computeRho());
155 conc_plusY.emplace_back(tpartners[0]->get(iX,iY+1).computeRho());
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;
163 if ((
int) materialNumber == 6 || (
int)materialNumber == 5){
164 physDiffusion_c = diffusion_S;
165 physDiffusion_phi = kappa_S;
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]));
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]));
192 T tau_c = (physDiffusion_c/ conversionDiffusivity * descriptors::invCs2<T,DESCRIPTOR>())+0.5;
193 T tmp_omega_c = 1.0/tau_c;
195 T tau_phi = (physDiffusion_phi/ conversionDiffusivity * descriptors::invCs2<T,DESCRIPTOR>())+0.5;
196 T tmp_omega_phi = 1.0/tau_phi;
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);
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);
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);