Actual source code: ex5.c
petsc-3.3-p6 2013-02-11
1: static char help[] = "Nonlinear, time-dependent. Developed from radiative_surface_balance.c - work with petsc-dev \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_sundials_monitor_steps -pc_type lu -ts_view
9: ./ex5 -da_grid_x 20 -da_grid_y 20 -log_summary
10: ./ex5 -drawcontours -draw_pause 0.1 -draw_fields 0,1,2,3,4
11: */
13: /*
14: -----------------------------------------------------------------------
16: Governing equations:
18: R = s*(Ea*Ta^4 - Es*Ts^4)
19: SH = p*Cp*Ch*wind*(Ta - Ts)
20: LH = p*L*Ch*wind*B(q(Ta) - q(Ts))
21: G = k*(Tgnd - Ts)/dz
23: Fnet = R + SH + LH + G
25: du/dt = -u*(du/dx) - v*(du/dy) - 2*omeg*sin(lat)*v - (1/p)*(dP/dx)
26: dv/dt = -u*(dv/dx) - v*(dv/dy) + 2*omeg*sin(lat)*u - (1/p)*(dP/dy)
27: dTs/dt = Fnet/(Cp*dz) - Div([u*Ts, v*Ts]) + D*Lap(Ts)
28: = Fnet/(Cs*dz) - u*(dTs/dx) - v*(dTs/dy) + D*(Ts_xx + Ts_yy)
29: dp/dt = -Div([u*p,v*p])
30: = - u*dp/dx - v*dp/dy
31: dTa/dt = Fnet/Cp
33: Equation of State:
35: P = p*R*Ts
37: -----------------------------------------------------------------------
39: Program considers the evolution of a two dimensional atmosphere from
40: sunset to sunrise. There are two components:
41: 1. Surface energy balance model to compute diabatic dT (Fnet)
42: 2. Dynamical model using simplified primitive equations
44: Program is to be initiated at sunset and run to sunrise.
46: Inputs are:
47: Surface temperature
48: Dew point temperature
49: Air temperature
50: Temperature at cloud base (if clouds are present)
51: Fraction of sky covered by clouds
52: Wind speed
53: Precipitable water in centimeters
54: Wind direction
56: Inputs are are read in from the text file ex5_control.txt. To change an
57: input value use ex5_control.txt.
59: Solvers:
60: Backward Euler = default solver
61: Sundials = fastest and most accurate, requires Sundials libraries
63: This model is under development and should be used only as an example
64: and not as a predictive weather model.
65: */
67: #include <petscts.h>
68: #include <petscdmda.h>
70: #define SIG 0.000000056703 //stefan-boltzmann constant
71: #define EMMSFC 1 //absorption-emission constant for surface
72: #define TIMESTEP 1 //amount of time(seconds) that passes before new flux is calculated
74: /* variables of interest to be solved at each grid point */
75: typedef struct {
76: PetscScalar Ts,Ta; // surface and air temperature
77: PetscScalar u,v; // wind speed
78: PetscScalar p; // density
79: } Field;
81: /* User defined variables. Used in solving for variables of interest */
82: typedef struct {
83: DM da; //grid
84: PetscScalar csoil; //heat constant for layer
85: PetscScalar dzlay; //thickness of top soil layer
86: PetscScalar emma; //emission parameter
87: PetscScalar wind; //wind speed
88: PetscScalar dewtemp; //dew point temperature (moisture in air)
89: PetscScalar pressure1; //sea level pressure
90: PetscScalar airtemp; //temperature of air near boundary layer inversion
91: PetscScalar Ts; //temperature at the surface
92: PetscScalar fract; //fraction of sky covered by clouds
93: PetscScalar Tc; //temperature at base of lowest cloud layer
94: PetscScalar lat; //Latitude in degrees
95: PetscScalar init; //initialization scenario
96: PetscScalar deep_grnd_temp;//temperature of ground under top soil surface layer
97: } AppCtx;
99: /* Struct for visualization */
100: typedef struct {
101: PetscBool drawcontours; /* flag - 1 indicates drawing contours */
102: PetscViewer drawviewer;
103: } MonitorCtx;
106: /* Inputs read in from text file */
107: struct in {
108: PetscScalar Ts; //surface temperature
109: PetscScalar Td; //dewpoint temperature
110: PetscScalar Tc; //temperature of cloud base
111: PetscScalar fr; //fraction of sky covered by clouds
112: PetscScalar wnd; //wind speed
113: PetscScalar Ta; //air temperature
114: PetscScalar pwt; //precipitable water
115: PetscScalar wndDir; //wind direction
116: PetscScalar lat; //latitude
117: PetscReal time; //time in hours
118: PetscScalar init;
119: };
121: //functions
122: extern PetscScalar emission(PetscScalar); //sets emission/absorption constant depending on water vapor content
123: extern PetscScalar calc_q(PetscScalar); //calculates specific humidity
124: extern PetscScalar mph2mpers(PetscScalar); //converts miles per hour to meters per second
125: 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.
126: extern PetscScalar fahr_to_cel(PetscScalar); //converts Fahrenheit to Celsius
127: extern PetscScalar cel_to_fahr(PetscScalar); //converts Celsius to Fahrenheit
128: extern PetscScalar calcmixingr(PetscScalar, PetscScalar); //calculates mixing ratio
129: extern PetscScalar cloud(PetscScalar); //cloud radiative parameterization
130: extern PetscErrorCode FormInitialSolution(DM,Vec,void*); //Specifies initial conditions for the system of equations (PETSc defined function)
131: extern PetscErrorCode RhsFunc(TS,PetscReal,Vec,Vec,void*); //Specifies the user defined functions (PETSc defined function)
132: extern PetscErrorCode Monitor(TS,PetscInt,PetscReal,Vec,void*); //Specifies output and visualization tools (PETSc defined function)
133: extern void readinput(struct in *put); //reads input from text file
134: extern PetscErrorCode calcfluxs(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*);//calculates upward IR from surface
135: extern PetscErrorCode calcfluxa(PetscScalar, PetscScalar, PetscScalar, PetscScalar*); //calculates downward IR from atmosphere
136: extern PetscErrorCode sensibleflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar*); //calculates sensible heat flux
137: extern PetscErrorCode potential_temperature(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*); //calculates potential temperature
138: extern PetscErrorCode latentflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar*); //calculates latent heat flux
139: extern PetscErrorCode calc_gflux(PetscScalar, PetscScalar, PetscScalar*); //calculates flux between top soil layer and underlying earth
143: int main(int argc,char **argv)
144: {
145: PetscErrorCode ierr;
146: int time; //amount of loops
147: struct in put;
148: PetscScalar rh; //relative humidity
149: PetscScalar x; //memory varialbe for relative humidity calculation
150: PetscScalar deep_grnd_temp; //temperature of ground under top soil surface layer
152: PetscScalar emma; //absorption-emission constant for air
153: PetscScalar pressure1 = 101300; //surface pressure
154: PetscScalar mixratio; //mixing ratio
155: PetscScalar airtemp; //temperature of air near boundary layer inversion
156: PetscScalar dewtemp; //dew point temperature
157: PetscScalar sfctemp; //temperature at surface
158: PetscScalar pwat; //total column precipitable water
159: PetscScalar cloudTemp; //temperature at base of cloud
160: AppCtx user; /* user-defined work context */
161: MonitorCtx usermonitor; /* user-defined monitor context */
162: PetscMPIInt rank,size;
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;
171: PetscInitialize(&argc,&argv,(char *)0,help);
172: MPI_Comm_size(PETSC_COMM_WORLD,&size);
173: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
175: /* Inputs */
176: readinput(&put);
178: sfctemp = put.Ts;
179: dewtemp = put.Td;
180: cloudTemp = put.Tc;
181: airtemp = put.Ta;
182: pwat = put.pwt;
184: if (!rank) PetscPrintf(PETSC_COMM_SELF,"Initial Temperature = %g\n",sfctemp); //input surface temperature
186: deep_grnd_temp = sfctemp - 10; //set underlying ground layer temperature
187: emma = emission(pwat); //accounts for radiative effects of water vapor
189: /* Converts from Fahrenheit to Celsuis */
190: sfctemp = fahr_to_cel(sfctemp);
191: airtemp = fahr_to_cel(airtemp);
192: dewtemp = fahr_to_cel(dewtemp);
193: cloudTemp = fahr_to_cel(cloudTemp);
194: deep_grnd_temp = fahr_to_cel(deep_grnd_temp);
196: /* Converts from Celsius to Kelvin */
197: sfctemp +=273;
198: airtemp +=273;
199: dewtemp +=273;
200: cloudTemp +=273;
201: deep_grnd_temp +=273;
203: /* Calculates initial relative humidity */
204: x = calcmixingr(dewtemp,pressure1);
205: mixratio = calcmixingr(sfctemp,pressure1);
206: rh = (x/mixratio)*100;
208: if (!rank){printf("Initial RH = %.1f percent\n\n",rh);} //prints initial relative humidity
210: time = 3600*put.time; //sets amount of timesteps to run model
212: /* Configure PETSc TS solver */
213: /*------------------------------------------*/
215: /*Create grid*/
216: DMDACreate2d(PETSC_COMM_WORLD,DMDA_BOUNDARY_PERIODIC,DMDA_BOUNDARY_PERIODIC,DMDA_STENCIL_STAR,-20,-20,
217: PETSC_DECIDE,PETSC_DECIDE,dof,1,PETSC_NULL,PETSC_NULL,&da);
219: /*Define output window for each variable of interest*/
220: DMDASetFieldName(da,0,"Ts");
221: DMDASetFieldName(da,1,"Ta");
222: DMDASetFieldName(da,2,"u");
223: DMDASetFieldName(da,3,"v");
224: DMDASetFieldName(da,4,"p");
226: /*set values for appctx*/
227: user.da = da;
228: user.Ts = sfctemp;
229: user.fract = put.fr; //fraction of sky covered by clouds
230: user.dewtemp = dewtemp; //dew point temperature (mositure in air)
231: user.csoil = 2000000; //heat constant for layer
232: user.dzlay = 0.08; //thickness of top soil layer
233: user.emma = emma; //emission parameter
234: user.wind = put.wnd; //wind spped
235: user.pressure1 = pressure1; //sea level pressure
236: user.airtemp = airtemp; //temperature of air near boundar layer inversion
237: user.Tc = cloudTemp; //temperature at base of lowest cloud layer
238: user.init = put.init; //user chosen initiation scenario
239: user.lat = 70*0.0174532; //converts latitude degrees to latitude in radians
240: user.deep_grnd_temp = deep_grnd_temp; //temp in lowest ground layer
242: /*set values for MonitorCtx*/
243: usermonitor.drawcontours = PETSC_FALSE;
244: PetscOptionsHasName(PETSC_NULL,"-drawcontours",&usermonitor.drawcontours);
245: if (usermonitor.drawcontours){
246: PetscReal bounds[] = {1000.0,-1000., -1000.,-1000., 1000.,-1000., 1000.,-1000., 1000,-1000, 100700,100800};
247: PetscViewerDrawOpen(PETSC_COMM_WORLD,0,0,0,0,300,300,&usermonitor.drawviewer);
248: PetscViewerDrawSetBounds(usermonitor.drawviewer,dof,bounds);
249: }
251: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
252: Extract global vectors from DA;
253: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
254: DMCreateGlobalVector(da,&T);
255: VecDuplicate(T,&rhs); //r: vector to put the computed right hand side
257: TSCreate(PETSC_COMM_WORLD,&ts);
258: TSSetProblemType(ts,TS_NONLINEAR);
259: TSSetType(ts,TSBEULER);
260: TSSetRHSFunction(ts,rhs,RhsFunc,&user);
262: /* Set Jacobian evaluation routine - use coloring to compute finite difference Jacobian efficiently */
263: PetscBool use_coloring=PETSC_TRUE;
264: MatFDColoring matfdcoloring=0;
265: DMCreateMatrix(da,MATAIJ,&J);
266: TSGetSNES(ts,&snes);
267: if (use_coloring){
268: ISColoring iscoloring;
269: DMCreateColoring(da,IS_COLORING_GLOBAL,MATAIJ,&iscoloring);
270: MatFDColoringCreate(J,iscoloring,&matfdcoloring);
271: MatFDColoringSetFromOptions(matfdcoloring);
272: ISColoringDestroy(&iscoloring);
273: MatFDColoringSetFunction(matfdcoloring,(PetscErrorCode (*)(void))SNESTSFormFunction,ts);
274: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobianColor,matfdcoloring);
275: } else {
276: SNESSetJacobian(snes,J,J,SNESDefaultComputeJacobian,PETSC_NULL);
277: }
279: /*Define what to print for ts_monitor option*/
280: TSMonitorSet(ts,Monitor,&usermonitor,PETSC_NULL);
281: FormInitialSolution(da,T,&user);
282: dt = TIMESTEP; /* initial time step */
283: ftime = TIMESTEP*time;
284: if (!rank){printf("time %d, ftime %g hour, TIMESTEP %g\n",time,ftime/3600,dt);}
285: TSSetInitialTimeStep(ts,0.0,dt);
286: TSSetDuration(ts,time,ftime);
287: TSSetSolution(ts,T);
289: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
290: Set runtime options
291: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
292: TSSetFromOptions(ts);
294: /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
295: Solve nonlinear system
296: - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */
297: TSSolve(ts,T,&ftime);
298: TSGetTimeStepNumber(ts,&steps);
299: if (!rank){PetscPrintf(PETSC_COMM_WORLD,"Solution T after %g hours %d steps\n",ftime/3600,steps);}
302: if (matfdcoloring){MatFDColoringDestroy(&matfdcoloring);}
303: if (usermonitor.drawcontours){
304: PetscViewerDestroy(&usermonitor.drawviewer);
305: }
306: MatDestroy(&J);
307: VecDestroy(&T);
308: VecDestroy(&rhs);
309: TSDestroy(&ts);
310: DMDestroy(&da);
312: PetscFinalize();
313: return 0;
314: }
315: /*****************************end main program********************************/
316: /*****************************************************************************/
317: /*****************************************************************************/
318: /*****************************************************************************/
321: PetscErrorCode calcfluxs(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar fract, PetscScalar cloudTemp, PetscScalar* flux)
322: {
324: *flux = SIG*((EMMSFC*emma*pow(airtemp,4)) + (EMMSFC*fract*(1 - emma)*pow(cloudTemp,4)) - (EMMSFC*pow(sfctemp,4))); //calculates flux using Stefan-Boltzmann relation
325: return(0);
326: }
330: PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar* flux) //this function is not currently called upon
331: {
332: PetscScalar emm = 0.001;
334: *flux = SIG*(- emm*(pow(airtemp,4))); //calculates flux usinge Stefan-Boltzmann relation
335: return(0);
336: }
339: PetscErrorCode sensibleflux(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar wind, PetscScalar* sheat)
340: {
341: PetscScalar density = 1; //air density
342: PetscScalar Cp = 1005; //heat capicity for dry air
343: PetscScalar wndmix; //temperature change from wind mixing: wind*Ch
347: wndmix = 0.0025 + 0.0042*wind; //regression equation valid for neutral and stable BL
348: *sheat = density*Cp*wndmix*(airtemp - sfctemp); //calculates sensible heat flux
350: return(0);
351: }
355: PetscErrorCode latentflux(PetscScalar sfctemp, PetscScalar dewtemp, PetscScalar wind, PetscScalar pressure1, PetscScalar* latentheat)
356: {
357: PetscScalar density = 1; //density of dry air
358: PetscScalar q; //actual specific humitity
359: PetscScalar qs; //saturation specific humidity
360: PetscScalar wndmix; //temperature change from wind mixing: wind*Ch
361: PetscScalar beta = .4; //moisture availability
362: PetscScalar mr ; //mixing ratio
363: PetscScalar lhcnst; //latent heat of vaporization constant = 2501000 J/kg at 0c
364: //latent heat of saturation const = 2834000 J/kg
365: //latent heat of fusion const = 333700 J/kg
368: wind = mph2mpers(wind); //converts wind from mph to meters per second
369: wndmix = 0.0025 + 0.0042*wind; //regression equation valid for neutral BL
370: lhcnst = Lconst(sfctemp); //calculates latent heat of evaporation
371: mr = calcmixingr(sfctemp,pressure1);//calculates saturation mixing ratio
372: qs = calc_q(mr); //calculates saturation specific humidty
373: mr = calcmixingr(dewtemp,pressure1); //calculates mixing ratio
374: q = calc_q(mr); //calculates specific humidty
376: *latentheat = density*wndmix*beta*lhcnst*(q - qs); //calculates latent heat flux
377: return(0);
378: }
382: PetscErrorCode potential_temperature(PetscScalar temp, PetscScalar pressure1, PetscScalar pressure2, PetscScalar sfctemp, PetscScalar* pottemp)
383: {
384: PetscScalar kdry; //poisson constant for dry atmosphere
385: PetscScalar kmoist; //poisson constant for moist atmosphere
386: PetscScalar pavg; //average atmospheric pressure
387: PetscScalar mixratio;//mixing ratio
390: mixratio = calcmixingr(sfctemp,pressure1);
392: /*initialize poisson constant */
393: kdry = 0.2854;
394: kmoist = 0.2854*(1 - 0.24*mixratio);
396: pavg = ((0.7*pressure1)+pressure2)/2; //calculates simple average press
397: *pottemp = temp*(pow((pressure1/pavg),kdry)); //calculates potential temperature
399: return(0);
400: }
401: extern PetscScalar calcmixingr(PetscScalar dtemp, PetscScalar pressure1)
402: {
403: PetscScalar e; //vapor pressure
404: PetscScalar mixratio; //mixing ratio
406: dtemp = dtemp - 273; //converts from Kelvin to Celsuis
407: e = 6.11*(pow(10,((7.5*dtemp)/(237.7+dtemp)))); //converts from dew point temp to vapor pressure
408: e = e*100; //converts from hPa to Pa
409: mixratio = (0.622*e)/(pressure1 - e); //computes mixing ratio
410: mixratio = mixratio*1; //convert to g/Kg
412: return mixratio;
413: }
414: extern PetscScalar calc_q(PetscScalar rv)
415: {
416: PetscScalar specific_humidity; //define specific humidity variable
417: specific_humidity = rv/(1 + rv); //calculates specific humidity
418: return specific_humidity;
419: }
423: PetscErrorCode calc_gflux(PetscScalar sfctemp, PetscScalar deep_grnd_temp, PetscScalar* Gflux)
424: {
425: PetscScalar k; //thermal conductivity parameter
426: PetscScalar n = 0.38; //value of soil porosity
427: PetscScalar dz = 1; //depth of layer between soil surface and deep soil layer
428: 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: }
438: extern PetscScalar emission(PetscScalar pwat)
439: {
440: PetscScalar emma;
442: emma = 0.725 + 0.17*log10(pwat);
444: return emma;
445: }
446: extern PetscScalar cloud(PetscScalar fract)
447: {
448: PetscScalar emma = 0;
450: /*modifies radiative balance depending on cloud cover */
451: if (fract >= 0.9)
452: emma = 1;
453: else if (0.9 > fract && fract >= 0.8)
454: emma = 0.9;
455: else if (0.8 > fract && fract >= 0.7)
456: emma = 0.85;
457: else if (0.7 > fract && fract >= 0.6)
458: emma = 0.75;
459: else if (0.6 > fract && fract >= 0.5)
460: emma = 0.65;
461: else if (0.4 > fract && fract >= 0.3)
462: emma = emma*1.086956;
463: return emma;
464: }
465: extern PetscScalar Lconst(PetscScalar sfctemp)
466: {
467: PetscScalar Lheat;
468: sfctemp -=273; //converts from kelvin to celsius
469: Lheat = 4186.8*(597.31 - 0.5625*sfctemp); //calculates latent heat constant
470: return Lheat;
471: }
472: extern PetscScalar mph2mpers(PetscScalar wind)
473: {
474: wind = ((wind*1.6*1000)/3600); //converts wind from mph to meters per second
475: return wind;
476: }
477: extern PetscScalar fahr_to_cel(PetscScalar temp)
478: {
479: temp = (5*(temp-32))/9; //converts from farhrenheit to celsuis
480: return temp;
481: }
482: extern PetscScalar cel_to_fahr(PetscScalar temp)
483: {
484: temp = ((temp*9)/5) + 32; //converts from celsuis to farhrenheit
485: return temp;
486: }
487: void readinput(struct in *put)
488: {
489: int i;
490: char x;
491: FILE *ifp;
493: ifp = fopen("ex5_control.txt", "r");
495: for (i=0;i<110;i++)
496: fscanf(ifp, "%c", &x);
497: fscanf(ifp, "%lf", &put->Ts);
499: for (i=0;i<43;i++)
500: fscanf(ifp, "%c", &x);
501: fscanf(ifp, "%lf", &put->Td);
503: for (i=0;i<43;i++)
504: fscanf(ifp, "%c", &x);
505: fscanf(ifp, "%lf", &put->Ta);
507: for (i=0;i<43;i++)
508: fscanf(ifp, "%c", &x);
509: fscanf(ifp, "%lf", &put->Tc);
511: for (i=0;i<43;i++)
512: fscanf(ifp, "%c", &x);
513: fscanf(ifp, "%lf", &put->fr);
515: for (i=0;i<43;i++)
516: fscanf(ifp, "%c", &x);
517: fscanf(ifp, "%lf", &put->wnd);
519: for (i=0;i<43;i++)
520: fscanf(ifp, "%c", &x);
521: fscanf(ifp, "%lf", &put->pwt);
523: for (i=0;i<43;i++)
524: fscanf(ifp, "%c", &x);
525: fscanf(ifp, "%lf", &put->wndDir);
527: for (i=0;i<43;i++)
528: fscanf(ifp, "%c", &x);
529: fscanf(ifp, "%lf", &put->time);
531: for (i=0;i<63;i++)
532: fscanf(ifp, "%c", &x);
533: fscanf(ifp, "%lf", &put->init);
534: }
536: /* ------------------------------------------------------------------- */
539: PetscErrorCode FormInitialSolution(DM da,Vec Xglobal,void *ctx)
540: {
542: AppCtx *user = (AppCtx*)ctx; /* user-defined application context */
543: PetscInt i,j,xs,ys,xm,ym,Mx,My;
544: Field **X;
545: PetscScalar deltT;
546: PetscReal hx,hy;
547: FILE *ifp;
548: FILE *ofp;
551: ofp = fopen("swing", "w");
552: ifp = fopen("grid.in", "r");
553: deltT = 0.8;
555: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
556: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
558: hx = 1/(PetscReal)(Mx-1);
559: hy = 1/(PetscReal)(My-1);
561: /* Get pointers to vector data */
562: DMDAVecGetArray(da,Xglobal,&X);
564: /* Get local grid boundaries */
565: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
567: /* Compute function over the locally owned part of the grid */
569: if (user->init == 1) {
570: for (j=ys; j<ys+ym; j++) {
571: for (i=xs; i<xs+xm; i++) {
572: X[j][i].Ts = user->Ts - i*0.0001;
573: X[j][i].Ta = X[j][i].Ts - 5;
574: X[j][i].u = 0;
575: X[j][i].v = 0;
576: X[j][i].p = 1.25;
577: if ((j == 5 || j == 6) && (i == 4 || i == 5)) X[j][i].p += 0.00001;
578: if ((j == 5 || j == 6) && (i == 12 || i == 13)) X[j][i].p += 0.00001;
579: }
580: }
581: }
583: else {
584: for (j=ys; j<ys+ym; j++) {
585: for (i=xs; i<xs+xm; i++) {
586: X[j][i].Ts = user->Ts;
587: X[j][i].Ta = X[j][i].Ts - 5;
588: X[j][i].u = 0;
589: X[j][i].v = 0;
590: X[j][i].p = 1.25;
591: }
592: }
593: }
595: /* Restore vectors */
596: DMDAVecRestoreArray(da,Xglobal,&X);
597: return(0);
598: }
602: /*
603: RhsFunc - Evaluates nonlinear function F(u).
605: Input Parameters:
606: . ts - the TS context
607: . t - current time
608: . Xglobal - input vector
609: . F - output vector
610: . ptr - optional user-defined context, as set by SNESSetFunction()
612: Output Parameter:
613: . F - rhs function vector
614: */
615: PetscErrorCode RhsFunc(TS ts,PetscReal t,Vec Xglobal,Vec F,void *ctx)
616: {
617: AppCtx *user = (AppCtx*)ctx; /* user-defined application context */
618: DM da = user->da;
620: PetscInt i,j,Mx,My,xs,ys,xm,ym;
621: PetscReal dhx,dhy;
622: Vec localT;
623: Field **X,**Frhs; //structures that contain variables of interest and left hand side of governing equations respectively
624: PetscScalar csoil = user->csoil; //heat constant for layer
625: PetscScalar dzlay = user->dzlay; //thickness of top soil layer
626: PetscScalar emma = user->emma; //emission parameter
627: PetscScalar wind = user->wind; //wind speed
628: PetscScalar dewtemp = user->dewtemp; //dew point temperature (moisture in air)
629: PetscScalar pressure1 = user->pressure1; //sea level pressure
630: PetscScalar airtemp = user->airtemp; //temperature of air near boundary layer inversion
631: PetscScalar fract = user->fract; //fraction of the sky covered by clouds
632: PetscScalar Tc = user->Tc; //temperature at base of lowest cloud layer
633: PetscScalar lat = user->lat; //latitude
634: PetscScalar Cp = 1005.7; //specific heat of air at constant pressure
635: PetscScalar Rd = 287.058; //gas constant for dry air
636: PetscScalar diffconst = 1000; //diffusion coefficient
637: PetscScalar f = 2*0.0000727*sin(lat); //coriolis force
638: PetscScalar deep_grnd_temp = user->deep_grnd_temp; //temp in lowest ground layer
639: PetscScalar Ts,u,v,p,P;
640: PetscScalar u_abs,u_plus,u_minus,v_abs,v_plus,v_minus;
642: PetscScalar sfctemp1,fsfc1,Ra;
643: PetscScalar sheat; //sensible heat flux
644: PetscScalar latentheat; //latent heat flux
645: PetscScalar groundflux; //flux from conduction of deep ground layer in contact with top soil
646: PetscInt xend,yend;
649: DMGetLocalVector(da,&localT);
650: DMDAGetInfo(da,PETSC_IGNORE,&Mx,&My,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,
651: PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE,PETSC_IGNORE);
653: dhx = (PetscReal)(Mx-1)/(5000*(Mx-1)); // dhx = 1/dx; assume 2D space domain: [0.0, 1.e5] x [0.0, 1.e5]
654: dhy = (PetscReal)(My-1)/(5000*(Mx-1)); // dhy = 1/dy;
657: /*
658: Scatter ghost points to local vector,using the 2-step process
659: DAGlobalToLocalBegin(),DAGlobalToLocalEnd().
660: By placing code between these two statements, computations can be
661: done while messages are in transition.
662: */
663: DMGlobalToLocalBegin(da,Xglobal,INSERT_VALUES,localT);
664: DMGlobalToLocalEnd(da,Xglobal,INSERT_VALUES,localT);
666: /* Get pointers to vector data */
667: DMDAVecGetArray(da,localT,&X);
668: DMDAVecGetArray(da,F,&Frhs);
670: /* Get local grid boundaries */
671: DMDAGetCorners(da,&xs,&ys,PETSC_NULL,&xm,&ym,PETSC_NULL);
673: /* Compute function over the locally owned part of the grid */
674: /* the interior points */
675: xend=xs+xm; yend=ys+ym;
676: for (j=ys; j<yend; j++) {
677: for (i=xs; i<xend; i++) {
678: 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;
680: sfctemp1 = (double)Ts;
681: sfctemp1 = (double)X[j][i].Ts;
682: calcfluxs(sfctemp1,airtemp,emma,fract,Tc,&fsfc1); //calculates surface net radiative flux
683: sensibleflux(sfctemp1,airtemp,wind,&sheat); //calculate sensible heat flux
684: latentflux(sfctemp1,dewtemp,wind,pressure1,&latentheat);//calculates latent heat flux
685: calc_gflux(sfctemp1,deep_grnd_temp,&groundflux); //calculates flux from earth below surface soil layer by conduction
686: calcfluxa(sfctemp1,airtemp,emma,&Ra); //Calculates the change in downward radiative flux
687: fsfc1 = fsfc1 + latentheat + sheat + groundflux; //adds radiative, sensible heat, latent heat, and ground heat flux yielding net flux
689: /* convective coefficients for upwinding */
690: u_abs = PetscAbsScalar(u);
691: u_plus = .5*(u + u_abs); // u if u>0; 0 if u<0
692: u_minus = .5*(u - u_abs); // u if u <0; 0 if u>0
694: v_abs = PetscAbsScalar(v);
695: v_plus = .5*(v + v_abs); // v if v>0; 0 if v<0
696: v_minus = .5*(v - v_abs); // v if v <0; 0 if v>0
698: /* Solve governing equations */
699: P = p*Rd*Ts;
701: /* du/dt -> time change of east-west component of the wind */
702: 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)
703: - v_plus*(u - X[j-1][i].u)*dhy - v_minus*(X[j+1][i].u - u)*dhy // - v(du/dy)
704: -(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)]
705: // -(1/p)*(X[j][i+1].P - X[j][i-1].P)*dhx
706: + f*v;
708: /* dv/dt -> time change of north-south component of the wind */
709: 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)
710: - v_plus*(v - X[j-1][i].v)*dhy - v_minus*(X[j+1][i].v - v)*dhy // - v(dv/dy)
711: -(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)]
712: // -(1/p)*(X[j+1][i].P - X[j-1][i].P)*dhy
713: -f*u;
715: /* dT/dt -> time change of temperature */
716: Frhs[j][i].Ts = (fsfc1/(csoil*dzlay)) // Fnet/(Cp*dz) diabatic change in T
717: -u_plus*(Ts - X[j][i-1].Ts)*dhx - u_minus*(X[j][i+1].Ts - Ts)*dhx // - u*(dTs/dx) advection x
718: -v_plus*(Ts - X[j-1][i].Ts)*dhy - v_minus*(X[j+1][i].Ts - Ts)*dhy // - v*(dTs/dy) advection y
719: + diffconst*((X[j][i+1].Ts - 2*Ts + X[j][i-1].Ts)*dhx*dhx // + D(Ts_xx + Ts_yy) diffusion
720: + (X[j+1][i].Ts - 2*Ts + X[j-1][i].Ts)*dhy*dhy);
722: /* dp/dt -> time change of */
723: 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)
724: -v_plus*(p - X[j-1][i].p)*dhy - v_minus*(X[j+1][i].p - p)*dhy; // - v*(dp/dy)
726: Frhs[j][i].Ta = Ra/Cp; //dTa/dt time change of air temperature
727: }
728: }
730: /* Restore vectors */
731: DMDAVecRestoreArray(da,localT,&X);
732: DMDAVecRestoreArray(da,F,&Frhs);
733: DMRestoreLocalVector(da,&localT);
734: return(0);
735: }
739: PetscErrorCode Monitor(TS ts,PetscInt step,PetscReal time,Vec T,void *ctx)
740: {
742: PetscScalar *array;
743: PetscInt itime=(PetscInt)(time);
744: const TSType type;
745: PetscBool sundials;
746: MonitorCtx *user = (MonitorCtx*)ctx;
747: PetscViewer viewer = user->drawviewer;
748: PetscMPIInt rank;
749: PetscReal norm;
752: MPI_Comm_rank(((PetscObject)ts)->comm,&rank);
754: TSGetType(ts,&type);
755: PetscObjectTypeCompare((PetscObject)ts,TSSUNDIALS,&sundials);
756: VecNorm(T,NORM_INFINITY,&norm);
757: if (sundials || itime%60 == 0){
758: VecGetArray(T,&array);
759: if (!rank){printf("step %4d, time %8.1f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f\n",step,time,(((array[0]-273)*9)/5 + 32),(((array[1]-273)*9)/5 + 32),array[2],array[3],array[4],array[5]);}
760: VecRestoreArray(T,&array);
762: if (user->drawcontours){
763: VecView(T,viewer);
764: }
765: }
766: return(0);
767: }