1 static char help[] = "Nonlinear, time-dependent. Developed from radiative_surface_balance.c \n"; 2 /* 3 Contributed by Steve Froehlich, Illinois Institute of Technology 4 5 Usage: 6 mpiexec -n <np> ./ex5 [options] 7 ./ex5 -help [view petsc options] 8 ./ex5 -ts_type sundials -ts_view 9 ./ex5 -da_grid_x 20 -da_grid_y 20 -log_view 10 ./ex5 -da_grid_x 20 -da_grid_y 20 -ts_type rosw -ts_atol 1.e-6 -ts_rtol 1.e-6 11 ./ex5 -drawcontours -draw_pause 0.1 -draw_fields 0,1,2,3,4 12 */ 13 14 /* 15 ----------------------------------------------------------------------- 16 17 Governing equations: 18 19 R = s*(Ea*Ta^4 - Es*Ts^4) 20 SH = p*Cp*Ch*wind*(Ta - Ts) 21 LH = p*L*Ch*wind*B(q(Ta) - q(Ts)) 22 G = k*(Tgnd - Ts)/dz 23 24 Fnet = R + SH + LH + G 25 26 du/dt = -u*(du/dx) - v*(du/dy) - 2*omeg*sin(lat)*v - (1/p)*(dP/dx) 27 dv/dt = -u*(dv/dx) - v*(dv/dy) + 2*omeg*sin(lat)*u - (1/p)*(dP/dy) 28 dTs/dt = Fnet/(Cp*dz) - Div([u*Ts, v*Ts]) + D*Lap(Ts) 29 = Fnet/(Cs*dz) - u*(dTs/dx) - v*(dTs/dy) + D*(Ts_xx + Ts_yy) 30 dp/dt = -Div([u*p,v*p]) 31 = - u*dp/dx - v*dp/dy 32 dTa/dt = Fnet/Cp 33 34 Equation of State: 35 36 P = p*R*Ts 37 38 ----------------------------------------------------------------------- 39 40 Program considers the evolution of a two dimensional atmosphere from 41 sunset to sunrise. There are two components: 42 1. Surface energy balance model to compute diabatic dT (Fnet) 43 2. Dynamical model using simplified primitive equations 44 45 Program is to be initiated at sunset and run to sunrise. 46 47 Inputs are: 48 Surface temperature 49 Dew point temperature 50 Air temperature 51 Temperature at cloud base (if clouds are present) 52 Fraction of sky covered by clouds 53 Wind speed 54 Precipitable water in centimeters 55 Wind direction 56 57 Inputs are are read in from the text file ex5_control.txt. To change an 58 input value use ex5_control.txt. 59 60 Solvers: 61 Backward Euler = default solver 62 Sundials = fastest and most accurate, requires Sundials libraries 63 64 This model is under development and should be used only as an example 65 and not as a predictive weather model. 66 */ 67 68 #include <petscts.h> 69 #include <petscdm.h> 70 #include <petscdmda.h> 71 72 /* stefan-boltzmann constant */ 73 #define SIG 0.000000056703 74 /* absorption-emission constant for surface */ 75 #define EMMSFC 1 76 /* amount of time (seconds) that passes before new flux is calculated */ 77 #define TIMESTEP 1 78 79 /* variables of interest to be solved at each grid point */ 80 typedef struct { 81 PetscScalar Ts, Ta; /* surface and air temperature */ 82 PetscScalar u, v; /* wind speed */ 83 PetscScalar p; /* density */ 84 } Field; 85 86 /* User defined variables. Used in solving for variables of interest */ 87 typedef struct { 88 DM da; /* grid */ 89 PetscScalar csoil; /* heat constant for layer */ 90 PetscScalar dzlay; /* thickness of top soil layer */ 91 PetscScalar emma; /* emission parameter */ 92 PetscScalar wind; /* wind speed */ 93 PetscScalar dewtemp; /* dew point temperature (moisture in air) */ 94 PetscScalar pressure1; /* sea level pressure */ 95 PetscScalar airtemp; /* temperature of air near boundary layer inversion */ 96 PetscScalar Ts; /* temperature at the surface */ 97 PetscScalar fract; /* fraction of sky covered by clouds */ 98 PetscScalar Tc; /* temperature at base of lowest cloud layer */ 99 PetscScalar lat; /* Latitude in degrees */ 100 PetscScalar init; /* initialization scenario */ 101 PetscScalar deep_grnd_temp; /* temperature of ground under top soil surface layer */ 102 } AppCtx; 103 104 /* Struct for visualization */ 105 typedef struct { 106 PetscBool drawcontours; /* flag - 1 indicates drawing contours */ 107 PetscViewer drawviewer; 108 PetscInt interval; 109 } MonitorCtx; 110 111 /* Inputs read in from text file */ 112 struct in { 113 PetscScalar Ts; /* surface temperature */ 114 PetscScalar Td; /* dewpoint temperature */ 115 PetscScalar Tc; /* temperature of cloud base */ 116 PetscScalar fr; /* fraction of sky covered by clouds */ 117 PetscScalar wnd; /* wind speed */ 118 PetscScalar Ta; /* air temperature */ 119 PetscScalar pwt; /* precipitable water */ 120 PetscScalar wndDir; /* wind direction */ 121 PetscScalar lat; /* latitude */ 122 PetscReal time; /* time in hours */ 123 PetscScalar init; 124 }; 125 126 /* functions */ 127 extern PetscScalar emission(PetscScalar); /* sets emission/absorption constant depending on water vapor content */ 128 extern PetscScalar calc_q(PetscScalar); /* calculates specific humidity */ 129 extern PetscScalar mph2mpers(PetscScalar); /* converts miles per hour to meters per second */ 130 extern PetscScalar Lconst(PetscScalar); /* calculates latent heat constant taken from Satellite estimates of wind speed and latent heat flux over the global oceans., Bentamy et al. */ 131 extern PetscScalar fahr_to_cel(PetscScalar); /* converts Fahrenheit to Celsius */ 132 extern PetscScalar cel_to_fahr(PetscScalar); /* converts Celsius to Fahrenheit */ 133 extern PetscScalar calcmixingr(PetscScalar, PetscScalar); /* calculates mixing ratio */ 134 extern PetscScalar cloud(PetscScalar); /* cloud radiative parameterization */ 135 extern PetscErrorCode FormInitialSolution(DM, Vec, void *); /* Specifies initial conditions for the system of equations (PETSc defined function) */ 136 extern PetscErrorCode RhsFunc(TS, PetscReal, Vec, Vec, void *); /* Specifies the user defined functions (PETSc defined function) */ 137 extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *); /* Specifies output and visualization tools (PETSc defined function) */ 138 extern PetscErrorCode readinput(struct in *put); /* reads input from text file */ 139 extern PetscErrorCode calcfluxs(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar *); /* calculates upward IR from surface */ 140 extern PetscErrorCode calcfluxa(PetscScalar, PetscScalar, PetscScalar, PetscScalar *); /* calculates downward IR from atmosphere */ 141 extern PetscErrorCode sensibleflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar *); /* calculates sensible heat flux */ 142 extern PetscErrorCode potential_temperature(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar *); /* calculates potential temperature */ 143 extern PetscErrorCode latentflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar *); /* calculates latent heat flux */ 144 extern PetscErrorCode calc_gflux(PetscScalar, PetscScalar, PetscScalar *); /* calculates flux between top soil layer and underlying earth */ 145 146 int main(int argc, char **argv) { 147 PetscInt time; /* amount of loops */ 148 struct in put; 149 PetscScalar rh; /* relative humidity */ 150 PetscScalar x; /* memory varialbe for relative humidity calculation */ 151 PetscScalar deep_grnd_temp; /* temperature of ground under top soil surface layer */ 152 PetscScalar emma; /* absorption-emission constant for air */ 153 PetscScalar pressure1 = 101300; /* surface pressure */ 154 PetscScalar mixratio; /* mixing ratio */ 155 PetscScalar airtemp; /* temperature of air near boundary layer inversion */ 156 PetscScalar dewtemp; /* dew point temperature */ 157 PetscScalar sfctemp; /* temperature at surface */ 158 PetscScalar pwat; /* total column precipitable water */ 159 PetscScalar cloudTemp; /* temperature at base of cloud */ 160 AppCtx user; /* user-defined work context */ 161 MonitorCtx usermonitor; /* user-defined monitor context */ 162 TS ts; 163 SNES snes; 164 DM da; 165 Vec T, rhs; /* solution vector */ 166 Mat J; /* Jacobian matrix */ 167 PetscReal ftime, dt; 168 PetscInt steps, dof = 5; 169 PetscBool use_coloring = PETSC_TRUE; 170 MatFDColoring matfdcoloring = 0; 171 PetscBool monitor_off = PETSC_FALSE; 172 173 PetscFunctionBeginUser; 174 PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 175 176 /* Inputs */ 177 PetscCall(readinput(&put)); 178 179 sfctemp = put.Ts; 180 dewtemp = put.Td; 181 cloudTemp = put.Tc; 182 airtemp = put.Ta; 183 pwat = put.pwt; 184 185 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial Temperature = %g\n", (double)sfctemp)); /* input surface temperature */ 186 187 deep_grnd_temp = sfctemp - 10; /* set underlying ground layer temperature */ 188 emma = emission(pwat); /* accounts for radiative effects of water vapor */ 189 190 /* Converts from Fahrenheit to Celsuis */ 191 sfctemp = fahr_to_cel(sfctemp); 192 airtemp = fahr_to_cel(airtemp); 193 dewtemp = fahr_to_cel(dewtemp); 194 cloudTemp = fahr_to_cel(cloudTemp); 195 deep_grnd_temp = fahr_to_cel(deep_grnd_temp); 196 197 /* Converts from Celsius to Kelvin */ 198 sfctemp += 273; 199 airtemp += 273; 200 dewtemp += 273; 201 cloudTemp += 273; 202 deep_grnd_temp += 273; 203 204 /* Calculates initial relative humidity */ 205 x = calcmixingr(dewtemp, pressure1); 206 mixratio = calcmixingr(sfctemp, pressure1); 207 rh = (x / mixratio) * 100; 208 209 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial RH = %.1f percent\n\n", (double)rh)); /* prints initial relative humidity */ 210 211 time = 3600 * put.time; /* sets amount of timesteps to run model */ 212 213 /* Configure PETSc TS solver */ 214 /*------------------------------------------*/ 215 216 /* Create grid */ 217 PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, 20, 20, PETSC_DECIDE, PETSC_DECIDE, dof, 1, NULL, NULL, &da)); 218 PetscCall(DMSetFromOptions(da)); 219 PetscCall(DMSetUp(da)); 220 PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 221 222 /* Define output window for each variable of interest */ 223 PetscCall(DMDASetFieldName(da, 0, "Ts")); 224 PetscCall(DMDASetFieldName(da, 1, "Ta")); 225 PetscCall(DMDASetFieldName(da, 2, "u")); 226 PetscCall(DMDASetFieldName(da, 3, "v")); 227 PetscCall(DMDASetFieldName(da, 4, "p")); 228 229 /* set values for appctx */ 230 user.da = da; 231 user.Ts = sfctemp; 232 user.fract = put.fr; /* fraction of sky covered by clouds */ 233 user.dewtemp = dewtemp; /* dew point temperature (mositure in air) */ 234 user.csoil = 2000000; /* heat constant for layer */ 235 user.dzlay = 0.08; /* thickness of top soil layer */ 236 user.emma = emma; /* emission parameter */ 237 user.wind = put.wnd; /* wind spped */ 238 user.pressure1 = pressure1; /* sea level pressure */ 239 user.airtemp = airtemp; /* temperature of air near boundar layer inversion */ 240 user.Tc = cloudTemp; /* temperature at base of lowest cloud layer */ 241 user.init = put.init; /* user chosen initiation scenario */ 242 user.lat = 70 * 0.0174532; /* converts latitude degrees to latitude in radians */ 243 user.deep_grnd_temp = deep_grnd_temp; /* temp in lowest ground layer */ 244 245 /* set values for MonitorCtx */ 246 usermonitor.drawcontours = PETSC_FALSE; 247 PetscCall(PetscOptionsHasName(NULL, NULL, "-drawcontours", &usermonitor.drawcontours)); 248 if (usermonitor.drawcontours) { 249 PetscReal bounds[] = {1000.0, -1000., -1000., -1000., 1000., -1000., 1000., -1000., 1000, -1000, 100700, 100800}; 250 PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 300, 300, &usermonitor.drawviewer)); 251 PetscCall(PetscViewerDrawSetBounds(usermonitor.drawviewer, dof, bounds)); 252 } 253 usermonitor.interval = 1; 254 PetscCall(PetscOptionsGetInt(NULL, NULL, "-monitor_interval", &usermonitor.interval, NULL)); 255 256 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 257 Extract global vectors from DA; 258 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 259 PetscCall(DMCreateGlobalVector(da, &T)); 260 PetscCall(VecDuplicate(T, &rhs)); /* r: vector to put the computed right hand side */ 261 262 PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 263 PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 264 PetscCall(TSSetType(ts, TSBEULER)); 265 PetscCall(TSSetRHSFunction(ts, rhs, RhsFunc, &user)); 266 267 /* Set Jacobian evaluation routine - use coloring to compute finite difference Jacobian efficiently */ 268 PetscCall(DMSetMatType(da, MATAIJ)); 269 PetscCall(DMCreateMatrix(da, &J)); 270 PetscCall(TSGetSNES(ts, &snes)); 271 if (use_coloring) { 272 ISColoring iscoloring; 273 PetscCall(DMCreateColoring(da, IS_COLORING_GLOBAL, &iscoloring)); 274 PetscCall(MatFDColoringCreate(J, iscoloring, &matfdcoloring)); 275 PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); 276 PetscCall(MatFDColoringSetUp(J, iscoloring, matfdcoloring)); 277 PetscCall(ISColoringDestroy(&iscoloring)); 278 PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, ts)); 279 PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring)); 280 } else { 281 PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault, NULL)); 282 } 283 284 /* Define what to print for ts_monitor option */ 285 PetscCall(PetscOptionsHasName(NULL, NULL, "-monitor_off", &monitor_off)); 286 if (!monitor_off) PetscCall(TSMonitorSet(ts, Monitor, &usermonitor, NULL)); 287 PetscCall(FormInitialSolution(da, T, &user)); 288 dt = TIMESTEP; /* initial time step */ 289 ftime = TIMESTEP * time; 290 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "time %" PetscInt_FMT ", ftime %g hour, TIMESTEP %g\n", time, (double)(ftime / 3600), (double)dt)); 291 292 PetscCall(TSSetTimeStep(ts, dt)); 293 PetscCall(TSSetMaxSteps(ts, time)); 294 PetscCall(TSSetMaxTime(ts, ftime)); 295 PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 296 PetscCall(TSSetSolution(ts, T)); 297 PetscCall(TSSetDM(ts, da)); 298 299 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 300 Set runtime options 301 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 302 PetscCall(TSSetFromOptions(ts)); 303 304 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 305 Solve nonlinear system 306 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 307 PetscCall(TSSolve(ts, T)); 308 PetscCall(TSGetSolveTime(ts, &ftime)); 309 PetscCall(TSGetStepNumber(ts, &steps)); 310 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solution T after %g hours %" PetscInt_FMT " steps\n", (double)(ftime / 3600), steps)); 311 312 if (matfdcoloring) PetscCall(MatFDColoringDestroy(&matfdcoloring)); 313 if (usermonitor.drawcontours) PetscCall(PetscViewerDestroy(&usermonitor.drawviewer)); 314 PetscCall(MatDestroy(&J)); 315 PetscCall(VecDestroy(&T)); 316 PetscCall(VecDestroy(&rhs)); 317 PetscCall(TSDestroy(&ts)); 318 PetscCall(DMDestroy(&da)); 319 320 PetscCall(PetscFinalize()); 321 return 0; 322 } 323 /*****************************end main program********************************/ 324 /*****************************************************************************/ 325 /*****************************************************************************/ 326 /*****************************************************************************/ 327 PetscErrorCode calcfluxs(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar fract, PetscScalar cloudTemp, PetscScalar *flux) { 328 PetscFunctionBeginUser; 329 *flux = SIG * ((EMMSFC * emma * PetscPowScalarInt(airtemp, 4)) + (EMMSFC * fract * (1 - emma) * PetscPowScalarInt(cloudTemp, 4)) - (EMMSFC * PetscPowScalarInt(sfctemp, 4))); /* calculates flux using Stefan-Boltzmann relation */ 330 PetscFunctionReturn(0); 331 } 332 333 PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar *flux) /* this function is not currently called upon */ 334 { 335 PetscScalar emm = 0.001; 336 337 PetscFunctionBeginUser; 338 *flux = SIG * (-emm * (PetscPowScalarInt(airtemp, 4))); /* calculates flux usinge Stefan-Boltzmann relation */ 339 PetscFunctionReturn(0); 340 } 341 PetscErrorCode sensibleflux(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar wind, PetscScalar *sheat) { 342 PetscScalar density = 1; /* air density */ 343 PetscScalar Cp = 1005; /* heat capicity for dry air */ 344 PetscScalar wndmix; /* temperature change from wind mixing: wind*Ch */ 345 346 PetscFunctionBeginUser; 347 wndmix = 0.0025 + 0.0042 * wind; /* regression equation valid for neutral and stable BL */ 348 *sheat = density * Cp * wndmix * (airtemp - sfctemp); /* calculates sensible heat flux */ 349 PetscFunctionReturn(0); 350 } 351 352 PetscErrorCode latentflux(PetscScalar sfctemp, PetscScalar dewtemp, PetscScalar wind, PetscScalar pressure1, PetscScalar *latentheat) { 353 PetscScalar density = 1; /* density of dry air */ 354 PetscScalar q; /* actual specific humitity */ 355 PetscScalar qs; /* saturation specific humidity */ 356 PetscScalar wndmix; /* temperature change from wind mixing: wind*Ch */ 357 PetscScalar beta = .4; /* moisture availability */ 358 PetscScalar mr; /* mixing ratio */ 359 PetscScalar lhcnst; /* latent heat of vaporization constant = 2501000 J/kg at 0c */ 360 /* latent heat of saturation const = 2834000 J/kg */ 361 /* latent heat of fusion const = 333700 J/kg */ 362 363 PetscFunctionBeginUser; 364 wind = mph2mpers(wind); /* converts wind from mph to meters per second */ 365 wndmix = 0.0025 + 0.0042 * wind; /* regression equation valid for neutral BL */ 366 lhcnst = Lconst(sfctemp); /* calculates latent heat of evaporation */ 367 mr = calcmixingr(sfctemp, pressure1); /* calculates saturation mixing ratio */ 368 qs = calc_q(mr); /* calculates saturation specific humidty */ 369 mr = calcmixingr(dewtemp, pressure1); /* calculates mixing ratio */ 370 q = calc_q(mr); /* calculates specific humidty */ 371 372 *latentheat = density * wndmix * beta * lhcnst * (q - qs); /* calculates latent heat flux */ 373 PetscFunctionReturn(0); 374 } 375 376 PetscErrorCode potential_temperature(PetscScalar temp, PetscScalar pressure1, PetscScalar pressure2, PetscScalar sfctemp, PetscScalar *pottemp) { 377 PetscScalar kdry; /* poisson constant for dry atmosphere */ 378 PetscScalar pavg; /* average atmospheric pressure */ 379 /* PetscScalar mixratio; mixing ratio */ 380 /* PetscScalar kmoist; poisson constant for moist atmosphere */ 381 382 PetscFunctionBeginUser; 383 /* mixratio = calcmixingr(sfctemp,pressure1); */ 384 385 /* initialize poisson constant */ 386 kdry = 0.2854; 387 /* kmoist = 0.2854*(1 - 0.24*mixratio); */ 388 389 pavg = ((0.7 * pressure1) + pressure2) / 2; /* calculates simple average press */ 390 *pottemp = temp * (PetscPowScalar((pressure1 / pavg), kdry)); /* calculates potential temperature */ 391 PetscFunctionReturn(0); 392 } 393 extern PetscScalar calcmixingr(PetscScalar dtemp, PetscScalar pressure1) { 394 PetscScalar e; /* vapor pressure */ 395 PetscScalar mixratio; /* mixing ratio */ 396 397 dtemp = dtemp - 273; /* converts from Kelvin to Celsuis */ 398 e = 6.11 * (PetscPowScalar(10, ((7.5 * dtemp) / (237.7 + dtemp)))); /* converts from dew point temp to vapor pressure */ 399 e = e * 100; /* converts from hPa to Pa */ 400 mixratio = (0.622 * e) / (pressure1 - e); /* computes mixing ratio */ 401 mixratio = mixratio * 1; /* convert to g/Kg */ 402 403 return mixratio; 404 } 405 extern PetscScalar calc_q(PetscScalar rv) { 406 PetscScalar specific_humidity; /* define specific humidity variable */ 407 specific_humidity = rv / (1 + rv); /* calculates specific humidity */ 408 return specific_humidity; 409 } 410 411 PetscErrorCode calc_gflux(PetscScalar sfctemp, PetscScalar deep_grnd_temp, PetscScalar *Gflux) { 412 PetscScalar k; /* thermal conductivity parameter */ 413 PetscScalar n = 0.38; /* value of soil porosity */ 414 PetscScalar dz = 1; /* depth of layer between soil surface and deep soil layer */ 415 PetscScalar unit_soil_weight = 2700; /* unit soil weight in kg/m^3 */ 416 417 PetscFunctionBeginUser; 418 k = ((0.135 * (1 - n) * unit_soil_weight) + 64.7) / (unit_soil_weight - (0.947 * (1 - n) * unit_soil_weight)); /* dry soil conductivity */ 419 *Gflux = (k * (deep_grnd_temp - sfctemp) / dz); /* calculates flux from deep ground layer */ 420 PetscFunctionReturn(0); 421 } 422 extern PetscScalar emission(PetscScalar pwat) { 423 PetscScalar emma; 424 425 emma = 0.725 + 0.17 * PetscLog10Real(PetscRealPart(pwat)); 426 427 return emma; 428 } 429 extern PetscScalar cloud(PetscScalar fract) { 430 PetscScalar emma = 0; 431 432 /* modifies radiative balance depending on cloud cover */ 433 if (fract >= 0.9) emma = 1; 434 else if (0.9 > fract && fract >= 0.8) emma = 0.9; 435 else if (0.8 > fract && fract >= 0.7) emma = 0.85; 436 else if (0.7 > fract && fract >= 0.6) emma = 0.75; 437 else if (0.6 > fract && fract >= 0.5) emma = 0.65; 438 else if (0.4 > fract && fract >= 0.3) emma = emma * 1.086956; 439 return emma; 440 } 441 extern PetscScalar Lconst(PetscScalar sfctemp) { 442 PetscScalar Lheat; 443 sfctemp -= 273; /* converts from kelvin to celsius */ 444 Lheat = 4186.8 * (597.31 - 0.5625 * sfctemp); /* calculates latent heat constant */ 445 return Lheat; 446 } 447 extern PetscScalar mph2mpers(PetscScalar wind) { 448 wind = ((wind * 1.6 * 1000) / 3600); /* converts wind from mph to meters per second */ 449 return wind; 450 } 451 extern PetscScalar fahr_to_cel(PetscScalar temp) { 452 temp = (5 * (temp - 32)) / 9; /* converts from farhrenheit to celsuis */ 453 return temp; 454 } 455 extern PetscScalar cel_to_fahr(PetscScalar temp) { 456 temp = ((temp * 9) / 5) + 32; /* converts from celsuis to farhrenheit */ 457 return temp; 458 } 459 PetscErrorCode readinput(struct in *put) { 460 int i; 461 char x; 462 FILE *ifp; 463 double tmp; 464 465 PetscFunctionBeginUser; 466 ifp = fopen("ex5_control.txt", "r"); 467 PetscCheck(ifp, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Unable to open input file"); 468 for (i = 0; i < 110; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 469 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 470 put->Ts = tmp; 471 472 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 473 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 474 put->Td = tmp; 475 476 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 477 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 478 put->Ta = tmp; 479 480 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 481 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 482 put->Tc = tmp; 483 484 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 485 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 486 put->fr = tmp; 487 488 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 489 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 490 put->wnd = tmp; 491 492 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 493 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 494 put->pwt = tmp; 495 496 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 497 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 498 put->wndDir = tmp; 499 500 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 501 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 502 put->time = tmp; 503 504 for (i = 0; i < 63; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 505 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 506 put->init = tmp; 507 PetscFunctionReturn(0); 508 } 509 510 /* ------------------------------------------------------------------- */ 511 PetscErrorCode FormInitialSolution(DM da, Vec Xglobal, void *ctx) { 512 AppCtx *user = (AppCtx *)ctx; /* user-defined application context */ 513 PetscInt i, j, xs, ys, xm, ym, Mx, My; 514 Field **X; 515 516 PetscFunctionBeginUser; 517 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 518 519 /* Get pointers to vector data */ 520 PetscCall(DMDAVecGetArray(da, Xglobal, &X)); 521 522 /* Get local grid boundaries */ 523 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 524 525 /* Compute function over the locally owned part of the grid */ 526 527 if (user->init == 1) { 528 for (j = ys; j < ys + ym; j++) { 529 for (i = xs; i < xs + xm; i++) { 530 X[j][i].Ts = user->Ts - i * 0.0001; 531 X[j][i].Ta = X[j][i].Ts - 5; 532 X[j][i].u = 0; 533 X[j][i].v = 0; 534 X[j][i].p = 1.25; 535 if ((j == 5 || j == 6) && (i == 4 || i == 5)) X[j][i].p += 0.00001; 536 if ((j == 5 || j == 6) && (i == 12 || i == 13)) X[j][i].p += 0.00001; 537 } 538 } 539 } else { 540 for (j = ys; j < ys + ym; j++) { 541 for (i = xs; i < xs + xm; i++) { 542 X[j][i].Ts = user->Ts; 543 X[j][i].Ta = X[j][i].Ts - 5; 544 X[j][i].u = 0; 545 X[j][i].v = 0; 546 X[j][i].p = 1.25; 547 } 548 } 549 } 550 551 /* Restore vectors */ 552 PetscCall(DMDAVecRestoreArray(da, Xglobal, &X)); 553 PetscFunctionReturn(0); 554 } 555 556 /* 557 RhsFunc - Evaluates nonlinear function F(u). 558 559 Input Parameters: 560 . ts - the TS context 561 . t - current time 562 . Xglobal - input vector 563 . F - output vector 564 . ptr - optional user-defined context, as set by SNESSetFunction() 565 566 Output Parameter: 567 . F - rhs function vector 568 */ 569 PetscErrorCode RhsFunc(TS ts, PetscReal t, Vec Xglobal, Vec F, void *ctx) { 570 AppCtx *user = (AppCtx *)ctx; /* user-defined application context */ 571 DM da = user->da; 572 PetscInt i, j, Mx, My, xs, ys, xm, ym; 573 PetscReal dhx, dhy; 574 Vec localT; 575 Field **X, **Frhs; /* structures that contain variables of interest and left hand side of governing equations respectively */ 576 PetscScalar csoil = user->csoil; /* heat constant for layer */ 577 PetscScalar dzlay = user->dzlay; /* thickness of top soil layer */ 578 PetscScalar emma = user->emma; /* emission parameter */ 579 PetscScalar wind = user->wind; /* wind speed */ 580 PetscScalar dewtemp = user->dewtemp; /* dew point temperature (moisture in air) */ 581 PetscScalar pressure1 = user->pressure1; /* sea level pressure */ 582 PetscScalar airtemp = user->airtemp; /* temperature of air near boundary layer inversion */ 583 PetscScalar fract = user->fract; /* fraction of the sky covered by clouds */ 584 PetscScalar Tc = user->Tc; /* temperature at base of lowest cloud layer */ 585 PetscScalar lat = user->lat; /* latitude */ 586 PetscScalar Cp = 1005.7; /* specific heat of air at constant pressure */ 587 PetscScalar Rd = 287.058; /* gas constant for dry air */ 588 PetscScalar diffconst = 1000; /* diffusion coefficient */ 589 PetscScalar f = 2 * 0.0000727 * PetscSinScalar(lat); /* coriolis force */ 590 PetscScalar deep_grnd_temp = user->deep_grnd_temp; /* temp in lowest ground layer */ 591 PetscScalar Ts, u, v, p; 592 PetscScalar u_abs, u_plus, u_minus, v_abs, v_plus, v_minus; 593 594 PetscScalar sfctemp1, fsfc1, Ra; 595 PetscScalar sheat; /* sensible heat flux */ 596 PetscScalar latentheat; /* latent heat flux */ 597 PetscScalar groundflux; /* flux from conduction of deep ground layer in contact with top soil */ 598 PetscInt xend, yend; 599 600 PetscFunctionBeginUser; 601 PetscCall(DMGetLocalVector(da, &localT)); 602 PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 603 604 dhx = (PetscReal)(Mx - 1) / (5000 * (Mx - 1)); /* dhx = 1/dx; assume 2D space domain: [0.0, 1.e5] x [0.0, 1.e5] */ 605 dhy = (PetscReal)(My - 1) / (5000 * (Mx - 1)); /* dhy = 1/dy; */ 606 607 /* 608 Scatter ghost points to local vector,using the 2-step process 609 DAGlobalToLocalBegin(),DAGlobalToLocalEnd(). 610 By placing code between these two statements, computations can be 611 done while messages are in transition. 612 */ 613 PetscCall(DMGlobalToLocalBegin(da, Xglobal, INSERT_VALUES, localT)); 614 PetscCall(DMGlobalToLocalEnd(da, Xglobal, INSERT_VALUES, localT)); 615 616 /* Get pointers to vector data */ 617 PetscCall(DMDAVecGetArrayRead(da, localT, &X)); 618 PetscCall(DMDAVecGetArray(da, F, &Frhs)); 619 620 /* Get local grid boundaries */ 621 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 622 623 /* Compute function over the locally owned part of the grid */ 624 /* the interior points */ 625 xend = xs + xm; 626 yend = ys + ym; 627 for (j = ys; j < yend; j++) { 628 for (i = xs; i < xend; i++) { 629 Ts = X[j][i].Ts; 630 u = X[j][i].u; 631 v = X[j][i].v; 632 p = X[j][i].p; /*P = X[j][i].P; */ 633 634 sfctemp1 = (double)Ts; 635 PetscCall(calcfluxs(sfctemp1, airtemp, emma, fract, Tc, &fsfc1)); /* calculates surface net radiative flux */ 636 PetscCall(sensibleflux(sfctemp1, airtemp, wind, &sheat)); /* calculate sensible heat flux */ 637 PetscCall(latentflux(sfctemp1, dewtemp, wind, pressure1, &latentheat)); /* calculates latent heat flux */ 638 PetscCall(calc_gflux(sfctemp1, deep_grnd_temp, &groundflux)); /* calculates flux from earth below surface soil layer by conduction */ 639 PetscCall(calcfluxa(sfctemp1, airtemp, emma, &Ra)); /* Calculates the change in downward radiative flux */ 640 fsfc1 = fsfc1 + latentheat + sheat + groundflux; /* adds radiative, sensible heat, latent heat, and ground heat flux yielding net flux */ 641 642 /* convective coefficients for upwinding */ 643 u_abs = PetscAbsScalar(u); 644 u_plus = .5 * (u + u_abs); /* u if u>0; 0 if u<0 */ 645 u_minus = .5 * (u - u_abs); /* u if u <0; 0 if u>0 */ 646 647 v_abs = PetscAbsScalar(v); 648 v_plus = .5 * (v + v_abs); /* v if v>0; 0 if v<0 */ 649 v_minus = .5 * (v - v_abs); /* v if v <0; 0 if v>0 */ 650 651 /* Solve governing equations */ 652 /* P = p*Rd*Ts; */ 653 654 /* du/dt -> time change of east-west component of the wind */ 655 Frhs[j][i].u = -u_plus * (u - X[j][i - 1].u) * dhx - u_minus * (X[j][i + 1].u - u) * dhx /* - u(du/dx) */ 656 - v_plus * (u - X[j - 1][i].u) * dhy - v_minus * (X[j + 1][i].u - u) * dhy /* - v(du/dy) */ 657 - (Rd / p) * (Ts * (X[j][i + 1].p - X[j][i - 1].p) * 0.5 * dhx + p * 0 * (X[j][i + 1].Ts - X[j][i - 1].Ts) * 0.5 * dhx) /* -(R/p)[Ts(dp/dx)+ p(dTs/dx)] */ 658 /* -(1/p)*(X[j][i+1].P - X[j][i-1].P)*dhx */ 659 + f * v; 660 661 /* dv/dt -> time change of north-south component of the wind */ 662 Frhs[j][i].v = -u_plus * (v - X[j][i - 1].v) * dhx - u_minus * (X[j][i + 1].v - v) * dhx /* - u(dv/dx) */ 663 - v_plus * (v - X[j - 1][i].v) * dhy - v_minus * (X[j + 1][i].v - v) * dhy /* - v(dv/dy) */ 664 - (Rd / p) * (Ts * (X[j + 1][i].p - X[j - 1][i].p) * 0.5 * dhy + p * 0 * (X[j + 1][i].Ts - X[j - 1][i].Ts) * 0.5 * dhy) /* -(R/p)[Ts(dp/dy)+ p(dTs/dy)] */ 665 /* -(1/p)*(X[j+1][i].P - X[j-1][i].P)*dhy */ 666 - f * u; 667 668 /* dT/dt -> time change of temperature */ 669 Frhs[j][i].Ts = (fsfc1 / (csoil * dzlay)) /* Fnet/(Cp*dz) diabatic change in T */ 670 - u_plus * (Ts - X[j][i - 1].Ts) * dhx - u_minus * (X[j][i + 1].Ts - Ts) * dhx /* - u*(dTs/dx) advection x */ 671 - v_plus * (Ts - X[j - 1][i].Ts) * dhy - v_minus * (X[j + 1][i].Ts - Ts) * dhy /* - v*(dTs/dy) advection y */ 672 + diffconst * ((X[j][i + 1].Ts - 2 * Ts + X[j][i - 1].Ts) * dhx * dhx /* + D(Ts_xx + Ts_yy) diffusion */ 673 + (X[j + 1][i].Ts - 2 * Ts + X[j - 1][i].Ts) * dhy * dhy); 674 675 /* dp/dt -> time change of */ 676 Frhs[j][i].p = -u_plus * (p - X[j][i - 1].p) * dhx - u_minus * (X[j][i + 1].p - p) * dhx /* - u*(dp/dx) */ 677 - v_plus * (p - X[j - 1][i].p) * dhy - v_minus * (X[j + 1][i].p - p) * dhy; /* - v*(dp/dy) */ 678 679 Frhs[j][i].Ta = Ra / Cp; /* dTa/dt time change of air temperature */ 680 } 681 } 682 683 /* Restore vectors */ 684 PetscCall(DMDAVecRestoreArrayRead(da, localT, &X)); 685 PetscCall(DMDAVecRestoreArray(da, F, &Frhs)); 686 PetscCall(DMRestoreLocalVector(da, &localT)); 687 PetscFunctionReturn(0); 688 } 689 690 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec T, void *ctx) { 691 const PetscScalar *array; 692 MonitorCtx *user = (MonitorCtx *)ctx; 693 PetscViewer viewer = user->drawviewer; 694 PetscReal norm; 695 696 PetscFunctionBeginUser; 697 PetscCall(VecNorm(T, NORM_INFINITY, &norm)); 698 699 if (step % user->interval == 0) { 700 PetscCall(VecGetArrayRead(T, &array)); 701 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step %" PetscInt_FMT ", time %8.1f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f\n", step, (double)time, (double)(((array[0] - 273) * 9) / 5 + 32), (double)(((array[1] - 273) * 9) / 5 + 32), (double)array[2], (double)array[3], (double)array[4], (double)array[5])); 702 PetscCall(VecRestoreArrayRead(T, &array)); 703 } 704 705 if (user->drawcontours) PetscCall(VecView(T, viewer)); 706 PetscFunctionReturn(0); 707 } 708 709 /*TEST 710 711 build: 712 requires: !complex !single 713 714 test: 715 args: -ts_max_steps 130 -monitor_interval 60 716 output_file: output/ex5.out 717 requires: !complex !single 718 localrunfiles: ex5_control.txt 719 720 test: 721 suffix: 2 722 nsize: 4 723 args: -ts_max_steps 130 -monitor_interval 60 724 output_file: output/ex5.out 725 localrunfiles: ex5_control.txt 726 requires: !complex !single 727 728 TEST*/ 729