Skip to content

Commit 51282ec

Browse files
glesurglesur
committed
fix the computation of velocity gradients in the viscosity module when coarsening is enabled.
1 parent 392d89d commit 51282ec

1 file changed

Lines changed: 126 additions & 48 deletions

File tree

‎src/fluid/viscosity.cpp‎

Lines changed: 126 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -111,9 +111,9 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
111111
IdefixArray1D<real> x1l = this->data->xl[IDIR];
112112
IdefixArray1D<real> x2l = this->data->xl[JDIR];
113113
IdefixArray1D<real> x3l = this->data->xl[KDIR];
114-
IdefixArray1D<real> dx1 = this->data->dx[IDIR];
115-
IdefixArray1D<real> dx2 = this->data->dx[JDIR];
116-
IdefixArray1D<real> dx3 = this->data->dx[KDIR];
114+
IdefixArray1D<real> dx1Array = this->data->dx[IDIR];
115+
IdefixArray1D<real> dx2Array = this->data->dx[JDIR];
116+
IdefixArray1D<real> dx3Array = this->data->dx[KDIR];
117117

118118
// Fargo variables use to correct energy fluxes
119119
#if HAVE_ENERGY
@@ -129,6 +129,14 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
129129
IdefixArray1D<real> tanx2m = this->data->tanx2m;
130130
#endif
131131

132+
// Coarsening parameters
133+
bool haveGridCoarsening = data->haveGridCoarsening != GridCoarsening::disabled;
134+
bool haveGridCoarseningX1 = data->coarseningDirection[IDIR];
135+
bool haveGridCoarseningX2 = data->coarseningDirection[JDIR];
136+
bool haveGridCoarseningX3 = data->coarseningDirection[KDIR];
137+
auto coarseningLevelX1 = data->coarseningLevel[IDIR];
138+
auto coarseningLevelX2 = data->coarseningLevel[JDIR];
139+
auto coarseningLevelX3 = data->coarseningLevel[KDIR];
132140

133141
HydroModuleStatus haveViscosity = this->status.status;
134142

@@ -201,18 +209,42 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
201209
etaC2 = eta2 = eta2Constant;
202210
}
203211

204-
EXPAND( dVxi = D_DX_I(Vc,VX1)/dx1(i); ,
205-
dVyi = D_DX_I(Vc,VX2)/dx1(i); ,
206-
dVzi = D_DX_I(Vc,VX3)/dx1(i); )
212+
// Compute spacing
213+
real dx1 = dx1Array(i);
214+
real dx2 = dx2Array(j);
215+
real dx3 = dx3Array(k);
216+
217+
// Correct spacing for grid coarsening
218+
219+
if(haveGridCoarsening) {
220+
if(haveGridCoarseningX1) {
221+
const int factor = 1 << (coarseningLevelX1(k,j) - 1);
222+
dx1 *= factor;
223+
}
224+
225+
if(haveGridCoarseningX2) {
226+
const int factor = 1 << (coarseningLevelX2(k,i) - 1);
227+
dx2 *= factor;
228+
}
229+
230+
if(haveGridCoarseningX3) {
231+
const int factor = 1 << (coarseningLevelX3(j,i) - 1);
232+
dx3 *= factor;
233+
}
234+
}
235+
236+
EXPAND( dVxi = D_DX_I(Vc,VX1)/dx1; ,
237+
dVyi = D_DX_I(Vc,VX2)/dx1; ,
238+
dVzi = D_DX_I(Vc,VX3)/dx1; )
207239

