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