Actual source code: ex5.c

  1: static char help[] = "Nonlinear, time-dependent. Developed from radiative_surface_balance.c \n";
  2: /*
  3:   Contributed by Steve Froehlich, Illinois Institute of Technology

  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: */

 14: /*
 15:    -----------------------------------------------------------------------

 17:    Governing equations:

 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

 24:         Fnet   = R + SH + LH + G

 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

 34:    Equation of State:

 36:         P = p*R*Ts

 38:    -----------------------------------------------------------------------

 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

 45:    Program is to be initiated at sunset and run to sunrise.

 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

 57:    Inputs are are read in from the text file ex5_control.txt. To change an
 58:    input value use ex5_control.txt.

 60:    Solvers:
 61:             Backward Euler = default solver
 62:             Sundials = fastest and most accurate, requires Sundials libraries

 64:    This model is under development and should be used only as an example
 65:    and not as a predictive weather model.
 66: */

 68: #include <petscts.h>
 69: #include <petscdm.h>
 70: #include <petscdmda.h>

 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

 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;

 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;

104: /* Struct for visualization */
105: typedef struct {
106:   PetscBool   drawcontours; /* flag - 1 indicates drawing contours */
107:   PetscViewer drawviewer;
108:   PetscInt    interval;
109: } MonitorCtx;

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: };

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 */

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;

175:   PetscInitialize(&argc, &argv, (char *)0, help);

177:   /* Inputs */
178:   readinput(&put);

180:   sfctemp   = put.Ts;
181:   dewtemp   = put.Td;
182:   cloudTemp = put.Tc;
183:   airtemp   = put.Ta;
184:   pwat      = put.pwt;

186:   PetscPrintf(PETSC_COMM_WORLD, "Initial Temperature = %g\n", (double)sfctemp); /* input surface temperature */

188:   deep_grnd_temp = sfctemp - 10;   /* set underlying ground layer temperature */
189:   emma           = emission(pwat); /* accounts for radiative effects of water vapor */

191:   /* Converts from Fahrenheit to Celsius */
192:   sfctemp        = fahr_to_cel(sfctemp);
193:   airtemp        = fahr_to_cel(airtemp);
194:   dewtemp        = fahr_to_cel(dewtemp);
195:   cloudTemp      = fahr_to_cel(cloudTemp);
196:   deep_grnd_temp = fahr_to_cel(deep_grnd_temp);

198:   /* Converts from Celsius to Kelvin */
199:   sfctemp += 273;
200:   airtemp += 273;
201:   dewtemp += 273;
202:   cloudTemp += 273;
203:   deep_grnd_temp += 273;

205:   /* Calculates initial relative humidity */
206:   x        = calcmixingr(dewtemp, pressure1);
207:   mixratio = calcmixingr(sfctemp, pressure1);
208:   rh       = (x / mixratio) * 100;

210:   PetscPrintf(PETSC_COMM_WORLD, "Initial RH = %.1f percent\n\n", (double)rh); /* prints initial relative humidity */

212:   time = 3600 * put.time; /* sets amount of timesteps to run model */

214:   /* Configure PETSc TS solver */
215:   /*------------------------------------------*/

217:   /* Create grid */
218:   DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, 20, 20, PETSC_DECIDE, PETSC_DECIDE, dof, 1, NULL, NULL, &da);
219:   DMSetFromOptions(da);
220:   DMSetUp(da);
221:   DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);

223:   /* Define output window for each variable of interest */
224:   DMDASetFieldName(da, 0, "Ts");
225:   DMDASetFieldName(da, 1, "Ta");
226:   DMDASetFieldName(da, 2, "u");
227:   DMDASetFieldName(da, 3, "v");
228:   DMDASetFieldName(da, 4, "p");