208240
#if DIMENSIONS >= 2
209-
EXPAND( dVxj = D_DY_I(Vc,VX1)/dx2(j); ,
210-
dVyj = D_DY_I(Vc,VX2)/dx2(j); ,
211-
dVzj = D_DY_I(Vc,VX3)/dx2(j); )
241+
EXPAND( dVxj = D_DY_I(Vc,VX1)/dx2; ,
242+
dVyj = D_DY_I(Vc,VX2)/dx2; ,
243+
dVzj = D_DY_I(Vc,VX3)/dx2; )
212244
#if DIMENSIONS == 3
213-
dVxk = D_DZ_I(Vc,VX1)/dx3(k);
214-
dVyk = D_DZ_I(Vc,VX2)/dx3(k);
215-
dVzk = D_DZ_I(Vc,VX3)/dx3(k);
245+
dVxk = D_DZ_I(Vc,VX1)/dx3;
246+
dVyk = D_DZ_I(Vc,VX2)/dx3;
247+
dVzk = D_DZ_I(Vc,VX3)/dx3;
216248
#endif
217249
#endif
218250

@@ -234,14 +266,14 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
234266
tau_xy = eta1*(dVxj + dVyi);
235267
SELECT( ,
236268
,
237-
tau_xz = eta1*0.5*(x1(i)+x1(i-1))/dx1(i)
269+
tau_xz = eta1*0.5*(x1(i)+x1(i-1))/dx1
238270
*(Vc(VX3,k,j,i)/x1(i) - Vc(VX3,k,j,i-1)/x1(i-1)); )
239271

240272
tau_xz = eta1*(dVxk + dVzi);
241273

242274
// compute tau_zz at cell center
243-
divV = D_EXPAND( 0.5*(Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i-1))/dx1(i) + Vc(VX1,k,j,i)/x1(i),
244-
+0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2(j) ,
275+
divV = D_EXPAND( 0.5*(Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i-1))/dx1 + Vc(VX1,k,j,i)/x1(i),
276+
+0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2 ,
245277
);
246278
tau_zz = 2.0*etaC1*Vc(VX1,k,j,i)/x1(i) + (etaC2 - (2.0/3.0)*etaC1)*divV;
247279

@@ -269,14 +301,14 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
269301

270302

271303
// compute tau_yy at cell center
272-
divV = D_EXPAND( 0.5*(Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i-1))/dx1(i) + Vc(VX1,k,j,i)/x1(i),
273-
+0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2(j)/x1(i) ,
274-
+0.5*(Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/dx3(k) );
304+
divV = D_EXPAND( 0.5*(Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i-1))/dx1 + Vc(VX1,k,j,i)/x1(i),
305+
+0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2/x1(i) ,
306+
+0.5*(Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/dx3 );
275307

