1 static char help[] = "Nonlinear, time-dependent. Developed from radiative_surface_balance.c \n";
2 /*
3 Contributed by Steve Froehlich, Illinois Institute of Technology
4
5 Usage:
6 mpiexec -n <np> ./ex5 [options]
7 ./ex5 -help [view petsc options]
8 ./ex5 -ts_type sundials -ts_view
9 ./ex5 -da_grid_x 20 -da_grid_y 20 -log_view
10 ./ex5 -da_grid_x 20 -da_grid_y 20 -ts_type rosw -ts_atol 1.e-6 -ts_rtol 1.e-6
11 ./ex5 -drawcontours -draw_pause 0.1 -draw_fields 0,1,2,3,4
12 */
13
14 /*
15 -----------------------------------------------------------------------
16
17 Governing equations:
18
19 R = s*(Ea*Ta^4 - Es*Ts^4)
20 SH = p*Cp*Ch*wind*(Ta - Ts)
21 LH = p*L*Ch*wind*B(q(Ta) - q(Ts))
22 G = k*(Tgnd - Ts)/dz
23
24 Fnet = R + SH + LH + G
25
26 du/dt = -u*(du/dx) - v*(du/dy) - 2*omeg*sin(lat)*v - (1/p)*(dP/dx)
27 dv/dt = -u*(dv/dx) - v*(dv/dy) + 2*omeg*sin(lat)*u - (1/p)*(dP/dy)
28 dTs/dt = Fnet/(Cp*dz) - Div([u*Ts, v*Ts]) + D*Lap(Ts)
29 = Fnet/(Cs*dz) - u*(dTs/dx) - v*(dTs/dy) + D*(Ts_xx + Ts_yy)
30 dp/dt = -Div([u*p,v*p])
31 = - u*dp/dx - v*dp/dy
32 dTa/dt = Fnet/Cp
33
34 Equation of State:
35
36 P = p*R*Ts
37
38 -----------------------------------------------------------------------
39
40 Program considers the evolution of a two dimensional atmosphere from
41 sunset to sunrise. There are two components:
42 1. Surface energy balance model to compute diabatic dT (Fnet)
43 2. Dynamical model using simplified primitive equations
44
45 Program is to be initiated at sunset and run to sunrise.
46
47 Inputs are:
48 Surface temperature
49 Dew point temperature
50 Air temperature
51 Temperature at cloud base (if clouds are present)
52 Fraction of sky covered by clouds
53 Wind speed
54 Precipitable water in centimeters
55 Wind direction
56
57 Inputs are are read in from the text file ex5_control.txt. To change an
58 input value use ex5_control.txt.
59
60 Solvers:
61 Backward Euler = default solver
62 Sundials = fastest and most accurate, requires Sundials libraries
63
64 This model is under development and should be used only as an example
65 and not as a predictive weather model.
66 */
67
68 #include <petscts.h>
69 #include <petscdm.h>
70 #include <petscdmda.h>
71
72 /* stefan-boltzmann constant */
73 #define SIG 0.000000056703
74 /* absorption-emission constant for surface */
75 #define EMMSFC 1
76 /* amount of time (seconds) that passes before new flux is calculated */
77 #define TIMESTEP 1
78
79 /* variables of interest to be solved at each grid point */
80 typedef struct {
81 PetscScalar Ts,Ta; /* surface and air temperature */
82 PetscScalar u,v; /* wind speed */
83 PetscScalar p; /* density */
84 } Field;
85
86 /* User defined variables. Used in solving for variables of interest */
87 typedef struct {
88 DM da; /* grid */
89 PetscScalar csoil; /* heat constant for layer */
90 PetscScalar dzlay; /* thickness of top soil layer */
91 PetscScalar emma; /* emission parameter */
92 PetscScalar wind; /* wind speed */
93 PetscScalar dewtemp; /* dew point temperature (moisture in air) */
94 PetscScalar pressure1; /* sea level pressure */
95 PetscScalar airtemp; /* temperature of air near boundary layer inversion */
96 PetscScalar Ts; /* temperature at the surface */
97 PetscScalar fract; /* fraction of sky covered by clouds */
98 PetscScalar Tc; /* temperature at base of lowest cloud layer */
99 PetscScalar lat; /* Latitude in degrees */
100 PetscScalar init; /* initialization scenario */
101 PetscScalar deep_grnd_temp; /* temperature of ground under top soil surface layer */
102 } AppCtx;
103
104 /* Struct for visualization */
105 typedef struct {
106 PetscBool drawcontours; /* flag - 1 indicates drawing contours */
107 PetscViewer drawviewer;
108 PetscInt interval;
109 } MonitorCtx;
110
111
112 /* Inputs read in from text file */
113 struct in {
114 PetscScalar Ts; /* surface temperature */
115 PetscScalar Td; /* dewpoint temperature */
116 PetscScalar Tc; /* temperature of cloud base */
117 PetscScalar fr; /* fraction of sky covered by clouds */
118 PetscScalar wnd; /* wind speed */
119 PetscScalar Ta; /* air temperature */
120 PetscScalar pwt; /* precipitable water */
121 PetscScalar wndDir; /* wind direction */
122 PetscScalar lat; /* latitude */
123 PetscReal time; /* time in hours */
124 PetscScalar init;
125 };
126
127 /* functions */
128 extern PetscScalar emission(PetscScalar); /* sets emission/absorption constant depending on water vapor content */
129 extern PetscScalar calc_q(PetscScalar); /* calculates specific humidity */
130 extern PetscScalar mph2mpers(PetscScalar); /* converts miles per hour to meters per second */
131 extern PetscScalar Lconst(PetscScalar); /* calculates latent heat constant taken from Satellite estimates of wind speed and latent heat flux over the global oceans., Bentamy et al. */
132 extern PetscScalar fahr_to_cel(PetscScalar); /* converts Fahrenheit to Celsius */
133 extern PetscScalar cel_to_fahr(PetscScalar); /* converts Celsius to Fahrenheit */
134 extern PetscScalar calcmixingr(PetscScalar, PetscScalar); /* calculates mixing ratio */
135 extern PetscScalar cloud(PetscScalar); /* cloud radiative parameterization */
136 extern PetscErrorCode FormInitialSolution(DM,Vec,void*); /* Specifies initial conditions for the system of equations (PETSc defined function) */
137 extern PetscErrorCode RhsFunc(TS,PetscReal,Vec,Vec,void*); /* Specifies the user defined functions (PETSc defined function) */
138 extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*); /* Specifies output and visualization tools (PETSc defined function) */
139 extern 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 */
146
main(int argc,char ** argv)147 int main(int argc,char **argv)
148 {
149 PetscErrorCode ierr;
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;
175
176 ierr = PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
177
178 /* Inputs */
179 ierr = readinput(&put);CHKERRQ(ierr);
180
181 sfctemp = put.Ts;
182 dewtemp = put.Td;
183 cloudTemp = put.Tc;
184 airtemp = put.Ta;
185 pwat = put.pwt;
186
187 ierr = PetscPrintf(PETSC_COMM_WORLD,"Initial Temperature = %g\n",(double)sfctemp);CHKERRQ(ierr); /* input surface temperature */
188
189 deep_grnd_temp = sfctemp - 10; /* set underlying ground layer temperature */
190 emma = emission(pwat); /* accounts for radiative effects of water vapor */
191
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);
198
199 /* Converts from Celsius to Kelvin */
200 sfctemp += 273;
201 airtemp += 273;
202 dewtemp += 273;
203 cloudTemp += 273;
204 deep_grnd_temp += 273;
205
206 /* Calculates initial relative humidity */
207 x = calcmixingr(dewtemp,pressure1);
208 mixratio = calcmixingr(sfctemp,pressure1);
209 rh = (x/mixratio)*100;
210
211 ierr = PetscPrintf(PETSC_COMM_WORLD,"Initial RH = %.1f percent\n\n",(double)rh);CHKERRQ(ierr); /* prints initial relative humidity */
212
213 time = 3600*put.time; /* sets amount of timesteps to run model */
214
215 /* Configure PETSc TS solver */
216 /*------------------------------------------*/
217
218 /* Create grid */
219 ierr = DMDACreate2d(PETSC_COMM_WORLD,DM_BOUNDARY_PERIODIC,DM_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,20,20,PETSC_DECIDE,PETSC_DECIDE,dof,1,NULL,NULL,&da);CHKERRQ(ierr);
220 ierr = DMSetFromOptions(da);CHKERRQ(ierr);
221 ierr = DMSetUp(da);CHKERRQ(ierr);
222 ierr = DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0);CHKERRQ(ierr);
223
224 /* Define output window for each variable of interest */
225 ierr = DMDASetFieldName(da,0,"Ts");CHKERRQ(ierr);
226 ierr = DMDASetFieldName(da,1,"Ta");CHKERRQ(ierr);
227 ierr = DMDASetFieldName(da,2,"u");CHKERRQ(ierr);
228 ierr = DMDASetFieldName(da,3,"v");CHKERRQ(ierr);
229 ierr = DMDASetFieldName(da,4,"p");CHKERRQ(ierr);
230
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 */
246
247 /* set values for MonitorCtx */
248 usermonitor.drawcontours = PETSC_FALSE;
249 ierr = PetscOptionsHasName(NULL,NULL,"-drawcontours",&usermonitor.drawcontours);CHKERRQ(ierr);
250 if (usermonitor.drawcontours) {
251 PetscReal bounds[] = {1000.0,-1000., -1000.,-1000., 1000.,-1000., 1000.,-1000., 1000,-1000, 100700,100800};
252 ierr = PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,300,300,&usermonitor.drawviewer);CHKERRQ(ierr);
253 ierr = PetscViewerDrawSetBounds(usermonitor.drawviewer,dof,bounds);CHKERRQ(ierr);
254 }
255 usermonitor.interval = 1;
256 ierr = PetscOptionsGetInt(NULL,NULL,"-monitor_interval",&usermonitor.interval,NULL);CHKERRQ(ierr);
257
258 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
259 Extract global vectors from DA;
260 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
261 ierr = DMCreateGlobalVector(da,&T);CHKERRQ(ierr);
262 ierr = VecDuplicate(T,&rhs);CHKERRQ(ierr); /* r: vector to put the computed right hand side */
263
264 ierr = TSCreate(PETSC_COMM_WORLD,&ts);CHKERRQ(ierr);
265 ierr = TSSetProblemType(ts,TS_NONLINEAR);CHKERRQ(ierr);
266 ierr = TSSetType(ts,TSBEULER);CHKERRQ(ierr);
267 ierr = TSSetRHSFunction(ts,rhs,RhsFunc,&user);CHKERRQ(ierr);
268
269 /* Set Jacobian evaluation routine - use coloring to compute finite difference Jacobian efficiently */
270 ierr = DMSetMatType(da,MATAIJ);CHKERRQ(ierr);
271 ierr = DMCreateMatrix(da,&J);CHKERRQ(ierr);
272 ierr = TSGetSNES(ts,&snes);CHKERRQ(ierr);
273 if (use_coloring) {
274 ISColoring iscoloring;
275 ierr = DMCreateColoring(da,IS_COLORING_GLOBAL,&iscoloring);CHKERRQ(ierr);
276 ierr = MatFDColoringCreate(J,iscoloring,&matfdcoloring);CHKERRQ(ierr);
277 ierr = MatFDColoringSetFromOptions(matfdcoloring);CHKERRQ(ierr);
278 ierr = MatFDColoringSetUp(J,iscoloring,matfdcoloring);CHKERRQ(ierr);
279 ierr = ISColoringDestroy(&iscoloring);CHKERRQ(ierr);
280 ierr = MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);CHKERRQ(ierr);
281 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefaultColor,matfdcoloring);CHKERRQ(ierr);
282 } else {
283 ierr = SNESSetJacobian(snes,J,J,SNESComputeJacobianDefault,NULL);CHKERRQ(ierr);
284 }
285
286 /* Define what to print for ts_monitor option */
287 ierr = PetscOptionsHasName(NULL,NULL,"-monitor_off",&monitor_off);CHKERRQ(ierr);
288 if (!monitor_off) {
289 ierr = TSMonitorSet(ts,Monitor,&usermonitor,NULL);CHKERRQ(ierr);
290 }
291 ierr = FormInitialSolution(da,T,&user);CHKERRQ(ierr);
292 dt = TIMESTEP; /* initial time step */
293 ftime = TIMESTEP*time;
294 ierr = PetscPrintf(PETSC_COMM_WORLD,"time %D, ftime %g hour, TIMESTEP %g\n",time,(double)(ftime/3600),(double)dt);CHKERRQ(ierr);
295
296 ierr = TSSetTimeStep(ts,dt);CHKERRQ(ierr);
297 ierr = TSSetMaxSteps(ts,time);CHKERRQ(ierr);
298 ierr = TSSetMaxTime(ts,ftime);CHKERRQ(ierr);
299 ierr = TSSetExactFinalTime(ts,TS_EXACTFINALTIME_STEPOVER);CHKERRQ(ierr);
300 ierr = TSSetSolution(ts,T);CHKERRQ(ierr);
301 ierr = TSSetDM(ts,da);CHKERRQ(ierr);
302
303 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
304 Set runtime options
305 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
306 ierr = TSSetFromOptions(ts);CHKERRQ(ierr);
307
308 /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
309 Solve nonlinear system
310 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
311 ierr = TSSolve(ts,T);CHKERRQ(ierr);
312 ierr = TSGetSolveTime(ts,&ftime);CHKERRQ(ierr);
313 ierr = TSGetStepNumber(ts,&steps);CHKERRQ(ierr);
314 ierr = PetscPrintf(PETSC_COMM_WORLD,"Solution T after %g hours %D steps\n",(double)(ftime/3600),steps);CHKERRQ(ierr);
315
316
317 if (matfdcoloring) {ierr = MatFDColoringDestroy(&matfdcoloring);CHKERRQ(ierr);}
318 if (usermonitor.drawcontours) {
319 ierr = PetscViewerDestroy(&usermonitor.drawviewer);CHKERRQ(ierr);
320 }
321 ierr = MatDestroy(&J);CHKERRQ(ierr);
322 ierr = VecDestroy(&T);CHKERRQ(ierr);
323 ierr = VecDestroy(&rhs);CHKERRQ(ierr);
324 ierr = TSDestroy(&ts);CHKERRQ(ierr);
325 ierr = DMDestroy(&da);CHKERRQ(ierr);
326
327 ierr = PetscFinalize();
328 return ierr;
329 }
330 /*****************************end main program********************************/
331 /*****************************************************************************/
332 /*****************************************************************************/
333 /*****************************************************************************/
calcfluxs(PetscScalar sfctemp,PetscScalar airtemp,PetscScalar emma,PetscScalar fract,PetscScalar cloudTemp,PetscScalar * flux)334 PetscErrorCode calcfluxs(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar fract, PetscScalar cloudTemp, PetscScalar *flux)
335 {
336 PetscFunctionBeginUser;
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 PetscFunctionReturn(0);
339 }
340
calcfluxa(PetscScalar sfctemp,PetscScalar airtemp,PetscScalar emma,PetscScalar * flux)341 PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar *flux) /* this function is not currently called upon */
342 {
343 PetscScalar emm = 0.001;
344
345 PetscFunctionBeginUser;
346 *flux = SIG*(-emm*(PetscPowScalarInt(airtemp,4))); /* calculates flux usinge Stefan-Boltzmann relation */
347 PetscFunctionReturn(0);
348 }
sensibleflux(PetscScalar sfctemp,PetscScalar airtemp,PetscScalar wind,PetscScalar * sheat)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 */
354
355 PetscFunctionBeginUser;
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 PetscFunctionReturn(0);
359 }
360
latentflux(PetscScalar sfctemp,PetscScalar dewtemp,PetscScalar wind,PetscScalar pressure1,PetscScalar * latentheat)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 */
372
373 PetscFunctionBeginUser;
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 */
381
382 *latentheat = density*wndmix*beta*lhcnst*(q - qs); /* calculates latent heat flux */
383 PetscFunctionReturn(0);
384 }
385
potential_temperature(PetscScalar temp,PetscScalar pressure1,PetscScalar pressure2,PetscScalar sfctemp,PetscScalar * pottemp)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 */
392
393 PetscFunctionBeginUser;
394 /* mixratio = calcmixingr(sfctemp,pressure1); */
395
396 /* initialize poisson constant */
397 kdry = 0.2854;
398 /* kmoist = 0.2854*(1 - 0.24*mixratio); */
399
400 pavg = ((0.7*pressure1)+pressure2)/2; /* calculates simple average press */
401 *pottemp = temp*(PetscPowScalar((pressure1/pavg),kdry)); /* calculates potential temperature */
402 PetscFunctionReturn(0);
403 }
calcmixingr(PetscScalar dtemp,PetscScalar pressure1)404 extern PetscScalar calcmixingr(PetscScalar dtemp, PetscScalar pressure1)
405 {
406 PetscScalar e; /* vapor pressure */
407 PetscScalar mixratio; /* mixing ratio */
408
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 */
414
415 return mixratio;
416 }
calc_q(PetscScalar rv)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 }
423
calc_gflux(PetscScalar sfctemp,PetscScalar deep_grnd_temp,PetscScalar * Gflux)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 */
430
431 PetscFunctionBeginUser;
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 PetscFunctionReturn(0);
435 }
emission(PetscScalar pwat)436 extern PetscScalar emission(PetscScalar pwat)
437 {
438 PetscScalar emma;
439
440 emma = 0.725 + 0.17*PetscLog10Real(PetscRealPart(pwat));
441
442 return emma;
443 }
cloud(PetscScalar fract)444 extern PetscScalar cloud(PetscScalar fract)
445 {
446 PetscScalar emma = 0;
447
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 }
Lconst(PetscScalar sfctemp)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 }
mph2mpers(PetscScalar wind)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 }
fahr_to_cel(PetscScalar temp)469 extern PetscScalar fahr_to_cel(PetscScalar temp)
470 {
471 temp = (5*(temp-32))/9; /* converts from farhrenheit to celsuis */
472 return temp;
473 }
cel_to_fahr(PetscScalar temp)474 extern PetscScalar cel_to_fahr(PetscScalar temp)
475 {
476 temp = ((temp*9)/5) + 32; /* converts from celsuis to farhrenheit */
477 return temp;
478 }
readinput(struct in * put)479 PetscErrorCode readinput(struct in *put)
480 {
481 int i;
482 char x;
483 FILE *ifp;
484 double tmp;
485
486 PetscFunctionBegin;
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;
492
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;
496
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;
500
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;
504
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;
508
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;
512
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;
516
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;
520
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;
524
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 PetscFunctionReturn(0);
529 }
530
531 /* ------------------------------------------------------------------- */
FormInitialSolution(DM da,Vec Xglobal,void * ctx)532 PetscErrorCode FormInitialSolution(DM da,Vec Xglobal,void *ctx)
533 {
534 PetscErrorCode ierr;
535 AppCtx *user = (AppCtx*)ctx; /* user-defined application context */
536 PetscInt i,j,xs,ys,xm,ym,Mx,My;
537 Field **X;
538
539 PetscFunctionBeginUser;
540 ierr = 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);CHKERRQ(ierr);
542
543 /* Get pointers to vector data */
544 ierr = DMDAVecGetArray(da,Xglobal,&X);CHKERRQ(ierr);
545
546 /* Get local grid boundaries */
547 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
548
549 /* Compute function over the locally owned part of the grid */
550
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 }
574
575 /* Restore vectors */
576 ierr = DMDAVecRestoreArray(da,Xglobal,&X);CHKERRQ(ierr);
577 PetscFunctionReturn(0);
578 }
579
580 /*
581 RhsFunc - Evaluates nonlinear function F(u).
582
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()
589
590 Output Parameter:
591 . F - rhs function vector
592 */
RhsFunc(TS ts,PetscReal t,Vec Xglobal,Vec F,void * ctx)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;
597 PetscErrorCode ierr;
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;
619
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;
625
626 PetscFunctionBeginUser;
627 ierr = DMGetLocalVector(da,&localT);CHKERRQ(ierr);
628 ierr = DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);CHKERRQ(ierr);
629
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; */
632
633
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 ierr = DMGlobalToLocalBegin(da,Xglobal,INSERT_VALUES,localT);CHKERRQ(ierr);
641 ierr = DMGlobalToLocalEnd(da,Xglobal,INSERT_VALUES,localT);CHKERRQ(ierr);
642
643 /* Get pointers to vector data */
644 ierr = DMDAVecGetArrayRead(da,localT,&X);CHKERRQ(ierr);
645 ierr = DMDAVecGetArray(da,F,&Frhs);CHKERRQ(ierr);
646
647 /* Get local grid boundaries */
648 ierr = DMDAGetCorners(da,&xs,&ys,NULL,&xm,&ym,NULL);CHKERRQ(ierr);
649
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; */
656
657 sfctemp1 = (double)Ts;
658 ierr = calcfluxs(sfctemp1,airtemp,emma,fract,Tc,&fsfc1);CHKERRQ(ierr); /* calculates surface net radiative flux */
659 ierr = sensibleflux(sfctemp1,airtemp,wind,&sheat);CHKERRQ(ierr); /* calculate sensible heat flux */
660 ierr = latentflux(sfctemp1,dewtemp,wind,pressure1,&latentheat);CHKERRQ(ierr); /* calculates latent heat flux */
661 ierr = calc_gflux(sfctemp1,deep_grnd_temp,&groundflux);CHKERRQ(ierr); /* calculates flux from earth below surface soil layer by conduction */
662 ierr = calcfluxa(sfctemp1,airtemp,emma,&Ra);CHKERRQ(ierr); /* 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 */
664
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 */
669
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 */
673
674 /* Solve governing equations */
675 /* P = p*Rd*Ts; */
676
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;
683
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;
690
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);
697
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) */
701
702 Frhs[j][i].Ta = Ra/Cp; /* dTa/dt time change of air temperature */
703 }
704 }
705
706 /* Restore vectors */
707 ierr = DMDAVecRestoreArrayRead(da,localT,&X);CHKERRQ(ierr);
708 ierr = DMDAVecRestoreArray(da,F,&Frhs);CHKERRQ(ierr);
709 ierr = DMRestoreLocalVector(da,&localT);CHKERRQ(ierr);
710 PetscFunctionReturn(0);
711 }
712
Monitor(TS ts,PetscInt step,PetscReal time,Vec T,void * ctx)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;
720
721 PetscFunctionBeginUser;
722 ierr = VecNorm(T,NORM_INFINITY,&norm);CHKERRQ(ierr);
723
724 if (step%user->interval == 0) {
725 ierr = VecGetArrayRead(T,&array);CHKERRQ(ierr);
726 ierr = 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 ierr = VecRestoreArrayRead(T,&array);CHKERRQ(ierr);
728 }
729
730 if (user->drawcontours) {
731 ierr = VecView(T,viewer);CHKERRQ(ierr);
732 }
733 PetscFunctionReturn(0);
734 }
735
736
737
738 /*TEST
739
740 build:
741 requires: !complex !single
742
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
748
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
756
757 TEST*/
758