230:   /* set values for appctx */
231:   user.da             = da;
232:   user.Ts             = sfctemp;
233:   user.fract          = put.fr;         /* fraction of sky covered by clouds */
234:   user.dewtemp        = dewtemp;        /* dew point temperature (mositure in air) */
235:   user.csoil          = 2000000;        /* heat constant for layer */
236:   user.dzlay          = 0.08;           /* thickness of top soil layer */
237:   user.emma           = emma;           /* emission parameter */
238:   user.wind           = put.wnd;        /* wind spped */
239:   user.pressure1      = pressure1;      /* sea level pressure */
240:   user.airtemp        = airtemp;        /* temperature of air near boundar layer inversion */
241:   user.Tc             = cloudTemp;      /* temperature at base of lowest cloud layer */
242:   user.init           = put.init;       /* user chosen initiation scenario */
243:   user.lat            = 70 * 0.0174532; /* converts latitude degrees to latitude in radians */
244:   user.deep_grnd_temp = deep_grnd_temp; /* temp in lowest ground layer */

246:   /* set values for MonitorCtx */
247:   usermonitor.drawcontours = PETSC_FALSE;
248:   PetscOptionsHasName(NULL, NULL, "-drawcontours", &usermonitor.drawcontours);
249:   if (usermonitor.drawcontours) {
250:     PetscReal bounds[] = {1000.0, -1000., -1000., -1000., 1000., -1000., 1000., -1000., 1000, -1000, 100700, 100800};
251:     PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 300, 300, &usermonitor.drawviewer);
252:     PetscViewerDrawSetBounds(usermonitor.drawviewer, dof, bounds);
253:   }
254:   usermonitor.interval = 1;
255:   PetscOptionsGetInt(NULL, NULL, "-monitor_interval", &usermonitor.interval, NULL);

257:   /*  - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
258:      Extract global vectors from DA;
259:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
260:   DMCreateGlobalVector(da, &T);
261:   VecDuplicate(T, &rhs); /* r: vector to put the computed right hand side */

263:   TSCreate(PETSC_COMM_WORLD, &ts);
264:   TSSetProblemType(ts, TS_NONLINEAR);
265:   TSSetType(ts, TSBEULER);
266:   TSSetRHSFunction(ts, rhs, RhsFunc, &user);

268:   /* Set Jacobian evaluation routine - use coloring to compute finite difference Jacobian efficiently */
269:   DMSetMatType(da, MATAIJ);
270:   DMCreateMatrix(da, &J);
271:   TSGetSNES(ts, &snes);
272:   if (use_coloring) {
273:     ISColoring iscoloring;
274:     DMCreateColoring(da, IS_COLORING_GLOBAL, &iscoloring);
275:     MatFDColoringCreate(J, iscoloring, &matfdcoloring);
276:     MatFDColoringSetFromOptions(matfdcoloring);
277:     MatFDColoringSetUp(J, iscoloring, matfdcoloring);
278:     ISColoringDestroy(&iscoloring);
279:     MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, ts);
280:     SNESSetJacobian(snes, J, J, SNESComputeJacobianDefaultColor, matfdcoloring);
281:   } else {
282:     SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault, NULL);
283:   }

285:   /* Define what to print for ts_monitor option */
286:   PetscOptionsHasName(NULL, NULL, "-monitor_off", &monitor_off);
287:   if (!monitor_off) TSMonitorSet(ts, Monitor, &usermonitor, NULL);
288:   FormInitialSolution(da, T, &user);
289:   dt    = TIMESTEP; /* initial time step */
290:   ftime = TIMESTEP * time;
291:   PetscPrintf(PETSC_COMM_WORLD, "time %" PetscInt_FMT ", ftime %g hour, TIMESTEP %g\n", time, (double)(ftime / 3600), (double)dt);

293:   TSSetTimeStep(ts, dt);
294:   TSSetMaxSteps(ts, time);
295:   TSSetMaxTime(ts, ftime);
296:   TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER);
297:   TSSetSolution(ts, T);
298:   TSSetDM(ts, da);

300:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
301:      Set runtime options
302:    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
303:   TSSetFromOptions(ts);

305:   /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
306:      Solve nonlinear system
307:      - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
308:   TSSolve(ts, T);
309:   TSGetSolveTime(ts, &ftime);
310:   TSGetStepNumber(ts, &steps);
311:   PetscPrintf(PETSC_COMM_WORLD, "Solution T after %g hours %" PetscInt_FMT " steps\n", (double)(ftime / 3600), steps);

313:   if (matfdcoloring) MatFDColoringDestroy(&matfdcoloring);
314:   if (usermonitor.drawcontours) PetscViewerDestroy(&usermonitor.drawviewer);
315:   MatDestroy(&J);
316:   VecDestroy(&T);
317:   VecDestroy(&rhs);
318:   TSDestroy(&ts);
319:   DMDestroy(&da);

