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