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