321:   PetscFinalize();
322:   return 0;
323: }
324: /*****************************end main program********************************/
325: /*****************************************************************************/
326: /*****************************************************************************/
327: /*****************************************************************************/
328: PetscErrorCode calcfluxs(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar fract, PetscScalar cloudTemp, PetscScalar *flux)
329: {
331:   *flux = SIG * ((EMMSFC * emma * PetscPowScalarInt(airtemp, 4)) + (EMMSFC * fract * (1 - emma) * PetscPowScalarInt(cloudTemp, 4)) - (EMMSFC * PetscPowScalarInt(sfctemp, 4))); /* calculates flux using Stefan-Boltzmann relation */
332:   return 0;
333: }

335: PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar *flux) /* this function is not currently called upon */
336: {
337:   PetscScalar emm = 0.001;

340:   *flux = SIG * (-emm * (PetscPowScalarInt(airtemp, 4))); /* calculates flux usinge Stefan-Boltzmann relation */
341:   return 0;
342: }
343: PetscErrorCode sensibleflux(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar wind, PetscScalar *sheat)
344: {
345:   PetscScalar density = 1;    /* air density */
346:   PetscScalar Cp      = 1005; /* heat capicity for dry air */
347:   PetscScalar wndmix;         /* temperature change from wind mixing: wind*Ch */

350:   wndmix = 0.0025 + 0.0042 * wind;                      /* regression equation valid for neutral and stable BL */
351:   *sheat = density * Cp * wndmix * (airtemp - sfctemp); /* calculates sensible heat flux */
352:   return 0;
353: }

355: PetscErrorCode latentflux(PetscScalar sfctemp, PetscScalar dewtemp, PetscScalar wind, PetscScalar pressure1, PetscScalar *latentheat)
356: {
357:   PetscScalar density = 1; /* density of dry air */
358:   PetscScalar q;           /* actual specific humitity */
359:   PetscScalar qs;          /* saturation specific humidity */
360:   PetscScalar wndmix;      /* temperature change from wind mixing: wind*Ch */
361:   PetscScalar beta = .4;   /* moisture availability */
362:   PetscScalar mr;          /* mixing ratio */
363:   PetscScalar lhcnst;      /* latent heat of vaporization constant = 2501000 J/kg at 0c */
364:                            /* latent heat of saturation const = 2834000 J/kg */
365:                            /* latent heat of fusion const = 333700 J/kg */

368:   wind   = mph2mpers(wind);                 /* converts wind from mph to meters per second */
369:   wndmix = 0.0025 + 0.0042 * wind;          /* regression equation valid for neutral BL */
370:   lhcnst = Lconst(sfctemp);                 /* calculates latent heat of evaporation */
371:   mr     = calcmixingr(sfctemp, pressure1); /* calculates saturation mixing ratio */
372:   qs     = calc_q(mr);                      /* calculates saturation specific humidty */
373:   mr     = calcmixingr(dewtemp, pressure1); /* calculates mixing ratio */
374:   q      = calc_q(mr);                      /* calculates specific humidty */

376:   *latentheat = density * wndmix * beta * lhcnst * (q - qs); /* calculates latent heat flux */
377:   return 0;
378: }

