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