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