380: PetscErrorCode potential_temperature(PetscScalar temp, PetscScalar pressure1, PetscScalar pressure2, PetscScalar sfctemp, PetscScalar *pottemp)
381: {
382:   PetscScalar kdry; /* poisson constant for dry atmosphere */
383:   PetscScalar pavg; /* average atmospheric pressure */
384:   /* PetscScalar mixratio; mixing ratio */
385:   /* PetscScalar kmoist;   poisson constant for moist atmosphere */

388:   /* mixratio = calcmixingr(sfctemp,pressure1); */

390:   /* initialize poisson constant */
391:   kdry = 0.2854;
392:   /* kmoist = 0.2854*(1 - 0.24*mixratio); */

394:   pavg     = ((0.7 * pressure1) + pressure2) / 2;               /* calculates simple average press */
395:   *pottemp = temp * (PetscPowScalar((pressure1 / pavg), kdry)); /* calculates potential temperature */
396:   return 0;
397: }
398: extern PetscScalar calcmixingr(PetscScalar dtemp, PetscScalar pressure1)
399: {
400:   PetscScalar e;        /* vapor pressure */
401:   PetscScalar mixratio; /* mixing ratio */

403:   dtemp    = dtemp - 273;                                                    /* converts from Kelvin to Celsius */
404:   e        = 6.11 * (PetscPowScalar(10, ((7.5 * dtemp) / (237.7 + dtemp)))); /* converts from dew point temp to vapor pressure */
405:   e        = e * 100;                                                        /* converts from hPa to Pa */
406:   mixratio = (0.622 * e) / (pressure1 - e);                                  /* computes mixing ratio */
407:   mixratio = mixratio * 1;                                                   /* convert to g/Kg */

409:   return mixratio;
410: }
411: extern PetscScalar calc_q(PetscScalar rv)
412: {
413:   PetscScalar specific_humidity;     /* define specific humidity variable */
414:   specific_humidity = rv / (1 + rv); /* calculates specific humidity */
415:   return specific_humidity;
416: }

418: PetscErrorCode calc_gflux(PetscScalar sfctemp, PetscScalar deep_grnd_temp, PetscScalar *Gflux)
419: {
420:   PetscScalar k;                       /* thermal conductivity parameter */
421:   PetscScalar n                = 0.38; /* value of soil porosity */
422:   PetscScalar dz               = 1;    /* depth of layer between soil surface and deep soil layer */
423:   PetscScalar unit_soil_weight = 2700; /* unit soil weight in kg/m^3 */

426:   k      = ((0.135 * (1 - n) * unit_soil_weight) + 64.7) / (unit_soil_weight - (0.947 * (1 - n) * unit_soil_weight)); /* dry soil conductivity */
427:   *Gflux = (k * (deep_grnd_temp - sfctemp) / dz);                                                                     /* calculates flux from deep ground layer */
428:   return 0;
429: }
430: extern PetscScalar emission(PetscScalar pwat)
431: {
432:   PetscScalar emma;

434:   emma = 0.725 + 0.17 * PetscLog10Real(PetscRealPart(pwat));

436:   return emma;
437: }
438: extern PetscScalar cloud(PetscScalar fract)
439: {
440:   PetscScalar emma = 0;

442:   /* modifies radiative balance depending on cloud cover */
443:   if (fract >= 0.9) emma = 1;
444:   else if (0.9 > fract && fract >= 0.8) emma = 0.9;
445:   else if (0.8 > fract && fract >= 0.7) emma = 0.85;
446:   else if (0.7 > fract && fract >= 0.6) emma = 0.75;
447:   else if (0.6 > fract && fract >= 0.5) emma = 0.65;
448:   else if (0.4 > fract && fract >= 0.3) emma = emma * 1.086956;
449:   return emma;
450: }
451: extern PetscScalar Lconst(PetscScalar sfctemp)
452: {
453:   PetscScalar Lheat;
454:   sfctemp -= 273;                               /* converts from kelvin to celsius */
455:   Lheat = 4186.8 * (597.31 - 0.5625 * sfctemp); /* calculates latent heat constant */
456:   return Lheat;
457: }
458: extern PetscScalar mph2mpers(PetscScalar wind)
459: {
460:   wind = ((wind * 1.6 * 1000) / 3600); /* converts wind from mph to meters per second */
461:   return wind;
462: }
463: extern PetscScalar fahr_to_cel(PetscScalar temp)
464: {
465:   temp = (5 * (temp - 32)) / 9; /* converts from farhrenheit to celsius */
466:   return temp;
467: }
468: extern PetscScalar cel_to_fahr(PetscScalar temp)
469: {
470:   temp = ((temp * 9) / 5) + 32; /* converts from celsius to farhrenheit */
471:   return temp;
472: }
473: PetscErrorCode readinput(struct in *put)
474: {
475:   int    i;
476:   char   x;
477:   FILE  *ifp;
478:   double tmp;

481:   ifp = fopen("ex5_control.txt", "r");
485:   put->Ts = tmp;

489:   put->Td = tmp;

493:   put->Ta = tmp;

497:   put->Tc = tmp;

501:   put->fr = tmp;

505:   put->wnd = tmp;

509:   put->pwt = tmp;

513:   put->wndDir = tmp;

517:   put->time = tmp;

521:   put->init = tmp;
522:   return 0;
523: }

