CG  Version 25
setTemperatureBC.h
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // Macro: Apply BC's on the Temperature
3 //
4 // There are 3 cases:
5 // (1) apply a dirichlet BC (OPTION=dirichlet)
6 // (2) extrapolate ghost pts on dirichlet BC's (OPTION=extrapolateGhost)
7 // (3) apply a mixed BC (OPTION=mixed)
8 //
9 // Macro args:
10 //
11 // tc : component to assign
12 // NAME : name of of the calling function (for comments)
13 // BCNAME : noSlipWall, inflowWithVelocityGiven etc.
14 // OPTION: dirichlet, mixed, extrapolateGhost
15 // ----------------------------------------------------------------------------
16 #beginMacro setTemperatureBC(tc,NAME,BCNAME,OPTION)
17  if( assignTemperature )
18  {
19  #If #OPTION == "dirichlet" || #OPTION == "mixed" || #OPTION == "extrapolateGhost"
20  #Else
21  Overture::abort("ERROR in calling setTemperatureBC macro with option=OPTION");
22  #End
23 
24  FILE *& debugFile = parameters.dbase.get<FILE* >("debugFile");
25  FILE *& pDebugFile = parameters.dbase.get<FILE* >("pDebugFile");
26  ForBoundary(side,axis)
27  {
28 
29  if( mg.boundaryCondition(side,axis)==BCNAME )
30  {
31 
32  if( interfaceType(side,axis,grid)!=Parameters::noInterface )
33  { // This is an interface between domains
34 
35  // for now we only know about interfaces at no-slip walls:
36  assert( mg.boundaryCondition(side,axis)==noSlipWall );
37 
38  // what about BC's applied at t=0 before the boundary data is set ??
39  // if( parameters.dbase.get<int >("globalStepNumber") < 2 ) continue; // ********************* TEMP *****
40 
41  // if this is an iterface we should turn off the TZ forcing for the boundary condition since we want
42  // to use the boundary data instead.
43  #ifdef USE_PPP
44  realSerialArray uLocal; getLocalArrayWithGhostBoundaries(u,uLocal);
45  #else
46  const realSerialArray & uLocal = u;
47  #endif
48 
49  Index Ib1,Ib2,Ib3;
50  getBoundaryIndex(mg.gridIndexRange(),side,axis,Ib1,Ib2,Ib3);
51  bool ok = ParallelUtility::getLocalArrayBounds(u,uLocal,Ib1,Ib2,Ib3);
52 
53  if( debug() & 4 )
54  {
55  printP("NAME:applyBC: apply a mixed BC on the interface (side,axis,grid)=(%i,%i,%i) %f*T + %f*T.n \n",
56  side,axis,grid,mixedCoeff(tc,side,axis,grid),mixedNormalCoeff(tc,side,axis,grid));
57 
58  fprintf(pDebugFile,"NAME:applyBC: apply a mixed BC on the interface (side,axis,grid)=(%i,%i,%i) %f*T + %f*T.n \n",
59  side,axis,grid,mixedCoeff(tc,side,axis,grid),mixedNormalCoeff(tc,side,axis,grid));
60  if( pBoundaryData[side][axis]==NULL )
61  {
62  if( !ok )
63  fprintf(pDebugFile," NAME:applyBC on T: pBoundaryData[side][axis] == NULL! \n");
64  else
65  fprintf(pDebugFile," NAME:applyBC on T: ERROR: pBoundaryData[side][axis] == NULL! \n");
66  }
67  else
68  {
69  // RealArray & bd = *pBoundaryData[side][axis];
70  // ::display(bd,"boundaryData",pDebugFile,"%8.2e ");
71  }
72 
73  }
74 
75  assert( mixedCoeff(tc,side,axis,grid)!=0. || mixedNormalCoeff(tc,side,axis,grid)!=0. );
76 
77  if( ok && pBoundaryData[side][axis]==NULL )
78  {
79  // At t=0 when the initial conditions are being set-up (and initial conditions projected) there may
80  // not be a boundaryData array created yet for the interface. We thus create one and fill in some default values
81  // based on the current solution
82  RealArray & bd = parameters.getBoundaryData(side,axis,grid,mg);
83  bd=0.;
84 
85  #ifdef USE_PPP
86  const RealArray & normal = mg.vertexBoundaryNormalArray(side,axis);
87  #else
88  const RealArray & normal = mg.vertexBoundaryNormal(side,axis);
89  #endif
90  real a0=mixedCoeff(tc,side,axis,grid);
91  real a1=mixedNormalCoeff(tc,side,axis,grid);
92 
93  // bd(Ib1,Ib2,Ib3,tc)=a0*uLocal(Ib1,Ib2,Ib3,tc); // Assume that T.n=0 at the interface at t=0 -- could do better --
94 
95  MappedGridOperators & op = *(u.getOperators());
96  Range N(tc,tc);
97  RealArray ux(Ib1,Ib2,Ib3,N), uy(Ib1,Ib2,Ib3,N);
98  op.derivative(MappedGridOperators::xDerivative,uLocal,ux,Ib1,Ib2,Ib3,N);
99  op.derivative(MappedGridOperators::yDerivative,uLocal,uy,Ib1,Ib2,Ib3,N);
100  if( mg.numberOfDimensions()==2 )
101  {
102  bd(Ib1,Ib2,Ib3,tc)=a1*(normal(Ib1,Ib2,Ib3,0)*ux(Ib1,Ib2,Ib3,tc)+
103  normal(Ib1,Ib2,Ib3,1)*uy(Ib1,Ib2,Ib3,tc)) + a0*uLocal(Ib1,Ib2,Ib3,tc);
104  }
105  else
106  {
107  RealArray uz(Ib1,Ib2,Ib3,N);
108  op.derivative(MappedGridOperators::zDerivative,uLocal,uz,Ib1,Ib2,Ib3,N);
109  bd(Ib1,Ib2,Ib3,tc)=a1*(normal(Ib1,Ib2,Ib3,0)*ux(Ib1,Ib2,Ib3,tc)+
110  normal(Ib1,Ib2,Ib3,1)*uy(Ib1,Ib2,Ib3,tc)+
111  normal(Ib1,Ib2,Ib3,2)*uz(Ib1,Ib2,Ib3,tc)) + a0*uLocal(Ib1,Ib2,Ib3,tc);
112  }
113 
114  const bool twilightZoneFlow = parameters.dbase.get<bool >("twilightZoneFlow");
115  if( false && twilightZoneFlow ) // *******************************************************
116  {
117  fprintf(pDebugFile," NAME:applyBC on T: ********** Set TRUE value for bd at t=%8.2e ************\n",t);
118  const bool rectangular= mg.isRectangular() && !twilightZoneFlow;
119 
120  realArray & x= mg.center();
121 #ifdef USE_PPP
122  realSerialArray xLocal;
123  if( !rectangular || twilightZoneFlow )
124  getLocalArrayWithGhostBoundaries(x,xLocal);
125 #else
126  const realSerialArray & xLocal = x;
127 #endif
128 
129  OGFunction & e = *(parameters.dbase.get<OGFunction* >("exactSolution"));
130  realSerialArray ue(Ib1,Ib2,Ib3), uex(Ib1,Ib2,Ib3), uey(Ib1,Ib2,Ib3);
131  e.gd( ue ,xLocal,mg.numberOfDimensions(),rectangular,0,0,0,0,Ib1,Ib2,Ib3,tc,t);
132  e.gd( uex,xLocal,mg.numberOfDimensions(),rectangular,0,1,0,0,Ib1,Ib2,Ib3,tc,t);
133  e.gd( uey,xLocal,mg.numberOfDimensions(),rectangular,0,0,1,0,Ib1,Ib2,Ib3,tc,t);
134 
135  real a0=mixedCoeff(tc,side,axis,grid), a1=mixedNormalCoeff(tc,side,axis,grid);
136  if( mg.numberOfDimensions()==2 )
137  {
138  bd(Ib1,Ib2,Ib3,tc)=a1*(normal(Ib1,Ib2,Ib3,0)*uex(Ib1,Ib2,Ib3)+
139  normal(Ib1,Ib2,Ib3,1)*uey(Ib1,Ib2,Ib3)) + a0*ue(Ib1,Ib2,Ib3);
140  }
141  else
142  {
143  realSerialArray uez(Ib1,Ib2,Ib3);
144  e.gd( uez,xLocal,mg.numberOfDimensions(),rectangular,0,0,0,1,Ib1,Ib2,Ib3,tc,t);
145  bd(Ib1,Ib2,Ib3,tc)=a1*(normal(Ib1,Ib2,Ib3,0)*uex(Ib1,Ib2,Ib3)+
146  normal(Ib1,Ib2,Ib3,1)*uey(Ib1,Ib2,Ib3)+
147  normal(Ib1,Ib2,Ib3,2)*uez(Ib1,Ib2,Ib3)) + a0*ue(Ib1,Ib2,Ib3);
148  }
149  } // *******************************************
150  }
151 
152 
153 
154  if( pBoundaryData[side][axis]!=NULL )
155  u.getOperators()->setTwilightZoneFlow( false );
156  else
157  {
158  if( t>0. || debug() & 4 )
159  {
160  if( ok )
161  printP("$$$$ NAME:applyBC:INFO:interface but no boundaryData, t=%9.3e\n",t);
162  }
163  }
164  }
165 
166 
167  if( debug() & 4 )
168  printF("++++NAME: BCNAME: (grid,side,axis)=(%i,%i,%i) : "
169  " BC for T: %3.2f*T+%3.2f*T.n=%3.2f, \n"
170  " BC: u: %3.2f*u+%3.2f*u.n=%3.2f, v: %3.2f*v+%3.2f*v.n=%3.2f \n",
171  grid,side,axis,
172  mixedCoeff(tc,side,axis,grid), mixedNormalCoeff(tc,side,axis,grid), mixedRHS(tc,side,axis,grid),
173  mixedCoeff(uc,side,axis,grid), mixedNormalCoeff(uc,side,axis,grid), mixedRHS(uc,side,axis,grid),
174  mixedCoeff(vc,side,axis,grid), mixedNormalCoeff(vc,side,axis,grid), mixedRHS(vc,side,axis,grid)
175  );
176 
177 // if( mixedNormalCoeff(tc,side,axis,grid)!=0. ) // coeff of T.n is non-zero
178 // {
179 // mixedBoundaryConditionOnTemperature=true;
180 
181 // if( debug() & 4 )
182 // printF("++++insBC: BCNAME:adiabaticWall: (grid,side,axis)=(%i,%i,%i) : "
183 // "Mixed BC for T: %3.2f*T+%3.2f*T.n=%3.2f, \n"
184 // " BC: u: %3.2f*u+%3.2f*u.n=%3.2f=%3.2f, v: %3.2f*v+%3.2f*v.n=%3.2f \n",
185 // grid,side,axis,
186 // mixedCoeff(tc,side,axis,grid), mixedNormalCoeff(tc,side,axis,grid), mixedRHS(tc,side,axis,grid),
187 // mixedCoeff(uc,side,axis,grid), mixedNormalCoeff(uc,side,axis,grid), mixedRHS(uc,side,axis,grid),
188 // mixedCoeff(vc,side,axis,grid), mixedNormalCoeff(vc,side,axis,grid), mixedRHS(vc,side,axis,grid)
189 // );
190 // }
191 
192  if( mixedNormalCoeff(tc,side,axis,grid)==0. ) // coeff of T.n
193  {
194  // Dirichlet
195 #If #OPTION == "dirichlet"
196  u.applyBoundaryCondition(tc,dirichlet,BCTypes::boundary(side,axis),bcData,pBoundaryData,t,
197  bcParams,grid);
198 #Elif #OPTION == "extrapolateGhost"
199  u.applyBoundaryCondition(tc,extrapolate,BCTypes::boundary(side,axis),0.,t);
200 // u.applyBoundaryCondition(tc,extrapolate,BCTypes::boundary(side,axis),bcData,pBoundaryData,t,
201 // bcParams,grid);
202 #End
203  }
204  else if( bd.hasVariableCoefficientBoundaryCondition(side,axis) )
205  {
206  // -- Variable Coefficient Temperature (const coeff.) BC ---
207 
208  printF("setTemperatureBC:INFO: grid=%i (side,axis)=(%i,%i) HAS a var coeff. temperature BC!\n",
209  grid,side,axis);
210 
211  // BC is : a0(x)*T + an(x)*T.n = g
212  // a0 = varCoeff(i1,i2,i3,0)
213  // an = varCoeff(i1,i2,i3,1)
214  RealArray & varCoeff = bd.getVariableCoefficientBoundaryConditionArray(
215  BoundaryData::variableCoefficientTemperatureBC,side,axis );
216 
217  bcParams.setVariableCoefficientsArray(&varCoeff);
218 
219  u.applyBoundaryCondition(tc,mixed,BCTypes::boundary(side,axis),bcData,pBoundaryData,t,
220  bcParams,grid);
221 
222  bcParams.setVariableCoefficientsArray(NULL); // reset
223 
224  }
225  else
226  {
227  // --- Mixed or Neumann Temperature (const coeff.) BC ---
228 
229 
230 #If #OPTION == "mixed"
231 
232  // Mixed BC or Neumann
233  real a0=mixedCoeff(tc,side,axis,grid);
234  real a1=mixedNormalCoeff(tc,side,axis,grid);
235  bcParams.a.redim(3);
236  if( a0==0. && a1==1. )
237  {
238  if( debug() & 4 )
239  printF("++++NAME: BCNAME:adiabaticWall: (grid,side,axis)=(%i,%i,%i) : apply neumannBC\n",
240  grid,side,axis);
241 
242 // real b0=bcData(tc+2,side,axis,grid);
243 // u.applyBoundaryCondition(tc,neumann,BCTypes::boundary(side,axis),b0,t); // b0 ignored??
244 
245  bcParams.a(0)=a0;
246  bcParams.a(1)=a1;
247  bcParams.a(2)=mixedRHS(tc,side,axis,grid); // this is not used -- this does not work
248  if( false )
249  { // **** TEMP FIX ****
250  u.applyBoundaryCondition(tc,extrapolate,BCTypes::boundary(side,axis),0.,t);
251  }
252  else
253  {
254  u.applyBoundaryCondition(tc,neumann,BCTypes::boundary(side,axis),bcData,pBoundaryData,t,
255  bcParams,grid);
256  }
257  }
258  else
259  {
260 
261  if( debug() & 4 )
262  {
263  fPrintF(pDebugFile,"++++NAME:BCNAME:adiabaticWall: (grid,side,axis)=(%i,%i,%i) : "
264  "Mixed BC for T: %3.2f*T+%3.2f*T.n=%3.2f (t=%8.2e)\n",
265  grid,side,axis,a0,a1,bcData(tc,side,axis,grid),t);
266  }
267  if( debug() & 4 )
268  {
269  #ifndef USE_PPP
270  Index Ib1,Ib2,Ib3;
271  getBoundaryIndex(mg.gridIndexRange(),side,axis,Ib1,Ib2,Ib3);
272  RealArray & bd = parameters.getBoundaryData(side,axis,grid,mg);
273  ::display(bd(Ib1,Ib2,Ib3,tc),"NAME:BCNAME:T: RHS for mixed BC: bd(Ib1,Ib2,Ib3,tc)",pDebugFile,"%5.2f ");
274  Index Ig1,Ig2,Ig3;
275  // getGhostIndex(mg.gridIndexRange(),side,axis,Ig1,Ig2,Ig3);
276  getIndex(mg.gridIndexRange(),Ig1,Ig2,Ig3,1);
277  ::display(u(Ig1,Ig2,Ig3,tc),"NAME:BCNAME:T: BEFORE mixed BC: u(I1,I2,I3,tc)",pDebugFile,"%5.2f ");
278  #endif
279  }
280 
281  bcParams.a(0)=a0;
282  bcParams.a(1)=a1;
283  bcParams.a(2)=mixedRHS(tc,side,axis,grid);
284  if( false )
285  { // **** TEMP FIX ****
286  u.applyBoundaryCondition(tc,extrapolate,BCTypes::boundary(side,axis),0.,t);
287  }
288  else
289  {
290  u.applyBoundaryCondition(tc,mixed,BCTypes::boundary(side,axis),bcData,pBoundaryData,t,
291  bcParams,grid);
292  }
293 
294  if( debug() & 4 )
295  {
296  #ifndef USE_PPP
297  Index Ig1,Ig2,Ig3;
298  // getGhostIndex(mg.gridIndexRange(),side,axis,Ig1,Ig2,Ig3);
299  getIndex(mg.gridIndexRange(),Ig1,Ig2,Ig3,1);
300  ::display(u(Ig1,Ig2,Ig3,tc),"NAME:BCNAME:T: AFTER mixed BC: u(I1,I2,I3,tc)",pDebugFile,"%5.2f ");
301  #endif
302  }
303 
304  }
305 
306 #End
307 
308  }
309 
310  if( interfaceType(side,axis,grid)!=Parameters::noInterface )
311  { // reset TZ
312  u.getOperators()->setTwilightZoneFlow( parameters.dbase.get<bool >("twilightZoneFlow") );
313  }
314 
315  } // end if bc = BCNAME
316 
317  } // end for boundary
318 
319  // ************ try this ********* 080909
320  // u.updateGhostBoundaries(); // this is done in finish boundary conditions now
321 
322  } // end if assignTemperature
323 #endMacro
Generated on Fri Jan 4 2013 10:19:27 for CG by   doxygen 1.8.3