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