525: /* ------------------------------------------------------------------- */
526: PetscErrorCode FormInitialSolution(DM da, Vec Xglobal, void *ctx)
527: {
528:   AppCtx  *user = (AppCtx *)ctx; /* user-defined application context */
529:   PetscInt i, j, xs, ys, xm, ym, Mx, My;
530:   Field  **X;

533:   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);

535:   /* Get pointers to vector data */
536:   DMDAVecGetArray(da, Xglobal, &X);

538:   /* Get local grid boundaries */
539:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);

541:   /* Compute function over the locally owned part of the grid */

543:   if (user->init == 1) {
544:     for (j = ys; j < ys + ym; j++) {
545:       for (i = xs; i < xs + xm; i++) {
546:         X[j][i].Ts = user->Ts - i * 0.0001;
547:         X[j][i].Ta = X[j][i].Ts - 5;
548:         X[j][i].u  = 0;
549:         X[j][i].v  = 0;
550:         X[j][i].p  = 1.25;
551:         if ((j == 5 || j == 6) && (i == 4 || i == 5)) X[j][i].p += 0.00001;
552:         if ((j == 5 || j == 6) && (i == 12 || i == 13)) X[j][i].p += 0.00001;
553:       }
554:     }
555:   } else {
556:     for (j = ys; j < ys + ym; j++) {
557:       for (i = xs; i < xs + xm; i++) {
558:         X[j][i].Ts = user->Ts;
559:         X[j][i].Ta = X[j][i].Ts - 5;
560:         X[j][i].u  = 0;
561:         X[j][i].v  = 0;
562:         X[j][i].p  = 1.25;
563:       }
564:     }
565:   }

567:   /* Restore vectors */
568:   DMDAVecRestoreArray(da, Xglobal, &X);
569:   return 0;
570: }

572: /*
573:    RhsFunc - Evaluates nonlinear function F(u).

575:    Input Parameters:
576: .  ts - the TS context
577: .  t - current time
578: .  Xglobal - input vector
579: .  F - output vector
580: .  ptr - optional user-defined context, as set by SNESSetFunction()

582:    Output Parameter:
583: .  F - rhs function vector
584:  */
585: PetscErrorCode RhsFunc(TS ts, PetscReal t, Vec Xglobal, Vec F, void *ctx)
586: {
587:   AppCtx     *user = (AppCtx *)ctx; /* user-defined application context */
588:   DM          da   = user->da;
589:   PetscInt    i, j, Mx, My, xs, ys, xm, ym;
590:   PetscReal   dhx, dhy;
591:   Vec         localT;
592:   Field     **X, **Frhs;                                            /* structures that contain variables of interest and left hand side of governing equations respectively */
593:   PetscScalar csoil          = user->csoil;                         /* heat constant for layer */
594:   PetscScalar dzlay          = user->dzlay;                         /* thickness of top soil layer */
595:   PetscScalar emma           = user->emma;                          /* emission parameter */
596:   PetscScalar wind           = user->wind;                          /* wind speed */
597:   PetscScalar dewtemp        = user->dewtemp;                       /* dew point temperature (moisture in air) */
598:   PetscScalar pressure1      = user->pressure1;                     /* sea level pressure */
599:   PetscScalar airtemp        = user->airtemp;                       /* temperature of air near boundary layer inversion */
600:   PetscScalar fract          = user->fract;                         /* fraction of the sky covered by clouds */
601:   PetscScalar Tc             = user->Tc;                            /* temperature at base of lowest cloud layer */
602:   PetscScalar lat            = user->lat;                           /* latitude */
603:   PetscScalar Cp             = 1005.7;                              /* specific heat of air at constant pressure */
604:   PetscScalar Rd             = 287.058;                             /* gas constant for dry air */
605:   PetscScalar diffconst      = 1000;                                /* diffusion coefficient */
606:   PetscScalar f              = 2 * 0.0000727 * PetscSinScalar(lat); /* coriolis force */
607:   PetscScalar deep_grnd_temp = user->deep_grnd_temp;                /* temp in lowest ground layer */
608:   PetscScalar Ts, u, v, p;
609:   PetscScalar u_abs, u_plus, u_minus, v_abs, v_plus, v_minus;

611:   PetscScalar sfctemp1, fsfc1, Ra;
612:   PetscScalar sheat;      /* sensible heat flux */
613:   PetscScalar latentheat; /* latent heat flux */
614:   PetscScalar groundflux; /* flux from conduction of deep ground layer in contact with top soil */
615:   PetscInt    xend, yend;

618:   DMGetLocalVector(da, &localT);
619:   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);