276308
#if DIMENSIONS == 1
277309
tau_yy = 2.0*etaC1*( Vc(VX1,k,j,i)/x1(i)) + (etaC2 - (2.0/3.0)*etaC1)*divV;
278310
#else
279-
tau_yy = 2.0*etaC1*( 0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2(j)/x1(i)
311+
tau_yy = 2.0*etaC1*( 0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2/x1(i)
280312
+ Vc(VX1,k,j,i)/x1(i))
281313
+ (etaC2 - (2.0/3.0)*etaC1)*divV;
282314
#endif
@@ -325,14 +357,14 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
325357

326358
// Compute tau_yy and tau_zz at cell center
327359

328-
divV = D_EXPAND( 0.5*(Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i-1))/dx1(i) + 2.0*Vc(VX1,k,j,i)/x1(i),
329-
+0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2(j)/x1(i)
360+
divV = D_EXPAND( 0.5*(Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i-1))/dx1 + 2.0*Vc(VX1,k,j,i)/x1(i),
361+
+0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2/x1(i)
330362
+Vc(VX2,k,j,i)/x1(i)*tan_1 ,
331-
+0.5*(Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/dx3(k)/x1(i)*s_1 );
363+
+0.5*(Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/dx3/x1(i)*s_1 );
332364

333365
tau_yy = Vc(VX1,k,j,i)/x1(i);
334366
#if DIMENSIONS >= 2
335-
tau_yy += 0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2(j)/x1(i);
367+
tau_yy += 0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2/x1(i);
336368
#endif
337369
tau_yy = 2.0*etaC1*tau_yy + (etaC2-2.0/3.0*etaC1)*divV;
338370

@@ -341,7 +373,7 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
341373
tau_zz += Vc(VX2,k,j,i)*tan_1/x1(i);
342374
#endif
343375
#if DIMENSIONS == 3
344-
tau_zz += 0.5*(Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/dx3(k)/x1(i)*s_1;
376+
tau_zz += 0.5*(Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/dx3/x1(i)*s_1;
345377
#endif
346378
tau_zz = 2.0*etaC1*tau_zz + (etaC2-2.0/3.0*etaC1)*divV;
347379

@@ -426,18 +458,41 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
426458
etaC2 = eta2 = eta2Constant;
427459
}
428460

429-
EXPAND( dVxi = D_DX_J(Vc,VX1)/dx1(i); ,
430-
dVyi = D_DX_J(Vc,VX2)/dx1(i); ,
431-
dVzi = D_DX_J(Vc,VX3)/dx1(i); )
461+
// Compute spacing
462+
real dx1 = dx1Array(i);
463+
real dx2 = dx2Array(j);
464+
real dx3 = dx3Array(k);
465+
466+
// Correct spacing for grid coarsening
467+
if(haveGridCoarsening) {
468+
if(haveGridCoarseningX1) {
469+
const int factor = 1 << (coarseningLevelX1(k,j) - 1);
470+
dx1 *= factor;
471+
}
432472

433-
EXPAND( dVxj = D_DY_J(Vc,VX1)/dx2(j); ,
434-
dVyj = D_DY_J(Vc,VX2)/dx2(j); ,
435-
dVzj = D_DY_J(Vc,VX3)/dx2(j); )
473+
if(haveGridCoarseningX2) {
474+
const int factor = 1 << (coarseningLevelX2(k,i) - 1);
475+
dx2 *= factor;
476+
}
477+
478+
if(haveGridCoarseningX3) {
479+
const int factor = 1 << (coarseningLevelX3(j,i) - 1);
480+
dx3 *= factor;
481+
}
482+
}
483+
484+
EXPAND( dVxi = D_DX_J(Vc,VX1)/dx1; ,
485+
dVyi = D_DX_J(Vc,VX2)/dx1; ,
486+
dVzi = D_DX_J(Vc,VX3)/dx1; )
487+
488+
EXPAND( dVxj = D_DY_J(Vc,VX1)/dx2; ,
489+
dVyj = D_DY_J(Vc,VX2)/dx2; ,
490+
dVzj = D_DY_J(Vc,VX3)/dx2; )
436491

437492
#if DIMENSIONS == 3
438-
dVxk = D_DZ_J(Vc,VX1)/dx3(k);
439-
dVyk = D_DZ_J(Vc,VX2)/dx3(k);
440-
dVzk = D_DZ_J(Vc,VX3)/dx3(k);
493+
dVxk = D_DZ_J(Vc,VX1)/dx3;
494+
dVyk = D_DZ_J(Vc,VX2)/dx3;
495+
dVzk = D_DZ_J(Vc,VX3)/dx3;
441496
#endif
442497

443498
#if GEOMETRY == CARTESIAN
@@ -508,8 +563,8 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
508563
+ dVzk/x1(i)*s_1 );
509564

510565
// tau_xy is initially cell centered since it is involved in the source term
511-
tau_xy = etaC1*( 0.5*(Vc(VX1,k,j+1,i)-Vc(VX1,k,j-1,i))/x1(i)/dx2(j)
512-
+0.5*(Vc(VX2,k,j,i+1)-Vc(VX2,k,j,i-1))/dx1(i)
566+
tau_xy = etaC1*( 0.5*(Vc(VX1,k,j+1,i)-Vc(VX1,k,j-1,i))/x1(i)/dx2
567+
+0.5*(Vc(VX2,k,j,i+1)-Vc(VX2,k,j,i-1))/dx1
513568
- Vc(VX2,k,j,i)/x1(i));
514569
tau_yy = 2.0*eta1*(dVyj/x1(i) + vx1i/x1(i)) + (eta2 - (2.0/3.0)*eta1)*divV;
515570

@@ -539,17 +594,17 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
539594
s_1 = ONE_F/s_1;
540595
}
541596

542-
divV = D_EXPAND( 0.5*(Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i-1))/dx1(i) + 2.0*Vc(VX1,k,j,i)/x1(i),
543-
+0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2(j)/x1(i)
597+
divV = D_EXPAND( 0.5*(Vc(VX1,k,j,i+1) - Vc(VX1,k,j,i-1))/dx1 + 2.0*Vc(VX1,k,j,i)/x1(i),
598+
+0.5*(Vc(VX2,k,j+1,i) - Vc(VX2,k,j-1,i))/dx2/x1(i)
544599
+Vc(VX2,k,j,i)/x1(i)*tan_1 ,
545-
+0.5*(Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/dx3(k)/x1(i)*s_1 );
600+
+0.5*(Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/dx3/x1(i)*s_1 );
546601

547602
tau_zz = Vc(VX1,k,j,i)/x1(i);
548603
#if COMPONENTS >= 2
549604
tau_zz += Vc(VX2,k,j,i)/x1(i)*tan_1;
550605
#endif
551606
#if DIMENSIONS == 3
552-
tau_zz += 0.5*(Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/dx3(k)/x1(i)*s_1;
607+
tau_zz += 0.5*(Vc(VX3,k+1,j,i) - Vc(VX3,k-1,j,i))/dx3/x1(i)*s_1;
553608
#endif
554609
tau_zz = 2*eta1*tau_zz + (eta2-2.0/3.0*eta1)*divV;
555610

@@ -639,15 +694,38 @@ void Viscosity::AddViscousFlux(int dir, const real t, const IdefixArray4D<real>
639694
etaC2 = eta2 = eta2Constant;
640695
}
641696

642-
dVxi = D_DX_K(Vc,VX1)/dx1(i);
643-
dVyi = D_DX_K(Vc,VX2)/dx1(i);
644-
dVzi = D_DX_K(Vc,VX3)/dx1(i);
645-
dVxj = D_DY_K(Vc,VX1)/dx2(j);
646-
dVyj = D_DY_K(Vc,VX2)/dx2(j);
647-
dVzj = D_DY_K(Vc,VX3)/dx2(j);
648-
dVxk = D_DZ_K(Vc,VX1)/dx3(k);
649-
dVyk = D_DZ_K(Vc,VX2)/dx3(k);
650-
dVzk = D_DZ_K(Vc,VX3)/dx3(k);
697+
// Compute spacing
698+
real dx1 = dx1Array(i);
699+
real dx2 = dx2Array(j);
700+
real dx3 = dx3Array(k);
701+
702+
// Correct spacing for grid coarsening
703+
if(haveGridCoarsening) {
704+
if(haveGridCoarseningX1) {
705+
const int factor = 1 << (coarseningLevelX1(k,j) - 1);
706+
dx1 *= factor;
707+
}
708+
709+
if(haveGridCoarseningX2) {
710+
const int factor = 1 << (coarseningLevelX2(k,i) - 1);
711+
dx2 *= factor;
712+
}
713+
714+
if(haveGridCoarseningX3) {
715+
const int factor = 1 << (coarseningLevelX3(j,i) - 1);
716+
dx3 *= factor;
717+
}
718+
}
719+
720+
dVxi = D_DX_K(Vc,VX1)/dx1;
721+
dVyi = D_DX_K(Vc,VX2)/dx1;
722+
dVzi = D_DX_K(Vc,VX3)/dx1;
723+
dVxj = D_DY_K(Vc,VX1)/dx2;
724+
dVyj = D_DY_K(Vc,VX2)/dx2;
725+
dVzj = D_DY_K(Vc,VX3)/dx2;
726+
dVxk = D_DZ_K(Vc,VX1)/dx3;
727+
dVyk = D_DZ_K(Vc,VX2)/dx3;
728+
dVzk = D_DZ_K(Vc,VX3)/dx3;
651729

652730
#if GEOMETRY == CARTESIAN
653731
divV = dVxi + dVyj + dVzk;

0 commit comments

Comments
 (0)