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 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 variable 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 prunejacobian = PETSC_FALSE; 174 175 PetscFunctionBeginUser; 176 PetscCall(PetscInitialize(&argc, &argv, NULL, 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 Celsius */ 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 speed */ 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, "-prune_jacobian", &prunejacobian, NULL)); 307 PetscCall(TSSetFromOptions(ts)); 308 if (prunejacobian && matfdcoloring) { 309 PetscRandom rctx; 310 Vec Tdot; 311 /* Compute the Jacobian with randomized vector values to capture the sparsity pattern for coloring */ 312 PetscCall(VecDuplicate(T, &Tdot)); 313 PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 314 PetscCall(PetscRandomSetInterval(rctx, 1.0, 2.0)); 315 PetscCall(VecSetRandom(T, rctx)); 316 PetscCall(VecSetRandom(Tdot, rctx)); 317 PetscCall(PetscRandomDestroy(&rctx)); 318 PetscCall(TSSetUp(ts)); 319 PetscCall(TSComputeIJacobian(ts, 0.0, T, Tdot, 0.12345, J, J, PETSC_FALSE)); 320 PetscCall(VecDestroy(&Tdot)); 321 PetscCall(MatFDColoringDestroy(&matfdcoloring)); 322 PetscCall(TSPruneIJacobianColor(ts, J, J)); 323 } 324 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 325 Solve nonlinear system 326 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 327 PetscCall(FormInitialSolution(da, T, &user)); 328 PetscCall(TSSolve(ts, T)); 329 PetscCall(TSGetSolveTime(ts, &ftime)); 330 PetscCall(TSGetStepNumber(ts, &steps)); 331 PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solution T after %g hours %" PetscInt_FMT " steps\n", (double)(ftime / 3600), steps)); 332 333 if (matfdcoloring) PetscCall(MatFDColoringDestroy(&matfdcoloring)); 334 if (usermonitor.drawcontours) PetscCall(PetscViewerDestroy(&usermonitor.drawviewer)); 335 PetscCall(MatDestroy(&J)); 336 PetscCall(VecDestroy(&T)); 337 PetscCall(VecDestroy(&rhs)); 338 PetscCall(TSDestroy(&ts)); 339 PetscCall(DMDestroy(&da)); 340 341 PetscCall(PetscFinalize()); 342 return 0; 343 } 344 /*****************************end main program********************************/ 345 /*****************************************************************************/ 346 /*****************************************************************************/ 347 /*****************************************************************************/ 348 PetscErrorCode calcfluxs(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar fract, PetscScalar cloudTemp, PetscScalar *flux) 349 { 350 PetscFunctionBeginUser; 351 *flux = SIG * ((EMMSFC * emma * PetscPowScalarInt(airtemp, 4)) + (EMMSFC * fract * (1 - emma) * PetscPowScalarInt(cloudTemp, 4)) - (EMMSFC * PetscPowScalarInt(sfctemp, 4))); /* calculates flux using Stefan-Boltzmann relation */ 352 PetscFunctionReturn(PETSC_SUCCESS); 353 } 354 355 PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar *flux) /* this function is not currently called upon */ 356 { 357 PetscScalar emm = 0.001; 358 359 PetscFunctionBeginUser; 360 *flux = SIG * (-emm * (PetscPowScalarInt(airtemp, 4))); /* calculates flux usinge Stefan-Boltzmann relation */ 361 PetscFunctionReturn(PETSC_SUCCESS); 362 } 363 PetscErrorCode sensibleflux(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar wind, PetscScalar *sheat) 364 { 365 PetscScalar density = 1; /* air density */ 366 PetscScalar Cp = 1005; /* heat capacity for dry air */ 367 PetscScalar wndmix; /* temperature change from wind mixing: wind*Ch */ 368 369 PetscFunctionBeginUser; 370 wndmix = 0.0025 + 0.0042 * wind; /* regression equation valid for neutral and stable BL */ 371 *sheat = density * Cp * wndmix * (airtemp - sfctemp); /* calculates sensible heat flux */ 372 PetscFunctionReturn(PETSC_SUCCESS); 373 } 374 375 PetscErrorCode latentflux(PetscScalar sfctemp, PetscScalar dewtemp, PetscScalar wind, PetscScalar pressure1, PetscScalar *latentheat) 376 { 377 PetscScalar density = 1; /* density of dry air */ 378 PetscScalar q; /* actual specific humitity */ 379 PetscScalar qs; /* saturation specific humidity */ 380 PetscScalar wndmix; /* temperature change from wind mixing: wind*Ch */ 381 PetscScalar beta = .4; /* moisture availability */ 382 PetscScalar mr; /* mixing ratio */ 383 PetscScalar lhcnst; /* latent heat of vaporization constant = 2501000 J/kg at 0c */ 384 /* latent heat of saturation const = 2834000 J/kg */ 385 /* latent heat of fusion const = 333700 J/kg */ 386 387 PetscFunctionBeginUser; 388 wind = mph2mpers(wind); /* converts wind from mph to meters per second */ 389 wndmix = 0.0025 + 0.0042 * wind; /* regression equation valid for neutral BL */ 390 lhcnst = Lconst(sfctemp); /* calculates latent heat of evaporation */ 391 mr = calcmixingr(sfctemp, pressure1); /* calculates saturation mixing ratio */ 392 qs = calc_q(mr); /* calculates saturation specific humidty */ 393 mr = calcmixingr(dewtemp, pressure1); /* calculates mixing ratio */ 394 q = calc_q(mr); /* calculates specific humidty */ 395 396 *latentheat = density * wndmix * beta * lhcnst * (q - qs); /* calculates latent heat flux */ 397 PetscFunctionReturn(PETSC_SUCCESS); 398 } 399 400 PetscErrorCode potential_temperature(PetscScalar temp, PetscScalar pressure1, PetscScalar pressure2, PetscScalar sfctemp, PetscScalar *pottemp) 401 { 402 PetscScalar kdry; /* poisson constant for dry atmosphere */ 403 PetscScalar pavg; /* average atmospheric pressure */ 404 /* PetscScalar mixratio; mixing ratio */ 405 /* PetscScalar kmoist; poisson constant for moist atmosphere */ 406 407 PetscFunctionBeginUser; 408 /* mixratio = calcmixingr(sfctemp,pressure1); */ 409 410 /* initialize poisson constant */ 411 kdry = 0.2854; 412 /* kmoist = 0.2854*(1 - 0.24*mixratio); */ 413 414 pavg = ((0.7 * pressure1) + pressure2) / 2; /* calculates simple average press */ 415 *pottemp = temp * (PetscPowScalar(pressure1 / pavg, kdry)); /* calculates potential temperature */ 416 PetscFunctionReturn(PETSC_SUCCESS); 417 } 418 extern PetscScalar calcmixingr(PetscScalar dtemp, PetscScalar pressure1) 419 { 420 PetscScalar e; /* vapor pressure */ 421 PetscScalar mixratio; /* mixing ratio */ 422 423 dtemp = dtemp - 273; /* converts from Kelvin to Celsius */ 424 e = 6.11 * (PetscPowScalar(10, (7.5 * dtemp) / (237.7 + dtemp))); /* converts from dew point temp to vapor pressure */ 425 e = e * 100; /* converts from hPa to Pa */ 426 mixratio = (0.622 * e) / (pressure1 - e); /* computes mixing ratio */ 427 mixratio = mixratio * 1; /* convert to g/Kg */ 428 429 return mixratio; 430 } 431 extern PetscScalar calc_q(PetscScalar rv) 432 { 433 PetscScalar specific_humidity; /* define specific humidity variable */ 434 specific_humidity = rv / (1 + rv); /* calculates specific humidity */ 435 return specific_humidity; 436 } 437 438 PetscErrorCode calc_gflux(PetscScalar sfctemp, PetscScalar deep_grnd_temp, PetscScalar *Gflux) 439 { 440 PetscScalar k; /* thermal conductivity parameter */ 441 PetscScalar n = 0.38; /* value of soil porosity */ 442 PetscScalar dz = 1; /* depth of layer between soil surface and deep soil layer */ 443 PetscScalar unit_soil_weight = 2700; /* unit soil weight in kg/m^3 */ 444 445 PetscFunctionBeginUser; 446 k = ((0.135 * (1 - n) * unit_soil_weight) + 64.7) / (unit_soil_weight - (0.947 * (1 - n) * unit_soil_weight)); /* dry soil conductivity */ 447 *Gflux = (k * (deep_grnd_temp - sfctemp) / dz); /* calculates flux from deep ground layer */ 448 PetscFunctionReturn(PETSC_SUCCESS); 449 } 450 extern PetscScalar emission(PetscScalar pwat) 451 { 452 PetscScalar emma; 453 454 emma = 0.725 + 0.17 * PetscLog10Real(PetscRealPart(pwat)); 455 456 return emma; 457 } 458 extern PetscScalar cloud(PetscScalar fract) 459 { 460 PetscScalar emma = 0; 461 462 /* modifies radiative balance depending on cloud cover */ 463 if (fract >= 0.9) emma = 1; 464 else if (0.9 > fract && fract >= 0.8) emma = 0.9; 465 else if (0.8 > fract && fract >= 0.7) emma = 0.85; 466 else if (0.7 > fract && fract >= 0.6) emma = 0.75; 467 else if (0.6 > fract && fract >= 0.5) emma = 0.65; 468 else if (0.4 > fract && fract >= 0.3) emma = emma * 1.086956; 469 return emma; 470 } 471 extern PetscScalar Lconst(PetscScalar sfctemp) 472 { 473 PetscScalar Lheat; 474 sfctemp -= 273; /* converts from kelvin to celsius */ 475 Lheat = 4186.8 * (597.31 - 0.5625 * sfctemp); /* calculates latent heat constant */ 476 return Lheat; 477 } 478 extern PetscScalar mph2mpers(PetscScalar wind) 479 { 480 wind = ((wind * 1.6 * 1000) / 3600); /* converts wind from mph to meters per second */ 481 return wind; 482 } 483 extern PetscScalar fahr_to_cel(PetscScalar temp) 484 { 485 temp = (5 * (temp - 32)) / 9; /* converts from farhrenheit to celsius */ 486 return temp; 487 } 488 extern PetscScalar cel_to_fahr(PetscScalar temp) 489 { 490 temp = ((temp * 9) / 5) + 32; /* converts from celsius to farhrenheit */ 491 return temp; 492 } 493 PetscErrorCode readinput(struct in *put) 494 { 495 int i; 496 char x; 497 FILE *ifp; 498 double tmp; 499 500 PetscFunctionBeginUser; 501 ifp = fopen("ex5_control.txt", "r"); 502 PetscCheck(ifp, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Unable to open input file"); 503 for (i = 0; i < 110; 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->Ts = 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->Td = 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->Ta = 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->Tc = tmp; 518 519 for (i = 0; i < 43; 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->fr = tmp; 522 523 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 524 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 525 put->wnd = tmp; 526 527 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 528 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 529 put->pwt = tmp; 530 531 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 532 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 533 put->wndDir = tmp; 534 535 for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 536 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 537 put->time = tmp; 538 539 for (i = 0; i < 63; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 540 PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 541 put->init = tmp; 542 PetscFunctionReturn(PETSC_SUCCESS); 543 } 544 545 /* ------------------------------------------------------------------- */ 546 PetscErrorCode FormInitialSolution(DM da, Vec Xglobal, PetscCtx ctx) 547 { 548 AppCtx *user = (AppCtx *)ctx; /* user-defined application context */ 549 PetscInt i, j, xs, ys, xm, ym, Mx, My; 550 Field **X; 551 552 PetscFunctionBeginUser; 553 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)); 554 555 /* Get pointers to vector data */ 556 PetscCall(DMDAVecGetArray(da, Xglobal, &X)); 557 558 /* Get local grid boundaries */ 559 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 560 561 /* Compute function over the locally owned part of the grid */ 562 563 if (user->init == 1) { 564 for (j = ys; j < ys + ym; j++) { 565 for (i = xs; i < xs + xm; i++) { 566 X[j][i].Ts = user->Ts - i * 0.0001; 567 X[j][i].Ta = X[j][i].Ts - 5; 568 X[j][i].u = 0; 569 X[j][i].v = 0; 570 X[j][i].p = 1.25; 571 if ((j == 5 || j == 6) && (i == 4 || i == 5)) X[j][i].p += 0.00001; 572 if ((j == 5 || j == 6) && (i == 12 || i == 13)) X[j][i].p += 0.00001; 573 } 574 } 575 } else { 576 for (j = ys; j < ys + ym; j++) { 577 for (i = xs; i < xs + xm; i++) { 578 X[j][i].Ts = user->Ts; 579 X[j][i].Ta = X[j][i].Ts - 5; 580 X[j][i].u = 0; 581 X[j][i].v = 0; 582 X[j][i].p = 1.25; 583 } 584 } 585 } 586 587 /* Restore vectors */ 588 PetscCall(DMDAVecRestoreArray(da, Xglobal, &X)); 589 PetscFunctionReturn(PETSC_SUCCESS); 590 } 591 592 /* 593 RhsFunc - Evaluates nonlinear function F(u). 594 595 Input Parameters: 596 . ts - the TS context 597 . t - current time 598 . Xglobal - input vector 599 . F - output vector 600 . ptr - optional user-defined context, as set by SNESSetFunction() 601 602 Output Parameter: 603 . F - rhs function vector 604 */ 605 PetscErrorCode RhsFunc(TS ts, PetscReal t, Vec Xglobal, Vec F, PetscCtx ctx) 606 { 607 AppCtx *user = (AppCtx *)ctx; /* user-defined application context */ 608 DM da = user->da; 609 PetscInt i, j, Mx, My, xs, ys, xm, ym; 610 PetscReal dhx, dhy; 611 Vec localT; 612 Field **X, **Frhs; /* structures that contain variables of interest and left hand side of governing equations respectively */ 613 PetscScalar csoil = user->csoil; /* heat constant for layer */ 614 PetscScalar dzlay = user->dzlay; /* thickness of top soil layer */ 615 PetscScalar emma = user->emma; /* emission parameter */ 616 PetscScalar wind = user->wind; /* wind speed */ 617 PetscScalar dewtemp = user->dewtemp; /* dew point temperature (moisture in air) */ 618 PetscScalar pressure1 = user->pressure1; /* sea level pressure */ 619 PetscScalar airtemp = user->airtemp; /* temperature of air near boundary layer inversion */ 620 PetscScalar fract = user->fract; /* fraction of the sky covered by clouds */ 621 PetscScalar Tc = user->Tc; /* temperature at base of lowest cloud layer */ 622 PetscScalar lat = user->lat; /* latitude */ 623 PetscScalar Cp = 1005.7; /* specific heat of air at constant pressure */ 624 PetscScalar Rd = 287.058; /* gas constant for dry air */ 625 PetscScalar diffconst = 1000; /* diffusion coefficient */ 626 PetscScalar f = 2 * 0.0000727 * PetscSinScalar(lat); /* coriolis force */ 627 PetscScalar deep_grnd_temp = user->deep_grnd_temp; /* temp in lowest ground layer */ 628 PetscScalar Ts, u, v, p; 629 PetscScalar u_abs, u_plus, u_minus, v_abs, v_plus, v_minus; 630 631 PetscScalar sfctemp1, fsfc1, Ra; 632 PetscScalar sheat; /* sensible heat flux */ 633 PetscScalar latentheat; /* latent heat flux */ 634 PetscScalar groundflux; /* flux from conduction of deep ground layer in contact with top soil */ 635 PetscInt xend, yend; 636 637 PetscFunctionBeginUser; 638 PetscCall(DMGetLocalVector(da, &localT)); 639 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)); 640 641 dhx = (PetscReal)(Mx - 1) / (5000 * (Mx - 1)); /* dhx = 1/dx; assume 2D space domain: [0.0, 1.e5] x [0.0, 1.e5] */ 642 dhy = (PetscReal)(My - 1) / (5000 * (Mx - 1)); /* dhy = 1/dy; */ 643 644 /* 645 Scatter ghost points to local vector,using the 2-step process 646 DAGlobalToLocalBegin(),DAGlobalToLocalEnd(). 647 By placing code between these two statements, computations can be 648 done while messages are in transition. 649 */ 650 PetscCall(DMGlobalToLocalBegin(da, Xglobal, INSERT_VALUES, localT)); 651 PetscCall(DMGlobalToLocalEnd(da, Xglobal, INSERT_VALUES, localT)); 652 653 /* Get pointers to vector data */ 654 PetscCall(DMDAVecGetArrayRead(da, localT, &X)); 655 PetscCall(DMDAVecGetArray(da, F, &Frhs)); 656 657 /* Get local grid boundaries */ 658 PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 659 660 /* Compute function over the locally owned part of the grid */ 661 /* the interior points */ 662 xend = xs + xm; 663 yend = ys + ym; 664 for (j = ys; j < yend; j++) { 665 for (i = xs; i < xend; i++) { 666 Ts = X[j][i].Ts; 667 u = X[j][i].u; 668 v = X[j][i].v; 669 p = X[j][i].p; /*P = X[j][i].P; */ 670 671 sfctemp1 = (double)Ts; 672 PetscCall(calcfluxs(sfctemp1, airtemp, emma, fract, Tc, &fsfc1)); /* calculates surface net radiative flux */ 673 PetscCall(sensibleflux(sfctemp1, airtemp, wind, &sheat)); /* calculate sensible heat flux */ 674 PetscCall(latentflux(sfctemp1, dewtemp, wind, pressure1, &latentheat)); /* calculates latent heat flux */ 675 PetscCall(calc_gflux(sfctemp1, deep_grnd_temp, &groundflux)); /* calculates flux from earth below surface soil layer by conduction */ 676 PetscCall(calcfluxa(sfctemp1, airtemp, emma, &Ra)); /* Calculates the change in downward radiative flux */ 677 fsfc1 = fsfc1 + latentheat + sheat + groundflux; /* adds radiative, sensible heat, latent heat, and ground heat flux yielding net flux */ 678 679 /* convective coefficients for upwinding */ 680 u_abs = PetscAbsScalar(u); 681 u_plus = .5 * (u + u_abs); /* u if u>0; 0 if u<0 */ 682 u_minus = .5 * (u - u_abs); /* u if u <0; 0 if u>0 */ 683 684 v_abs = PetscAbsScalar(v); 685 v_plus = .5 * (v + v_abs); /* v if v>0; 0 if v<0 */ 686 v_minus = .5 * (v - v_abs); /* v if v <0; 0 if v>0 */ 687 688 /* Solve governing equations */ 689 /* P = p*Rd*Ts; */ 690 691 /* du/dt -> time change of east-west component of the wind */ 692 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) */ 693 - v_plus * (u - X[j - 1][i].u) * dhy - v_minus * (X[j + 1][i].u - u) * dhy /* - v(du/dy) */ 694 - (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)] */ 695 /* -(1/p)*(X[j][i+1].P - X[j][i-1].P)*dhx */ 696 + f * v; 697 698 /* dv/dt -> time change of north-south component of the wind */ 699 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) */ 700 - v_plus * (v - X[j - 1][i].v) * dhy - v_minus * (X[j + 1][i].v - v) * dhy /* - v(dv/dy) */ 701 - (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)] */ 702 /* -(1/p)*(X[j+1][i].P - X[j-1][i].P)*dhy */ 703 - f * u; 704 705 /* dT/dt -> time change of temperature */ 706 Frhs[j][i].Ts = (fsfc1 / (csoil * dzlay)) /* Fnet/(Cp*dz) diabatic change in T */ 707 - u_plus * (Ts - X[j][i - 1].Ts) * dhx - u_minus * (X[j][i + 1].Ts - Ts) * dhx /* - u*(dTs/dx) advection x */ 708 - v_plus * (Ts - X[j - 1][i].Ts) * dhy - v_minus * (X[j + 1][i].Ts - Ts) * dhy /* - v*(dTs/dy) advection y */ 709 + diffconst * ((X[j][i + 1].Ts - 2 * Ts + X[j][i - 1].Ts) * dhx * dhx /* + D(Ts_xx + Ts_yy) diffusion */ 710 + (X[j + 1][i].Ts - 2 * Ts + X[j - 1][i].Ts) * dhy * dhy); 711 712 /* dp/dt -> time change of */ 713 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) */ 714 - v_plus * (p - X[j - 1][i].p) * dhy - v_minus * (X[j + 1][i].p - p) * dhy; /* - v*(dp/dy) */ 715 716 Frhs[j][i].Ta = Ra / Cp; /* dTa/dt time change of air temperature */ 717 } 718 } 719 720 /* Restore vectors */ 721 PetscCall(DMDAVecRestoreArrayRead(da, localT, &X)); 722 PetscCall(DMDAVecRestoreArray(da, F, &Frhs)); 723 PetscCall(DMRestoreLocalVector(da, &localT)); 724 PetscFunctionReturn(PETSC_SUCCESS); 725 } 726 727 PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec T, PetscCtx ctx) 728 { 729 const PetscScalar *array; 730 MonitorCtx *user = (MonitorCtx *)ctx; 731 PetscViewer viewer = user->drawviewer; 732 PetscReal norm; 733 734 PetscFunctionBeginUser; 735 PetscCall(VecNorm(T, NORM_INFINITY, &norm)); 736 737 if (step % user->interval == 0) { 738 PetscCall(VecGetArrayRead(T, &array)); 739 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])); 740 PetscCall(VecRestoreArrayRead(T, &array)); 741 } 742 743 if (user->drawcontours) PetscCall(VecView(T, viewer)); 744 PetscFunctionReturn(PETSC_SUCCESS); 745 } 746 747 /*TEST 748 749 build: 750 requires: !complex !single 751 752 test: 753 args: -ts_max_steps 130 -monitor_interval 60 754 output_file: output/ex5.out 755 requires: !complex !single 756 localrunfiles: ex5_control.txt 757 758 test: 759 suffix: 2 760 nsize: 4 761 args: -ts_max_steps 130 -monitor_interval 60 762 output_file: output/ex5.out 763 localrunfiles: ex5_control.txt 764 requires: !complex !single 765 766 # Test TSPruneIJacobianColor() for improved FD coloring 767 test: 768 suffix: 3 769 nsize: 4 770 args: -ts_max_steps 130 -monitor_interval 60 -prune_jacobian -mat_coloring_view 771 requires: !complex !single 772 localrunfiles: ex5_control.txt 773 774 TEST*/ 775