621:   dhx = (PetscReal)(Mx - 1) / (5000 * (Mx - 1)); /* dhx = 1/dx; assume 2D space domain: [0.0, 1.e5] x [0.0, 1.e5] */
622:   dhy = (PetscReal)(My - 1) / (5000 * (Mx - 1)); /* dhy = 1/dy; */

624:   /*
625:      Scatter ghost points to local vector,using the 2-step process
626:         DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
627:      By placing code between these two statements, computations can be
628:      done while messages are in transition.
629:   */
630:   DMGlobalToLocalBegin(da, Xglobal, INSERT_VALUES, localT);
631:   DMGlobalToLocalEnd(da, Xglobal, INSERT_VALUES, localT);

633:   /* Get pointers to vector data */
634:   DMDAVecGetArrayRead(da, localT, &X);
635:   DMDAVecGetArray(da, F, &Frhs);

637:   /* Get local grid boundaries */
638:   DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL);

640:   /* Compute function over the locally owned part of the grid */
641:   /* the interior points */
642:   xend = xs + xm;
643:   yend = ys + ym;
644:   for (j = ys; j < yend; j++) {
645:     for (i = xs; i < xend; i++) {
646:       Ts = X[j][i].Ts;
647:       u  = X[j][i].u;
648:       v  = X[j][i].v;
649:       p  = X[j][i].p; /*P = X[j][i].P; */

651:       sfctemp1 = (double)Ts;
652:       calcfluxs(sfctemp1, airtemp, emma, fract, Tc, &fsfc1);       /* calculates surface net radiative flux */
653:       sensibleflux(sfctemp1, airtemp, wind, &sheat);               /* calculate sensible heat flux */
654:       latentflux(sfctemp1, dewtemp, wind, pressure1, &latentheat); /* calculates latent heat flux */
655:       calc_gflux(sfctemp1, deep_grnd_temp, &groundflux);           /* calculates flux from earth below surface soil layer by conduction */
656:       calcfluxa(sfctemp1, airtemp, emma, &Ra);                     /* Calculates the change in downward radiative flux */
657:       fsfc1 = fsfc1 + latentheat + sheat + groundflux;                        /* adds radiative, sensible heat, latent heat, and ground heat flux yielding net flux */

659:       /* convective coefficients for upwinding */
660:       u_abs   = PetscAbsScalar(u);
661:       u_plus  = .5 * (u + u_abs); /* u if u>0; 0 if u<0 */
662:       u_minus = .5 * (u - u_abs); /* u if u <0; 0 if u>0 */

664:       v_abs   = PetscAbsScalar(v);
665:       v_plus  = .5 * (v + v_abs); /* v if v>0; 0 if v<0 */
666:       v_minus = .5 * (v - v_abs); /* v if v <0; 0 if v>0 */

668:       /* Solve governing equations */
669:       /* P = p*Rd*Ts; */

671:       /* du/dt -> time change of east-west component of the wind */
672:       Frhs[j][i].u = -u_plus * (u - X[j][i - 1].u) * dhx - u_minus * (X[j][i + 1].u - u) * dhx                                             /* - u(du/dx) */
673:                    - v_plus * (u - X[j - 1][i].u) * dhy - v_minus * (X[j + 1][i].u - u) * dhy                                              /* - v(du/dy) */
674:                    - (Rd / p) * (Ts * (X[j][i + 1].p - X[j][i - 1].p) * 0.5 * dhx + p * 0 * (X[j][i + 1].Ts - X[j][i - 1].Ts) * 0.5 * dhx) /* -(R/p)[Ts(dp/dx)+ p(dTs/dx)] */
675:                                                                                                                                            /*                     -(1/p)*(X[j][i+1].P - X[j][i-1].P)*dhx */
676:                    + f * v;

678:       /* dv/dt -> time change of north-south component of the wind */
679:       Frhs[j][i].v = -u_plus * (v - X[j][i - 1].v) * dhx - u_minus * (X[j][i + 1].v - v) * dhx                                             /* - u(dv/dx) */
680:                    - v_plus * (v - X[j - 1][i].v) * dhy - v_minus * (X[j + 1][i].v - v) * dhy                                              /* - v(dv/dy) */
681:                    - (Rd / p) * (Ts * (X[j + 1][i].p - X[j - 1][i].p) * 0.5 * dhy + p * 0 * (X[j + 1][i].Ts - X[j - 1][i].Ts) * 0.5 * dhy) /* -(R/p)[Ts(dp/dy)+ p(dTs/dy)] */
682:                                                                                                                                            /*                     -(1/p)*(X[j+1][i].P - X[j-1][i].P)*dhy */
683:                    - f * u;

685:       /* dT/dt -> time change of temperature */
686:       Frhs[j][i].Ts = (fsfc1 / (csoil * dzlay))                                                    /* Fnet/(Cp*dz)  diabatic change in T */
687:                     - u_plus * (Ts - X[j][i - 1].Ts) * dhx - u_minus * (X[j][i + 1].Ts - Ts) * dhx /* - u*(dTs/dx)  advection x */
688:                     - v_plus * (Ts - X[j - 1][i].Ts) * dhy - v_minus * (X[j + 1][i].Ts - Ts) * dhy /* - v*(dTs/dy)  advection y */
689:                     + diffconst * ((X[j][i + 1].Ts - 2 * Ts + X[j][i - 1].Ts) * dhx * dhx          /* + D(Ts_xx + Ts_yy)  diffusion */
690:                                    + (X[j + 1][i].Ts - 2 * Ts + X[j - 1][i].Ts) * dhy * dhy);

692:       /* dp/dt -> time change of */
693:       Frhs[j][i].p = -u_plus * (p - X[j][i - 1].p) * dhx - u_minus * (X[j][i + 1].p - p) * dhx /* - u*(dp/dx) */
694:                    - v_plus * (p - X[j - 1][i].p) * dhy - v_minus * (X[j + 1][i].p - p) * dhy; /* - v*(dp/dy) */

696:       Frhs[j][i].Ta = Ra / Cp; /* dTa/dt time change of air temperature */
697:     }
698:   }

