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