1c4762a1bSJed Brown static char help[] = "Nonlinear, time-dependent. Developed from radiative_surface_balance.c \n"; 2c4762a1bSJed Brown /* 3c4762a1bSJed Brown Contributed by Steve Froehlich, Illinois Institute of Technology 4c4762a1bSJed Brown 5c4762a1bSJed Brown Usage: 6c4762a1bSJed Brown mpiexec -n <np> ./ex5 [options] 7c4762a1bSJed Brown ./ex5 -help [view petsc options] 8c4762a1bSJed Brown ./ex5 -ts_type sundials -ts_view 9c4762a1bSJed Brown ./ex5 -da_grid_x 20 -da_grid_y 20 -log_view 10c4762a1bSJed Brown ./ex5 -da_grid_x 20 -da_grid_y 20 -ts_type rosw -ts_atol 1.e-6 -ts_rtol 1.e-6 11c4762a1bSJed Brown ./ex5 -drawcontours -draw_pause 0.1 -draw_fields 0,1,2,3,4 12c4762a1bSJed Brown */ 13c4762a1bSJed Brown 14c4762a1bSJed Brown /* 15c4762a1bSJed Brown ----------------------------------------------------------------------- 16c4762a1bSJed Brown 17c4762a1bSJed Brown Governing equations: 18c4762a1bSJed Brown 19c4762a1bSJed Brown R = s*(Ea*Ta^4 - Es*Ts^4) 20c4762a1bSJed Brown SH = p*Cp*Ch*wind*(Ta - Ts) 21c4762a1bSJed Brown LH = p*L*Ch*wind*B(q(Ta) - q(Ts)) 22c4762a1bSJed Brown G = k*(Tgnd - Ts)/dz 23c4762a1bSJed Brown 24c4762a1bSJed Brown Fnet = R + SH + LH + G 25c4762a1bSJed Brown 26c4762a1bSJed Brown du/dt = -u*(du/dx) - v*(du/dy) - 2*omeg*sin(lat)*v - (1/p)*(dP/dx) 27c4762a1bSJed Brown dv/dt = -u*(dv/dx) - v*(dv/dy) + 2*omeg*sin(lat)*u - (1/p)*(dP/dy) 28c4762a1bSJed Brown dTs/dt = Fnet/(Cp*dz) - Div([u*Ts, v*Ts]) + D*Lap(Ts) 29c4762a1bSJed Brown = Fnet/(Cs*dz) - u*(dTs/dx) - v*(dTs/dy) + D*(Ts_xx + Ts_yy) 30c4762a1bSJed Brown dp/dt = -Div([u*p,v*p]) 31c4762a1bSJed Brown = - u*dp/dx - v*dp/dy 32c4762a1bSJed Brown dTa/dt = Fnet/Cp 33c4762a1bSJed Brown 34c4762a1bSJed Brown Equation of State: 35c4762a1bSJed Brown 36c4762a1bSJed Brown P = p*R*Ts 37c4762a1bSJed Brown 38c4762a1bSJed Brown ----------------------------------------------------------------------- 39c4762a1bSJed Brown 40c4762a1bSJed Brown Program considers the evolution of a two dimensional atmosphere from 41c4762a1bSJed Brown sunset to sunrise. There are two components: 42c4762a1bSJed Brown 1. Surface energy balance model to compute diabatic dT (Fnet) 43c4762a1bSJed Brown 2. Dynamical model using simplified primitive equations 44c4762a1bSJed Brown 45c4762a1bSJed Brown Program is to be initiated at sunset and run to sunrise. 46c4762a1bSJed Brown 47c4762a1bSJed Brown Inputs are: 48c4762a1bSJed Brown Surface temperature 49c4762a1bSJed Brown Dew point temperature 50c4762a1bSJed Brown Air temperature 51c4762a1bSJed Brown Temperature at cloud base (if clouds are present) 52c4762a1bSJed Brown Fraction of sky covered by clouds 53c4762a1bSJed Brown Wind speed 54c4762a1bSJed Brown Precipitable water in centimeters 55c4762a1bSJed Brown Wind direction 56c4762a1bSJed Brown 57c4762a1bSJed Brown Inputs are are read in from the text file ex5_control.txt. To change an 58c4762a1bSJed Brown input value use ex5_control.txt. 59c4762a1bSJed Brown 60c4762a1bSJed Brown Solvers: 61c4762a1bSJed Brown Backward Euler = default solver 62c4762a1bSJed Brown Sundials = fastest and most accurate, requires Sundials libraries 63c4762a1bSJed Brown 64c4762a1bSJed Brown This model is under development and should be used only as an example 65c4762a1bSJed Brown and not as a predictive weather model. 66c4762a1bSJed Brown */ 67c4762a1bSJed Brown 68c4762a1bSJed Brown #include <petscts.h> 69c4762a1bSJed Brown #include <petscdm.h> 70c4762a1bSJed Brown #include <petscdmda.h> 71c4762a1bSJed Brown 72c4762a1bSJed Brown /* stefan-boltzmann constant */ 73c4762a1bSJed Brown #define SIG 0.000000056703 74c4762a1bSJed Brown /* absorption-emission constant for surface */ 75c4762a1bSJed Brown #define EMMSFC 1 76c4762a1bSJed Brown /* amount of time (seconds) that passes before new flux is calculated */ 77c4762a1bSJed Brown #define TIMESTEP 1 78c4762a1bSJed Brown 79c4762a1bSJed Brown /* variables of interest to be solved at each grid point */ 80c4762a1bSJed Brown typedef struct { 81c4762a1bSJed Brown PetscScalar Ts, Ta; /* surface and air temperature */ 82c4762a1bSJed Brown PetscScalar u, v; /* wind speed */ 83c4762a1bSJed Brown PetscScalar p; /* density */ 84c4762a1bSJed Brown } Field; 85c4762a1bSJed Brown 86c4762a1bSJed Brown /* User defined variables. Used in solving for variables of interest */ 87c4762a1bSJed Brown typedef struct { 88c4762a1bSJed Brown DM da; /* grid */ 89c4762a1bSJed Brown PetscScalar csoil; /* heat constant for layer */ 90c4762a1bSJed Brown PetscScalar dzlay; /* thickness of top soil layer */ 91c4762a1bSJed Brown PetscScalar emma; /* emission parameter */ 92c4762a1bSJed Brown PetscScalar wind; /* wind speed */ 93c4762a1bSJed Brown PetscScalar dewtemp; /* dew point temperature (moisture in air) */ 94c4762a1bSJed Brown PetscScalar pressure1; /* sea level pressure */ 95c4762a1bSJed Brown PetscScalar airtemp; /* temperature of air near boundary layer inversion */ 96c4762a1bSJed Brown PetscScalar Ts; /* temperature at the surface */ 97c4762a1bSJed Brown PetscScalar fract; /* fraction of sky covered by clouds */ 98c4762a1bSJed Brown PetscScalar Tc; /* temperature at base of lowest cloud layer */ 99c4762a1bSJed Brown PetscScalar lat; /* Latitude in degrees */ 100c4762a1bSJed Brown PetscScalar init; /* initialization scenario */ 101c4762a1bSJed Brown PetscScalar deep_grnd_temp; /* temperature of ground under top soil surface layer */ 102c4762a1bSJed Brown } AppCtx; 103c4762a1bSJed Brown 104c4762a1bSJed Brown /* Struct for visualization */ 105c4762a1bSJed Brown typedef struct { 106c4762a1bSJed Brown PetscBool drawcontours; /* flag - 1 indicates drawing contours */ 107c4762a1bSJed Brown PetscViewer drawviewer; 108c4762a1bSJed Brown PetscInt interval; 109c4762a1bSJed Brown } MonitorCtx; 110c4762a1bSJed Brown 111c4762a1bSJed Brown /* Inputs read in from text file */ 112c4762a1bSJed Brown struct in { 113c4762a1bSJed Brown PetscScalar Ts; /* surface temperature */ 114c4762a1bSJed Brown PetscScalar Td; /* dewpoint temperature */ 115c4762a1bSJed Brown PetscScalar Tc; /* temperature of cloud base */ 116c4762a1bSJed Brown PetscScalar fr; /* fraction of sky covered by clouds */ 117c4762a1bSJed Brown PetscScalar wnd; /* wind speed */ 118c4762a1bSJed Brown PetscScalar Ta; /* air temperature */ 119c4762a1bSJed Brown PetscScalar pwt; /* precipitable water */ 120c4762a1bSJed Brown PetscScalar wndDir; /* wind direction */ 121c4762a1bSJed Brown PetscScalar lat; /* latitude */ 122c4762a1bSJed Brown PetscReal time; /* time in hours */ 123c4762a1bSJed Brown PetscScalar init; 124c4762a1bSJed Brown }; 125c4762a1bSJed Brown 126c4762a1bSJed Brown /* functions */ 127c4762a1bSJed Brown extern PetscScalar emission(PetscScalar); /* sets emission/absorption constant depending on water vapor content */ 128c4762a1bSJed Brown extern PetscScalar calc_q(PetscScalar); /* calculates specific humidity */ 129c4762a1bSJed Brown extern PetscScalar mph2mpers(PetscScalar); /* converts miles per hour to meters per second */ 130c4762a1bSJed Brown 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. */ 131c4762a1bSJed Brown extern PetscScalar fahr_to_cel(PetscScalar); /* converts Fahrenheit to Celsius */ 132c4762a1bSJed Brown extern PetscScalar cel_to_fahr(PetscScalar); /* converts Celsius to Fahrenheit */ 133c4762a1bSJed Brown extern PetscScalar calcmixingr(PetscScalar, PetscScalar); /* calculates mixing ratio */ 134c4762a1bSJed Brown extern PetscScalar cloud(PetscScalar); /* cloud radiative parameterization */ 135c4762a1bSJed Brown extern PetscErrorCode FormInitialSolution(DM, Vec, void *); /* Specifies initial conditions for the system of equations (PETSc defined function) */ 136c4762a1bSJed Brown extern PetscErrorCode RhsFunc(TS, PetscReal, Vec, Vec, void *); /* Specifies the user defined functions (PETSc defined function) */ 137c4762a1bSJed Brown extern PetscErrorCode Monitor(TS, PetscInt, PetscReal, Vec, void *); /* Specifies output and visualization tools (PETSc defined function) */ 138303a5415SBarry Smith extern PetscErrorCode readinput(struct in *put); /* reads input from text file */ 139c4762a1bSJed Brown extern PetscErrorCode calcfluxs(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar *); /* calculates upward IR from surface */ 140c4762a1bSJed Brown extern PetscErrorCode calcfluxa(PetscScalar, PetscScalar, PetscScalar, PetscScalar *); /* calculates downward IR from atmosphere */ 141c4762a1bSJed Brown extern PetscErrorCode sensibleflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar *); /* calculates sensible heat flux */ 142c4762a1bSJed Brown extern PetscErrorCode potential_temperature(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar *); /* calculates potential temperature */ 143c4762a1bSJed Brown extern PetscErrorCode latentflux(PetscScalar, PetscScalar, PetscScalar, PetscScalar, PetscScalar *); /* calculates latent heat flux */ 144c4762a1bSJed Brown extern PetscErrorCode calc_gflux(PetscScalar, PetscScalar, PetscScalar *); /* calculates flux between top soil layer and underlying earth */ 145c4762a1bSJed Brown 146d71ae5a4SJacob Faibussowitsch int main(int argc, char **argv) 147d71ae5a4SJacob Faibussowitsch { 148303a5415SBarry Smith PetscInt time; /* amount of loops */ 149c4762a1bSJed Brown struct in put; 150c4762a1bSJed Brown PetscScalar rh; /* relative humidity */ 151c4762a1bSJed Brown PetscScalar x; /* memory varialbe for relative humidity calculation */ 152c4762a1bSJed Brown PetscScalar deep_grnd_temp; /* temperature of ground under top soil surface layer */ 153c4762a1bSJed Brown PetscScalar emma; /* absorption-emission constant for air */ 154c4762a1bSJed Brown PetscScalar pressure1 = 101300; /* surface pressure */ 155c4762a1bSJed Brown PetscScalar mixratio; /* mixing ratio */ 156c4762a1bSJed Brown PetscScalar airtemp; /* temperature of air near boundary layer inversion */ 157c4762a1bSJed Brown PetscScalar dewtemp; /* dew point temperature */ 158c4762a1bSJed Brown PetscScalar sfctemp; /* temperature at surface */ 159c4762a1bSJed Brown PetscScalar pwat; /* total column precipitable water */ 160c4762a1bSJed Brown PetscScalar cloudTemp; /* temperature at base of cloud */ 161c4762a1bSJed Brown AppCtx user; /* user-defined work context */ 162c4762a1bSJed Brown MonitorCtx usermonitor; /* user-defined monitor context */ 163c4762a1bSJed Brown TS ts; 164c4762a1bSJed Brown SNES snes; 165c4762a1bSJed Brown DM da; 166c4762a1bSJed Brown Vec T, rhs; /* solution vector */ 167c4762a1bSJed Brown Mat J; /* Jacobian matrix */ 168c4762a1bSJed Brown PetscReal ftime, dt; 169c4762a1bSJed Brown PetscInt steps, dof = 5; 170c4762a1bSJed Brown PetscBool use_coloring = PETSC_TRUE; 171c4762a1bSJed Brown MatFDColoring matfdcoloring = 0; 172c4762a1bSJed Brown PetscBool monitor_off = PETSC_FALSE; 173*252b6434SHong Zhang PetscBool removezero = PETSC_FALSE; 174c4762a1bSJed Brown 175327415f7SBarry Smith PetscFunctionBeginUser; 1769566063dSJacob Faibussowitsch PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); 177c4762a1bSJed Brown 178c4762a1bSJed Brown /* Inputs */ 1799566063dSJacob Faibussowitsch PetscCall(readinput(&put)); 180c4762a1bSJed Brown 181c4762a1bSJed Brown sfctemp = put.Ts; 182c4762a1bSJed Brown dewtemp = put.Td; 183c4762a1bSJed Brown cloudTemp = put.Tc; 184c4762a1bSJed Brown airtemp = put.Ta; 185c4762a1bSJed Brown pwat = put.pwt; 186c4762a1bSJed Brown 1879566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial Temperature = %g\n", (double)sfctemp)); /* input surface temperature */ 188c4762a1bSJed Brown 189c4762a1bSJed Brown deep_grnd_temp = sfctemp - 10; /* set underlying ground layer temperature */ 190c4762a1bSJed Brown emma = emission(pwat); /* accounts for radiative effects of water vapor */ 191c4762a1bSJed Brown 192c4762a1bSJed Brown /* Converts from Fahrenheit to Celsuis */ 193c4762a1bSJed Brown sfctemp = fahr_to_cel(sfctemp); 194c4762a1bSJed Brown airtemp = fahr_to_cel(airtemp); 195c4762a1bSJed Brown dewtemp = fahr_to_cel(dewtemp); 196c4762a1bSJed Brown cloudTemp = fahr_to_cel(cloudTemp); 197c4762a1bSJed Brown deep_grnd_temp = fahr_to_cel(deep_grnd_temp); 198c4762a1bSJed Brown 199c4762a1bSJed Brown /* Converts from Celsius to Kelvin */ 200c4762a1bSJed Brown sfctemp += 273; 201c4762a1bSJed Brown airtemp += 273; 202c4762a1bSJed Brown dewtemp += 273; 203c4762a1bSJed Brown cloudTemp += 273; 204c4762a1bSJed Brown deep_grnd_temp += 273; 205c4762a1bSJed Brown 206c4762a1bSJed Brown /* Calculates initial relative humidity */ 207c4762a1bSJed Brown x = calcmixingr(dewtemp, pressure1); 208c4762a1bSJed Brown mixratio = calcmixingr(sfctemp, pressure1); 209c4762a1bSJed Brown rh = (x / mixratio) * 100; 210c4762a1bSJed Brown 2119566063dSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Initial RH = %.1f percent\n\n", (double)rh)); /* prints initial relative humidity */ 212c4762a1bSJed Brown 213c4762a1bSJed Brown time = 3600 * put.time; /* sets amount of timesteps to run model */ 214c4762a1bSJed Brown 215c4762a1bSJed Brown /* Configure PETSc TS solver */ 216c4762a1bSJed Brown /*------------------------------------------*/ 217c4762a1bSJed Brown 218c4762a1bSJed Brown /* Create grid */ 2199566063dSJacob Faibussowitsch PetscCall(DMDACreate2d(PETSC_COMM_WORLD, DM_BOUNDARY_PERIODIC, DM_BOUNDARY_PERIODIC, DMDA_STENCIL_STAR, 20, 20, PETSC_DECIDE, PETSC_DECIDE, dof, 1, NULL, NULL, &da)); 2209566063dSJacob Faibussowitsch PetscCall(DMSetFromOptions(da)); 2219566063dSJacob Faibussowitsch PetscCall(DMSetUp(da)); 2229566063dSJacob Faibussowitsch PetscCall(DMDASetUniformCoordinates(da, 0.0, 1.0, 0.0, 1.0, 0.0, 1.0)); 223c4762a1bSJed Brown 224c4762a1bSJed Brown /* Define output window for each variable of interest */ 2259566063dSJacob Faibussowitsch PetscCall(DMDASetFieldName(da, 0, "Ts")); 2269566063dSJacob Faibussowitsch PetscCall(DMDASetFieldName(da, 1, "Ta")); 2279566063dSJacob Faibussowitsch PetscCall(DMDASetFieldName(da, 2, "u")); 2289566063dSJacob Faibussowitsch PetscCall(DMDASetFieldName(da, 3, "v")); 2299566063dSJacob Faibussowitsch PetscCall(DMDASetFieldName(da, 4, "p")); 230c4762a1bSJed Brown 231c4762a1bSJed Brown /* set values for appctx */ 232c4762a1bSJed Brown user.da = da; 233c4762a1bSJed Brown user.Ts = sfctemp; 234c4762a1bSJed Brown user.fract = put.fr; /* fraction of sky covered by clouds */ 235c4762a1bSJed Brown user.dewtemp = dewtemp; /* dew point temperature (mositure in air) */ 236c4762a1bSJed Brown user.csoil = 2000000; /* heat constant for layer */ 237c4762a1bSJed Brown user.dzlay = 0.08; /* thickness of top soil layer */ 238c4762a1bSJed Brown user.emma = emma; /* emission parameter */ 239c4762a1bSJed Brown user.wind = put.wnd; /* wind spped */ 240c4762a1bSJed Brown user.pressure1 = pressure1; /* sea level pressure */ 241c4762a1bSJed Brown user.airtemp = airtemp; /* temperature of air near boundar layer inversion */ 242c4762a1bSJed Brown user.Tc = cloudTemp; /* temperature at base of lowest cloud layer */ 243c4762a1bSJed Brown user.init = put.init; /* user chosen initiation scenario */ 244c4762a1bSJed Brown user.lat = 70 * 0.0174532; /* converts latitude degrees to latitude in radians */ 245c4762a1bSJed Brown user.deep_grnd_temp = deep_grnd_temp; /* temp in lowest ground layer */ 246c4762a1bSJed Brown 247c4762a1bSJed Brown /* set values for MonitorCtx */ 248c4762a1bSJed Brown usermonitor.drawcontours = PETSC_FALSE; 2499566063dSJacob Faibussowitsch PetscCall(PetscOptionsHasName(NULL, NULL, "-drawcontours", &usermonitor.drawcontours)); 250c4762a1bSJed Brown if (usermonitor.drawcontours) { 251c4762a1bSJed Brown PetscReal bounds[] = {1000.0, -1000., -1000., -1000., 1000., -1000., 1000., -1000., 1000, -1000, 100700, 100800}; 2529566063dSJacob Faibussowitsch PetscCall(PetscViewerDrawOpen(PETSC_COMM_WORLD, 0, 0, 0, 0, 300, 300, &usermonitor.drawviewer)); 2539566063dSJacob Faibussowitsch PetscCall(PetscViewerDrawSetBounds(usermonitor.drawviewer, dof, bounds)); 254c4762a1bSJed Brown } 255c4762a1bSJed Brown usermonitor.interval = 1; 2569566063dSJacob Faibussowitsch PetscCall(PetscOptionsGetInt(NULL, NULL, "-monitor_interval", &usermonitor.interval, NULL)); 257c4762a1bSJed Brown 258c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 259c4762a1bSJed Brown Extract global vectors from DA; 260c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 2619566063dSJacob Faibussowitsch PetscCall(DMCreateGlobalVector(da, &T)); 2629566063dSJacob Faibussowitsch PetscCall(VecDuplicate(T, &rhs)); /* r: vector to put the computed right hand side */ 263c4762a1bSJed Brown 2649566063dSJacob Faibussowitsch PetscCall(TSCreate(PETSC_COMM_WORLD, &ts)); 2659566063dSJacob Faibussowitsch PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); 2669566063dSJacob Faibussowitsch PetscCall(TSSetType(ts, TSBEULER)); 2679566063dSJacob Faibussowitsch PetscCall(TSSetRHSFunction(ts, rhs, RhsFunc, &user)); 268c4762a1bSJed Brown 269c4762a1bSJed Brown /* Set Jacobian evaluation routine - use coloring to compute finite difference Jacobian efficiently */ 2709566063dSJacob Faibussowitsch PetscCall(DMSetMatType(da, MATAIJ)); 2719566063dSJacob Faibussowitsch PetscCall(DMCreateMatrix(da, &J)); 272c4762a1bSJed Brown if (use_coloring) { 273c4762a1bSJed Brown ISColoring iscoloring; 274*252b6434SHong Zhang PetscInt ncolors; 275*252b6434SHong Zhang 2769566063dSJacob Faibussowitsch PetscCall(DMCreateColoring(da, IS_COLORING_GLOBAL, &iscoloring)); 2779566063dSJacob Faibussowitsch PetscCall(MatFDColoringCreate(J, iscoloring, &matfdcoloring)); 2789566063dSJacob Faibussowitsch PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); 2799566063dSJacob Faibussowitsch PetscCall(MatFDColoringSetUp(J, iscoloring, matfdcoloring)); 280*252b6434SHong Zhang PetscCall(ISColoringGetColors(iscoloring, NULL, &ncolors, NULL)); 281*252b6434SHong Zhang PetscCall(PetscPrintf(PETSC_COMM_WORLD, "DMColoring generates %" PetscInt_FMT " colors\n", ncolors)); 2829566063dSJacob Faibussowitsch PetscCall(ISColoringDestroy(&iscoloring)); 283*252b6434SHong Zhang PetscCall(TSSetIJacobian(ts, J, J, TSComputeIJacobianDefaultColor, NULL)); 284c4762a1bSJed Brown } else { 285*252b6434SHong Zhang PetscCall(TSGetSNES(ts, &snes)); 2869566063dSJacob Faibussowitsch PetscCall(SNESSetJacobian(snes, J, J, SNESComputeJacobianDefault, NULL)); 287c4762a1bSJed Brown } 288c4762a1bSJed Brown 289c4762a1bSJed Brown /* Define what to print for ts_monitor option */ 2909566063dSJacob Faibussowitsch PetscCall(PetscOptionsHasName(NULL, NULL, "-monitor_off", &monitor_off)); 29148a46eb9SPierre Jolivet if (!monitor_off) PetscCall(TSMonitorSet(ts, Monitor, &usermonitor, NULL)); 292c4762a1bSJed Brown dt = TIMESTEP; /* initial time step */ 293c4762a1bSJed Brown ftime = TIMESTEP * time; 29463a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "time %" PetscInt_FMT ", ftime %g hour, TIMESTEP %g\n", time, (double)(ftime / 3600), (double)dt)); 295c4762a1bSJed Brown 2969566063dSJacob Faibussowitsch PetscCall(TSSetTimeStep(ts, dt)); 2979566063dSJacob Faibussowitsch PetscCall(TSSetMaxSteps(ts, time)); 2989566063dSJacob Faibussowitsch PetscCall(TSSetMaxTime(ts, ftime)); 2999566063dSJacob Faibussowitsch PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); 3009566063dSJacob Faibussowitsch PetscCall(TSSetSolution(ts, T)); 3019566063dSJacob Faibussowitsch PetscCall(TSSetDM(ts, da)); 302c4762a1bSJed Brown 303c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 304c4762a1bSJed Brown Set runtime options 305c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 306*252b6434SHong Zhang PetscCall(PetscOptionsGetBool(NULL, NULL, "-removezero", &removezero, NULL)); 3079566063dSJacob Faibussowitsch PetscCall(TSSetFromOptions(ts)); 308*252b6434SHong Zhang if (removezero && matfdcoloring) { 309*252b6434SHong Zhang PetscRandom rctx; 310*252b6434SHong Zhang Vec Tdot; 311*252b6434SHong Zhang MatColoring mc; 312*252b6434SHong Zhang ISColoring iscoloring; 313*252b6434SHong Zhang PetscInt ncolors; 314*252b6434SHong Zhang /* Compute the Jacobian with randomized vector values to capture the sparsity pattern for coloring */ 315*252b6434SHong Zhang PetscCall(VecDuplicate(T, &Tdot)); 316*252b6434SHong Zhang PetscCall(PetscRandomCreate(PETSC_COMM_WORLD, &rctx)); 317*252b6434SHong Zhang PetscCall(PetscRandomSetInterval(rctx, 1.0, 2.0)); 318*252b6434SHong Zhang PetscCall(VecSetRandom(T, rctx)); 319*252b6434SHong Zhang PetscCall(VecSetRandom(Tdot, rctx)); 320*252b6434SHong Zhang PetscCall(PetscRandomDestroy(&rctx)); 321*252b6434SHong Zhang PetscCall(TSSetUp(ts)); 322*252b6434SHong Zhang PetscCall(TSComputeIJacobian(ts, 0.0, T, Tdot, 0.12345, J, J, PETSC_FALSE)); 323*252b6434SHong Zhang PetscCall(VecDestroy(&Tdot)); 324*252b6434SHong Zhang /* Generate new coloring after eliminating zeros in the matrix */ 325*252b6434SHong Zhang PetscCall(MatEliminateZeros(J)); 326*252b6434SHong Zhang PetscCall(MatColoringCreate(J, &mc)); 327*252b6434SHong Zhang PetscCall(MatColoringSetDistance(mc, 2)); 328*252b6434SHong Zhang PetscCall(MatColoringSetType(mc, MATCOLORINGSL)); 329*252b6434SHong Zhang PetscCall(MatColoringSetFromOptions(mc)); 330*252b6434SHong Zhang PetscCall(MatColoringApply(mc, &iscoloring)); 331*252b6434SHong Zhang PetscCall(MatColoringDestroy(&mc)); 332*252b6434SHong Zhang PetscCall(ISColoringGetColors(iscoloring, NULL, &ncolors, NULL)); 333*252b6434SHong Zhang PetscCall(PetscPrintf(PETSC_COMM_WORLD, "MatColoring after eliminating zeros generates %" PetscInt_FMT " colors\n", ncolors)); 334*252b6434SHong Zhang /* Replace the old coloring with the new one */ 335*252b6434SHong Zhang PetscCall(MatFDColoringDestroy(&matfdcoloring)); 336*252b6434SHong Zhang PetscCall(MatFDColoringCreate(J, iscoloring, &matfdcoloring)); 337*252b6434SHong Zhang PetscCall(MatFDColoringSetFunction(matfdcoloring, (PetscErrorCode(*)(void))SNESTSFormFunction, (void *)ts)); 338*252b6434SHong Zhang PetscCall(MatFDColoringSetFromOptions(matfdcoloring)); 339*252b6434SHong Zhang PetscCall(MatFDColoringSetUp(J, iscoloring, matfdcoloring)); 340*252b6434SHong Zhang PetscCall(PetscObjectCompose((PetscObject)J, "TSMatFDColoring", (PetscObject)matfdcoloring)); 341*252b6434SHong Zhang PetscCall(ISColoringDestroy(&iscoloring)); 342*252b6434SHong Zhang } 343c4762a1bSJed Brown /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 344c4762a1bSJed Brown Solve nonlinear system 345c4762a1bSJed Brown - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ 346*252b6434SHong Zhang PetscCall(FormInitialSolution(da, T, &user)); 3479566063dSJacob Faibussowitsch PetscCall(TSSolve(ts, T)); 3489566063dSJacob Faibussowitsch PetscCall(TSGetSolveTime(ts, &ftime)); 3499566063dSJacob Faibussowitsch PetscCall(TSGetStepNumber(ts, &steps)); 35063a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solution T after %g hours %" PetscInt_FMT " steps\n", (double)(ftime / 3600), steps)); 351c4762a1bSJed Brown 3529566063dSJacob Faibussowitsch if (matfdcoloring) PetscCall(MatFDColoringDestroy(&matfdcoloring)); 35348a46eb9SPierre Jolivet if (usermonitor.drawcontours) PetscCall(PetscViewerDestroy(&usermonitor.drawviewer)); 3549566063dSJacob Faibussowitsch PetscCall(MatDestroy(&J)); 3559566063dSJacob Faibussowitsch PetscCall(VecDestroy(&T)); 3569566063dSJacob Faibussowitsch PetscCall(VecDestroy(&rhs)); 3579566063dSJacob Faibussowitsch PetscCall(TSDestroy(&ts)); 3589566063dSJacob Faibussowitsch PetscCall(DMDestroy(&da)); 359c4762a1bSJed Brown 3609566063dSJacob Faibussowitsch PetscCall(PetscFinalize()); 361b122ec5aSJacob Faibussowitsch return 0; 362c4762a1bSJed Brown } 363c4762a1bSJed Brown /*****************************end main program********************************/ 364c4762a1bSJed Brown /*****************************************************************************/ 365c4762a1bSJed Brown /*****************************************************************************/ 366c4762a1bSJed Brown /*****************************************************************************/ 367d71ae5a4SJacob Faibussowitsch PetscErrorCode calcfluxs(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar fract, PetscScalar cloudTemp, PetscScalar *flux) 368d71ae5a4SJacob Faibussowitsch { 369c4762a1bSJed Brown PetscFunctionBeginUser; 370c4762a1bSJed Brown *flux = SIG * ((EMMSFC * emma * PetscPowScalarInt(airtemp, 4)) + (EMMSFC * fract * (1 - emma) * PetscPowScalarInt(cloudTemp, 4)) - (EMMSFC * PetscPowScalarInt(sfctemp, 4))); /* calculates flux using Stefan-Boltzmann relation */ 371c4762a1bSJed Brown PetscFunctionReturn(0); 372c4762a1bSJed Brown } 373c4762a1bSJed Brown 374c4762a1bSJed Brown PetscErrorCode calcfluxa(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar emma, PetscScalar *flux) /* this function is not currently called upon */ 375c4762a1bSJed Brown { 376c4762a1bSJed Brown PetscScalar emm = 0.001; 377c4762a1bSJed Brown 378c4762a1bSJed Brown PetscFunctionBeginUser; 379c4762a1bSJed Brown *flux = SIG * (-emm * (PetscPowScalarInt(airtemp, 4))); /* calculates flux usinge Stefan-Boltzmann relation */ 380c4762a1bSJed Brown PetscFunctionReturn(0); 381c4762a1bSJed Brown } 382d71ae5a4SJacob Faibussowitsch PetscErrorCode sensibleflux(PetscScalar sfctemp, PetscScalar airtemp, PetscScalar wind, PetscScalar *sheat) 383d71ae5a4SJacob Faibussowitsch { 384c4762a1bSJed Brown PetscScalar density = 1; /* air density */ 385c4762a1bSJed Brown PetscScalar Cp = 1005; /* heat capicity for dry air */ 386c4762a1bSJed Brown PetscScalar wndmix; /* temperature change from wind mixing: wind*Ch */ 387c4762a1bSJed Brown 388c4762a1bSJed Brown PetscFunctionBeginUser; 389c4762a1bSJed Brown wndmix = 0.0025 + 0.0042 * wind; /* regression equation valid for neutral and stable BL */ 390c4762a1bSJed Brown *sheat = density * Cp * wndmix * (airtemp - sfctemp); /* calculates sensible heat flux */ 391c4762a1bSJed Brown PetscFunctionReturn(0); 392c4762a1bSJed Brown } 393c4762a1bSJed Brown 394d71ae5a4SJacob Faibussowitsch PetscErrorCode latentflux(PetscScalar sfctemp, PetscScalar dewtemp, PetscScalar wind, PetscScalar pressure1, PetscScalar *latentheat) 395d71ae5a4SJacob Faibussowitsch { 396c4762a1bSJed Brown PetscScalar density = 1; /* density of dry air */ 397c4762a1bSJed Brown PetscScalar q; /* actual specific humitity */ 398c4762a1bSJed Brown PetscScalar qs; /* saturation specific humidity */ 399c4762a1bSJed Brown PetscScalar wndmix; /* temperature change from wind mixing: wind*Ch */ 400c4762a1bSJed Brown PetscScalar beta = .4; /* moisture availability */ 401c4762a1bSJed Brown PetscScalar mr; /* mixing ratio */ 402c4762a1bSJed Brown PetscScalar lhcnst; /* latent heat of vaporization constant = 2501000 J/kg at 0c */ 403c4762a1bSJed Brown /* latent heat of saturation const = 2834000 J/kg */ 404c4762a1bSJed Brown /* latent heat of fusion const = 333700 J/kg */ 405c4762a1bSJed Brown 406c4762a1bSJed Brown PetscFunctionBeginUser; 407c4762a1bSJed Brown wind = mph2mpers(wind); /* converts wind from mph to meters per second */ 408c4762a1bSJed Brown wndmix = 0.0025 + 0.0042 * wind; /* regression equation valid for neutral BL */ 409c4762a1bSJed Brown lhcnst = Lconst(sfctemp); /* calculates latent heat of evaporation */ 410c4762a1bSJed Brown mr = calcmixingr(sfctemp, pressure1); /* calculates saturation mixing ratio */ 411c4762a1bSJed Brown qs = calc_q(mr); /* calculates saturation specific humidty */ 412c4762a1bSJed Brown mr = calcmixingr(dewtemp, pressure1); /* calculates mixing ratio */ 413c4762a1bSJed Brown q = calc_q(mr); /* calculates specific humidty */ 414c4762a1bSJed Brown 415c4762a1bSJed Brown *latentheat = density * wndmix * beta * lhcnst * (q - qs); /* calculates latent heat flux */ 416c4762a1bSJed Brown PetscFunctionReturn(0); 417c4762a1bSJed Brown } 418c4762a1bSJed Brown 419d71ae5a4SJacob Faibussowitsch PetscErrorCode potential_temperature(PetscScalar temp, PetscScalar pressure1, PetscScalar pressure2, PetscScalar sfctemp, PetscScalar *pottemp) 420d71ae5a4SJacob Faibussowitsch { 421c4762a1bSJed Brown PetscScalar kdry; /* poisson constant for dry atmosphere */ 422c4762a1bSJed Brown PetscScalar pavg; /* average atmospheric pressure */ 423c4762a1bSJed Brown /* PetscScalar mixratio; mixing ratio */ 424c4762a1bSJed Brown /* PetscScalar kmoist; poisson constant for moist atmosphere */ 425c4762a1bSJed Brown 426c4762a1bSJed Brown PetscFunctionBeginUser; 427c4762a1bSJed Brown /* mixratio = calcmixingr(sfctemp,pressure1); */ 428c4762a1bSJed Brown 429c4762a1bSJed Brown /* initialize poisson constant */ 430c4762a1bSJed Brown kdry = 0.2854; 431c4762a1bSJed Brown /* kmoist = 0.2854*(1 - 0.24*mixratio); */ 432c4762a1bSJed Brown 433c4762a1bSJed Brown pavg = ((0.7 * pressure1) + pressure2) / 2; /* calculates simple average press */ 434c4762a1bSJed Brown *pottemp = temp * (PetscPowScalar((pressure1 / pavg), kdry)); /* calculates potential temperature */ 435c4762a1bSJed Brown PetscFunctionReturn(0); 436c4762a1bSJed Brown } 437d71ae5a4SJacob Faibussowitsch extern PetscScalar calcmixingr(PetscScalar dtemp, PetscScalar pressure1) 438d71ae5a4SJacob Faibussowitsch { 439c4762a1bSJed Brown PetscScalar e; /* vapor pressure */ 440c4762a1bSJed Brown PetscScalar mixratio; /* mixing ratio */ 441c4762a1bSJed Brown 442c4762a1bSJed Brown dtemp = dtemp - 273; /* converts from Kelvin to Celsuis */ 443c4762a1bSJed Brown e = 6.11 * (PetscPowScalar(10, ((7.5 * dtemp) / (237.7 + dtemp)))); /* converts from dew point temp to vapor pressure */ 444c4762a1bSJed Brown e = e * 100; /* converts from hPa to Pa */ 445c4762a1bSJed Brown mixratio = (0.622 * e) / (pressure1 - e); /* computes mixing ratio */ 446c4762a1bSJed Brown mixratio = mixratio * 1; /* convert to g/Kg */ 447c4762a1bSJed Brown 448c4762a1bSJed Brown return mixratio; 449c4762a1bSJed Brown } 450d71ae5a4SJacob Faibussowitsch extern PetscScalar calc_q(PetscScalar rv) 451d71ae5a4SJacob Faibussowitsch { 452c4762a1bSJed Brown PetscScalar specific_humidity; /* define specific humidity variable */ 453c4762a1bSJed Brown specific_humidity = rv / (1 + rv); /* calculates specific humidity */ 454c4762a1bSJed Brown return specific_humidity; 455c4762a1bSJed Brown } 456c4762a1bSJed Brown 457d71ae5a4SJacob Faibussowitsch PetscErrorCode calc_gflux(PetscScalar sfctemp, PetscScalar deep_grnd_temp, PetscScalar *Gflux) 458d71ae5a4SJacob Faibussowitsch { 459c4762a1bSJed Brown PetscScalar k; /* thermal conductivity parameter */ 460c4762a1bSJed Brown PetscScalar n = 0.38; /* value of soil porosity */ 461c4762a1bSJed Brown PetscScalar dz = 1; /* depth of layer between soil surface and deep soil layer */ 462c4762a1bSJed Brown PetscScalar unit_soil_weight = 2700; /* unit soil weight in kg/m^3 */ 463c4762a1bSJed Brown 464c4762a1bSJed Brown PetscFunctionBeginUser; 465c4762a1bSJed Brown k = ((0.135 * (1 - n) * unit_soil_weight) + 64.7) / (unit_soil_weight - (0.947 * (1 - n) * unit_soil_weight)); /* dry soil conductivity */ 466c4762a1bSJed Brown *Gflux = (k * (deep_grnd_temp - sfctemp) / dz); /* calculates flux from deep ground layer */ 467c4762a1bSJed Brown PetscFunctionReturn(0); 468c4762a1bSJed Brown } 469d71ae5a4SJacob Faibussowitsch extern PetscScalar emission(PetscScalar pwat) 470d71ae5a4SJacob Faibussowitsch { 471c4762a1bSJed Brown PetscScalar emma; 472c4762a1bSJed Brown 473c4762a1bSJed Brown emma = 0.725 + 0.17 * PetscLog10Real(PetscRealPart(pwat)); 474c4762a1bSJed Brown 475c4762a1bSJed Brown return emma; 476c4762a1bSJed Brown } 477d71ae5a4SJacob Faibussowitsch extern PetscScalar cloud(PetscScalar fract) 478d71ae5a4SJacob Faibussowitsch { 479c4762a1bSJed Brown PetscScalar emma = 0; 480c4762a1bSJed Brown 481c4762a1bSJed Brown /* modifies radiative balance depending on cloud cover */ 482c4762a1bSJed Brown if (fract >= 0.9) emma = 1; 483c4762a1bSJed Brown else if (0.9 > fract && fract >= 0.8) emma = 0.9; 484c4762a1bSJed Brown else if (0.8 > fract && fract >= 0.7) emma = 0.85; 485c4762a1bSJed Brown else if (0.7 > fract && fract >= 0.6) emma = 0.75; 486c4762a1bSJed Brown else if (0.6 > fract && fract >= 0.5) emma = 0.65; 487c4762a1bSJed Brown else if (0.4 > fract && fract >= 0.3) emma = emma * 1.086956; 488c4762a1bSJed Brown return emma; 489c4762a1bSJed Brown } 490d71ae5a4SJacob Faibussowitsch extern PetscScalar Lconst(PetscScalar sfctemp) 491d71ae5a4SJacob Faibussowitsch { 492c4762a1bSJed Brown PetscScalar Lheat; 493c4762a1bSJed Brown sfctemp -= 273; /* converts from kelvin to celsius */ 494c4762a1bSJed Brown Lheat = 4186.8 * (597.31 - 0.5625 * sfctemp); /* calculates latent heat constant */ 495c4762a1bSJed Brown return Lheat; 496c4762a1bSJed Brown } 497d71ae5a4SJacob Faibussowitsch extern PetscScalar mph2mpers(PetscScalar wind) 498d71ae5a4SJacob Faibussowitsch { 499c4762a1bSJed Brown wind = ((wind * 1.6 * 1000) / 3600); /* converts wind from mph to meters per second */ 500c4762a1bSJed Brown return wind; 501c4762a1bSJed Brown } 502d71ae5a4SJacob Faibussowitsch extern PetscScalar fahr_to_cel(PetscScalar temp) 503d71ae5a4SJacob Faibussowitsch { 504c4762a1bSJed Brown temp = (5 * (temp - 32)) / 9; /* converts from farhrenheit to celsuis */ 505c4762a1bSJed Brown return temp; 506c4762a1bSJed Brown } 507d71ae5a4SJacob Faibussowitsch extern PetscScalar cel_to_fahr(PetscScalar temp) 508d71ae5a4SJacob Faibussowitsch { 509c4762a1bSJed Brown temp = ((temp * 9) / 5) + 32; /* converts from celsuis to farhrenheit */ 510c4762a1bSJed Brown return temp; 511c4762a1bSJed Brown } 512d71ae5a4SJacob Faibussowitsch PetscErrorCode readinput(struct in *put) 513d71ae5a4SJacob Faibussowitsch { 514c4762a1bSJed Brown int i; 515c4762a1bSJed Brown char x; 516c4762a1bSJed Brown FILE *ifp; 517c4762a1bSJed Brown double tmp; 518c4762a1bSJed Brown 5197510d9b0SBarry Smith PetscFunctionBeginUser; 520c4762a1bSJed Brown ifp = fopen("ex5_control.txt", "r"); 5213c633725SBarry Smith PetscCheck(ifp, PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Unable to open input file"); 522ad540459SPierre Jolivet for (i = 0; i < 110; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 5233c633725SBarry Smith PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 524c4762a1bSJed Brown put->Ts = tmp; 525c4762a1bSJed Brown 526ad540459SPierre Jolivet for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 5273c633725SBarry Smith PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 528c4762a1bSJed Brown put->Td = tmp; 529c4762a1bSJed Brown 530ad540459SPierre Jolivet for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 5313c633725SBarry Smith PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 532c4762a1bSJed Brown put->Ta = tmp; 533c4762a1bSJed Brown 534ad540459SPierre Jolivet for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 5353c633725SBarry Smith PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 536c4762a1bSJed Brown put->Tc = tmp; 537c4762a1bSJed Brown 538ad540459SPierre Jolivet for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 5393c633725SBarry Smith PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 540c4762a1bSJed Brown put->fr = tmp; 541c4762a1bSJed Brown 542ad540459SPierre Jolivet for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 5433c633725SBarry Smith PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 544c4762a1bSJed Brown put->wnd = tmp; 545c4762a1bSJed Brown 546ad540459SPierre Jolivet for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 5473c633725SBarry Smith PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 548c4762a1bSJed Brown put->pwt = tmp; 549c4762a1bSJed Brown 550ad540459SPierre Jolivet for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 5513c633725SBarry Smith PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 552c4762a1bSJed Brown put->wndDir = tmp; 553c4762a1bSJed Brown 554ad540459SPierre Jolivet for (i = 0; i < 43; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 5553c633725SBarry Smith PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 556c4762a1bSJed Brown put->time = tmp; 557c4762a1bSJed Brown 558ad540459SPierre Jolivet for (i = 0; i < 63; i++) PetscCheck(fscanf(ifp, "%c", &x) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 5593c633725SBarry Smith PetscCheck(fscanf(ifp, "%lf", &tmp) == 1, PETSC_COMM_SELF, PETSC_ERR_FILE_READ, "Unable to read file"); 560c4762a1bSJed Brown put->init = tmp; 561303a5415SBarry Smith PetscFunctionReturn(0); 562c4762a1bSJed Brown } 563c4762a1bSJed Brown 564c4762a1bSJed Brown /* ------------------------------------------------------------------- */ 565d71ae5a4SJacob Faibussowitsch PetscErrorCode FormInitialSolution(DM da, Vec Xglobal, void *ctx) 566d71ae5a4SJacob Faibussowitsch { 567c4762a1bSJed Brown AppCtx *user = (AppCtx *)ctx; /* user-defined application context */ 568c4762a1bSJed Brown PetscInt i, j, xs, ys, xm, ym, Mx, My; 569c4762a1bSJed Brown Field **X; 570c4762a1bSJed Brown 571c4762a1bSJed Brown PetscFunctionBeginUser; 5729371c9d4SSatish Balay PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 573c4762a1bSJed Brown 574c4762a1bSJed Brown /* Get pointers to vector data */ 5759566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, Xglobal, &X)); 576c4762a1bSJed Brown 577c4762a1bSJed Brown /* Get local grid boundaries */ 5789566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 579c4762a1bSJed Brown 580c4762a1bSJed Brown /* Compute function over the locally owned part of the grid */ 581c4762a1bSJed Brown 582c4762a1bSJed Brown if (user->init == 1) { 583c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) { 584c4762a1bSJed Brown for (i = xs; i < xs + xm; i++) { 585c4762a1bSJed Brown X[j][i].Ts = user->Ts - i * 0.0001; 586c4762a1bSJed Brown X[j][i].Ta = X[j][i].Ts - 5; 587c4762a1bSJed Brown X[j][i].u = 0; 588c4762a1bSJed Brown X[j][i].v = 0; 589c4762a1bSJed Brown X[j][i].p = 1.25; 590c4762a1bSJed Brown if ((j == 5 || j == 6) && (i == 4 || i == 5)) X[j][i].p += 0.00001; 591c4762a1bSJed Brown if ((j == 5 || j == 6) && (i == 12 || i == 13)) X[j][i].p += 0.00001; 592c4762a1bSJed Brown } 593c4762a1bSJed Brown } 594c4762a1bSJed Brown } else { 595c4762a1bSJed Brown for (j = ys; j < ys + ym; j++) { 596c4762a1bSJed Brown for (i = xs; i < xs + xm; i++) { 597c4762a1bSJed Brown X[j][i].Ts = user->Ts; 598c4762a1bSJed Brown X[j][i].Ta = X[j][i].Ts - 5; 599c4762a1bSJed Brown X[j][i].u = 0; 600c4762a1bSJed Brown X[j][i].v = 0; 601c4762a1bSJed Brown X[j][i].p = 1.25; 602c4762a1bSJed Brown } 603c4762a1bSJed Brown } 604c4762a1bSJed Brown } 605c4762a1bSJed Brown 606c4762a1bSJed Brown /* Restore vectors */ 6079566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, Xglobal, &X)); 608c4762a1bSJed Brown PetscFunctionReturn(0); 609c4762a1bSJed Brown } 610c4762a1bSJed Brown 611c4762a1bSJed Brown /* 612c4762a1bSJed Brown RhsFunc - Evaluates nonlinear function F(u). 613c4762a1bSJed Brown 614c4762a1bSJed Brown Input Parameters: 615c4762a1bSJed Brown . ts - the TS context 616c4762a1bSJed Brown . t - current time 617c4762a1bSJed Brown . Xglobal - input vector 618c4762a1bSJed Brown . F - output vector 619c4762a1bSJed Brown . ptr - optional user-defined context, as set by SNESSetFunction() 620c4762a1bSJed Brown 621c4762a1bSJed Brown Output Parameter: 622c4762a1bSJed Brown . F - rhs function vector 623c4762a1bSJed Brown */ 624d71ae5a4SJacob Faibussowitsch PetscErrorCode RhsFunc(TS ts, PetscReal t, Vec Xglobal, Vec F, void *ctx) 625d71ae5a4SJacob Faibussowitsch { 626c4762a1bSJed Brown AppCtx *user = (AppCtx *)ctx; /* user-defined application context */ 627c4762a1bSJed Brown DM da = user->da; 628c4762a1bSJed Brown PetscInt i, j, Mx, My, xs, ys, xm, ym; 629c4762a1bSJed Brown PetscReal dhx, dhy; 630c4762a1bSJed Brown Vec localT; 631c4762a1bSJed Brown Field **X, **Frhs; /* structures that contain variables of interest and left hand side of governing equations respectively */ 632c4762a1bSJed Brown PetscScalar csoil = user->csoil; /* heat constant for layer */ 633c4762a1bSJed Brown PetscScalar dzlay = user->dzlay; /* thickness of top soil layer */ 634c4762a1bSJed Brown PetscScalar emma = user->emma; /* emission parameter */ 635c4762a1bSJed Brown PetscScalar wind = user->wind; /* wind speed */ 636c4762a1bSJed Brown PetscScalar dewtemp = user->dewtemp; /* dew point temperature (moisture in air) */ 637c4762a1bSJed Brown PetscScalar pressure1 = user->pressure1; /* sea level pressure */ 638c4762a1bSJed Brown PetscScalar airtemp = user->airtemp; /* temperature of air near boundary layer inversion */ 639c4762a1bSJed Brown PetscScalar fract = user->fract; /* fraction of the sky covered by clouds */ 640c4762a1bSJed Brown PetscScalar Tc = user->Tc; /* temperature at base of lowest cloud layer */ 641c4762a1bSJed Brown PetscScalar lat = user->lat; /* latitude */ 642c4762a1bSJed Brown PetscScalar Cp = 1005.7; /* specific heat of air at constant pressure */ 643c4762a1bSJed Brown PetscScalar Rd = 287.058; /* gas constant for dry air */ 644c4762a1bSJed Brown PetscScalar diffconst = 1000; /* diffusion coefficient */ 645c4762a1bSJed Brown PetscScalar f = 2 * 0.0000727 * PetscSinScalar(lat); /* coriolis force */ 646c4762a1bSJed Brown PetscScalar deep_grnd_temp = user->deep_grnd_temp; /* temp in lowest ground layer */ 647c4762a1bSJed Brown PetscScalar Ts, u, v, p; 648c4762a1bSJed Brown PetscScalar u_abs, u_plus, u_minus, v_abs, v_plus, v_minus; 649c4762a1bSJed Brown 650c4762a1bSJed Brown PetscScalar sfctemp1, fsfc1, Ra; 651c4762a1bSJed Brown PetscScalar sheat; /* sensible heat flux */ 652c4762a1bSJed Brown PetscScalar latentheat; /* latent heat flux */ 653c4762a1bSJed Brown PetscScalar groundflux; /* flux from conduction of deep ground layer in contact with top soil */ 654c4762a1bSJed Brown PetscInt xend, yend; 655c4762a1bSJed Brown 656c4762a1bSJed Brown PetscFunctionBeginUser; 6579566063dSJacob Faibussowitsch PetscCall(DMGetLocalVector(da, &localT)); 6589566063dSJacob Faibussowitsch PetscCall(DMDAGetInfo(da, PETSC_IGNORE, &Mx, &My, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE, PETSC_IGNORE)); 659c4762a1bSJed Brown 660c4762a1bSJed Brown dhx = (PetscReal)(Mx - 1) / (5000 * (Mx - 1)); /* dhx = 1/dx; assume 2D space domain: [0.0, 1.e5] x [0.0, 1.e5] */ 661c4762a1bSJed Brown dhy = (PetscReal)(My - 1) / (5000 * (Mx - 1)); /* dhy = 1/dy; */ 662c4762a1bSJed Brown 663c4762a1bSJed Brown /* 664c4762a1bSJed Brown Scatter ghost points to local vector,using the 2-step process 665c4762a1bSJed Brown DAGlobalToLocalBegin(),DAGlobalToLocalEnd(). 666c4762a1bSJed Brown By placing code between these two statements, computations can be 667c4762a1bSJed Brown done while messages are in transition. 668c4762a1bSJed Brown */ 6699566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalBegin(da, Xglobal, INSERT_VALUES, localT)); 6709566063dSJacob Faibussowitsch PetscCall(DMGlobalToLocalEnd(da, Xglobal, INSERT_VALUES, localT)); 671c4762a1bSJed Brown 672c4762a1bSJed Brown /* Get pointers to vector data */ 6739566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArrayRead(da, localT, &X)); 6749566063dSJacob Faibussowitsch PetscCall(DMDAVecGetArray(da, F, &Frhs)); 675c4762a1bSJed Brown 676c4762a1bSJed Brown /* Get local grid boundaries */ 6779566063dSJacob Faibussowitsch PetscCall(DMDAGetCorners(da, &xs, &ys, NULL, &xm, &ym, NULL)); 678c4762a1bSJed Brown 679c4762a1bSJed Brown /* Compute function over the locally owned part of the grid */ 680c4762a1bSJed Brown /* the interior points */ 6819371c9d4SSatish Balay xend = xs + xm; 6829371c9d4SSatish Balay yend = ys + ym; 683c4762a1bSJed Brown for (j = ys; j < yend; j++) { 684c4762a1bSJed Brown for (i = xs; i < xend; i++) { 6859371c9d4SSatish Balay Ts = X[j][i].Ts; 6869371c9d4SSatish Balay u = X[j][i].u; 6879371c9d4SSatish Balay v = X[j][i].v; 6889371c9d4SSatish Balay p = X[j][i].p; /*P = X[j][i].P; */ 689c4762a1bSJed Brown 690c4762a1bSJed Brown sfctemp1 = (double)Ts; 6919566063dSJacob Faibussowitsch PetscCall(calcfluxs(sfctemp1, airtemp, emma, fract, Tc, &fsfc1)); /* calculates surface net radiative flux */ 6929566063dSJacob Faibussowitsch PetscCall(sensibleflux(sfctemp1, airtemp, wind, &sheat)); /* calculate sensible heat flux */ 6939566063dSJacob Faibussowitsch PetscCall(latentflux(sfctemp1, dewtemp, wind, pressure1, &latentheat)); /* calculates latent heat flux */ 6949566063dSJacob Faibussowitsch PetscCall(calc_gflux(sfctemp1, deep_grnd_temp, &groundflux)); /* calculates flux from earth below surface soil layer by conduction */ 6959566063dSJacob Faibussowitsch PetscCall(calcfluxa(sfctemp1, airtemp, emma, &Ra)); /* Calculates the change in downward radiative flux */ 696c4762a1bSJed Brown fsfc1 = fsfc1 + latentheat + sheat + groundflux; /* adds radiative, sensible heat, latent heat, and ground heat flux yielding net flux */ 697c4762a1bSJed Brown 698c4762a1bSJed Brown /* convective coefficients for upwinding */ 699c4762a1bSJed Brown u_abs = PetscAbsScalar(u); 700c4762a1bSJed Brown u_plus = .5 * (u + u_abs); /* u if u>0; 0 if u<0 */ 701c4762a1bSJed Brown u_minus = .5 * (u - u_abs); /* u if u <0; 0 if u>0 */ 702c4762a1bSJed Brown 703c4762a1bSJed Brown v_abs = PetscAbsScalar(v); 704c4762a1bSJed Brown v_plus = .5 * (v + v_abs); /* v if v>0; 0 if v<0 */ 705c4762a1bSJed Brown v_minus = .5 * (v - v_abs); /* v if v <0; 0 if v>0 */ 706c4762a1bSJed Brown 707c4762a1bSJed Brown /* Solve governing equations */ 708c4762a1bSJed Brown /* P = p*Rd*Ts; */ 709c4762a1bSJed Brown 710c4762a1bSJed Brown /* du/dt -> time change of east-west component of the wind */ 711c4762a1bSJed Brown 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) */ 712c4762a1bSJed Brown - v_plus * (u - X[j - 1][i].u) * dhy - v_minus * (X[j + 1][i].u - u) * dhy /* - v(du/dy) */ 713c4762a1bSJed Brown - (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)] */ 714c4762a1bSJed Brown /* -(1/p)*(X[j][i+1].P - X[j][i-1].P)*dhx */ 715c4762a1bSJed Brown + f * v; 716c4762a1bSJed Brown 717c4762a1bSJed Brown /* dv/dt -> time change of north-south component of the wind */ 718c4762a1bSJed Brown 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) */ 719c4762a1bSJed Brown - v_plus * (v - X[j - 1][i].v) * dhy - v_minus * (X[j + 1][i].v - v) * dhy /* - v(dv/dy) */ 720c4762a1bSJed Brown - (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)] */ 721c4762a1bSJed Brown /* -(1/p)*(X[j+1][i].P - X[j-1][i].P)*dhy */ 722c4762a1bSJed Brown - f * u; 723c4762a1bSJed Brown 724c4762a1bSJed Brown /* dT/dt -> time change of temperature */ 725c4762a1bSJed Brown Frhs[j][i].Ts = (fsfc1 / (csoil * dzlay)) /* Fnet/(Cp*dz) diabatic change in T */ 726c4762a1bSJed Brown - u_plus * (Ts - X[j][i - 1].Ts) * dhx - u_minus * (X[j][i + 1].Ts - Ts) * dhx /* - u*(dTs/dx) advection x */ 727c4762a1bSJed Brown - v_plus * (Ts - X[j - 1][i].Ts) * dhy - v_minus * (X[j + 1][i].Ts - Ts) * dhy /* - v*(dTs/dy) advection y */ 728c4762a1bSJed Brown + diffconst * ((X[j][i + 1].Ts - 2 * Ts + X[j][i - 1].Ts) * dhx * dhx /* + D(Ts_xx + Ts_yy) diffusion */ 729c4762a1bSJed Brown + (X[j + 1][i].Ts - 2 * Ts + X[j - 1][i].Ts) * dhy * dhy); 730c4762a1bSJed Brown 731c4762a1bSJed Brown /* dp/dt -> time change of */ 732c4762a1bSJed Brown 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) */ 733c4762a1bSJed Brown - v_plus * (p - X[j - 1][i].p) * dhy - v_minus * (X[j + 1][i].p - p) * dhy; /* - v*(dp/dy) */ 734c4762a1bSJed Brown 735c4762a1bSJed Brown Frhs[j][i].Ta = Ra / Cp; /* dTa/dt time change of air temperature */ 736c4762a1bSJed Brown } 737c4762a1bSJed Brown } 738c4762a1bSJed Brown 739c4762a1bSJed Brown /* Restore vectors */ 7409566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArrayRead(da, localT, &X)); 7419566063dSJacob Faibussowitsch PetscCall(DMDAVecRestoreArray(da, F, &Frhs)); 7429566063dSJacob Faibussowitsch PetscCall(DMRestoreLocalVector(da, &localT)); 743c4762a1bSJed Brown PetscFunctionReturn(0); 744c4762a1bSJed Brown } 745c4762a1bSJed Brown 746d71ae5a4SJacob Faibussowitsch PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec T, void *ctx) 747d71ae5a4SJacob Faibussowitsch { 748c4762a1bSJed Brown const PetscScalar *array; 749c4762a1bSJed Brown MonitorCtx *user = (MonitorCtx *)ctx; 750c4762a1bSJed Brown PetscViewer viewer = user->drawviewer; 751c4762a1bSJed Brown PetscReal norm; 752c4762a1bSJed Brown 753c4762a1bSJed Brown PetscFunctionBeginUser; 7549566063dSJacob Faibussowitsch PetscCall(VecNorm(T, NORM_INFINITY, &norm)); 755c4762a1bSJed Brown 756c4762a1bSJed Brown if (step % user->interval == 0) { 7579566063dSJacob Faibussowitsch PetscCall(VecGetArrayRead(T, &array)); 75863a3b9bcSJacob Faibussowitsch PetscCall(PetscPrintf(PETSC_COMM_WORLD, "step %" PetscInt_FMT ", time %8.1f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f, %6.4f\n", step, (double)time, (double)(((array[0] - 273) * 9) / 5 + 32), (double)(((array[1] - 273) * 9) / 5 + 32), (double)array[2], (double)array[3], (double)array[4], (double)array[5])); 7599566063dSJacob Faibussowitsch PetscCall(VecRestoreArrayRead(T, &array)); 760c4762a1bSJed Brown } 761c4762a1bSJed Brown 7621baa6e33SBarry Smith if (user->drawcontours) PetscCall(VecView(T, viewer)); 763c4762a1bSJed Brown PetscFunctionReturn(0); 764c4762a1bSJed Brown } 765c4762a1bSJed Brown 766c4762a1bSJed Brown /*TEST 767c4762a1bSJed Brown 768c4762a1bSJed Brown build: 769c4762a1bSJed Brown requires: !complex !single 770c4762a1bSJed Brown 771c4762a1bSJed Brown test: 772c4762a1bSJed Brown args: -ts_max_steps 130 -monitor_interval 60 773c4762a1bSJed Brown output_file: output/ex5.out 774c4762a1bSJed Brown requires: !complex !single 775c4762a1bSJed Brown localrunfiles: ex5_control.txt 776c4762a1bSJed Brown 777c4762a1bSJed Brown test: 778c4762a1bSJed Brown suffix: 2 779c4762a1bSJed Brown nsize: 4 780c4762a1bSJed Brown args: -ts_max_steps 130 -monitor_interval 60 781c4762a1bSJed Brown output_file: output/ex5.out 782c4762a1bSJed Brown localrunfiles: ex5_control.txt 783c4762a1bSJed Brown requires: !complex !single 784c4762a1bSJed Brown 785*252b6434SHong Zhang # Test MatEliminateZeros() for improved FD coloring 786*252b6434SHong Zhang test: 787*252b6434SHong Zhang suffix: 3 788*252b6434SHong Zhang args: -ts_max_steps 130 -monitor_interval 60 -removezero 789*252b6434SHong Zhang requires: !complex !single 790*252b6434SHong Zhang localrunfiles: ex5_control.txt 791*252b6434SHong Zhang 792c4762a1bSJed Brown TEST*/ 793