700:   /* Restore vectors */
701:   DMDAVecRestoreArrayRead(da, localT, &X);
702:   DMDAVecRestoreArray(da, F, &Frhs);
703:   DMRestoreLocalVector(da, &localT);
704:   return 0;
705: }

707: PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec T, void *ctx)
708: {
709:   const PetscScalar *array;
710:   MonitorCtx        *user   = (MonitorCtx *)ctx;
711:   PetscViewer        viewer = user->drawviewer;
712:   PetscReal          norm;

715:   VecNorm(T, NORM_INFINITY, &norm);

717:   if (step % user->interval == 0) {
718:     VecGetArrayRead(T, &array);
719:     PetscPrintf(PETSC_COMM_WORLD, "step %" PetscInt_FMT ", time %8.1f,  %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f\n", step, (double)time, (double)(((array[0] - 273) * 9) / 5 + 32), (double)(((array[1] - 273) * 9) / 5 + 32), (double)array[2], (double)array[3], (double)array[4], (double)array[5]);
720:     VecRestoreArrayRead(T, &array);
721:   }

723:   if (user->drawcontours) VecView(T, viewer);
724:   return 0;
725: }

727: /*TEST

729:    build:
730:       requires: !complex !single

732:    test:
733:       args: -ts_max_steps 130 -monitor_interval 60
734:       output_file: output/ex5.out
735:       requires: !complex !single
736:       localrunfiles: ex5_control.txt

738:    test:
739:       suffix: 2
740:       nsize: 4
741:       args: -ts_max_steps 130 -monitor_interval 60
742:       output_file: output/ex5.out
743:       localrunfiles: ex5_control.txt
744:       requires: !complex !single

746: TEST*/