PIHM: pihm.c Source File

00001 /*******************************************************************************/
00002 /*                                                                             */
00003 /*                 8888888b. 8888888 888    888 888b     d888                  */
00004 /*                 888   Y88b  888   888    888 8888b   d8888                  */
00005 /*                 888    888  888   888    888 88888b.d88888                  */
00006 /*                 888   d88P  888   8888888888 888Y88888P888                  */
00007 /*                 8888888P"   888   888    888 888 Y888P 888                  */
00008 /*                 888         888   888    888 888  Y8P  888                  */
00009 /*                 888         888   888    888 888   "   888                  */
00010 /*                 888       8888888 888    888 888       888                  */
00011 /*                                                                             */
00012 /* Version     : v2.0 (July 10, 2007)                                          */
00013 /*                                                                             */
00014 /* File        : pihm.c                                                        */
00015 /* Function    : This is main entrance of PIHM                                 */
00016 /* Programmers : Yizhong Qu   @ Pennsylvania State Univeristy                  */
00017 /*               Mukesh Kumar @ Pennsylvania State Univeristy                  */
00018 /*               Gopal Bhatt  @ Pennsylvania State Univeristy                  */
00019 /*-----------------------------------------------------------------------------*/
00020 /*                                                                             */
00021 /* PIHM is an integrated finite volume hydrologic model. It simulates channel  */
00022 /* routing, overland flow and groundwater flow in full coupled scheme. It uses */
00023 /* semi-discrete approach to discretize PDE into ODE, and solved it with       */
00024 /* SUNDIAL (CVODE 2.2.0) package [http://www.llnl.gov/CASC/sundials/]          */
00025 /*                                                                             */
00026 /* References:                                                                 */
00027 /*                                                                             */
00028 /* 1.Yizhong Qu, An Integrated Hydrological Model Using Semi-Discrete Finite   */
00029 /*               Volume Formulation, PhD Thesis, the Pennsylvanian State       */
00030 /*               University, 2004                                              */
00031 /* 2.Yizhong Qu, Christopher Duffy, Toward An Integrate Hydrological Model For */
00032 /*               River Basins: A Semi-Discrete, Finite Volume Formulation      */
00033 /*                                                                             */
00034 /*                                                                             */
00035 /* This code is free for users with research purpose only, if appropriate      */
00036 /* citation is refered. However, there is no warranty in any format for this   */
00037 /* product.                                                                    */
00038 /*                                                                             */
00039 /* For questions or comments, please contact the authors of the reference.     */
00040 /* One who want to use it for other consideration may also contact Dr. Duffy   */
00041 /* at cxd11@psu.edu.                                                           */
00042 /*                                                                             */
00043 /*******************************************************************************/
00044 
00067 /* C Header Files */
00068 #include <stdio.h>
00069 #include <stdlib.h>
00070 #include <math.h>
00071 #include <string.h>
00072 #include <time.h>
00073 
00074 
00075 
00076 /* SUNDIALS Header Files */
00077 #include "sundials_types.h"             /* contains the definition of the type realtype         */
00078 #include "sundials_dense.h"             /* defines the DenseMat type and  accessor macros       */
00079 #include "sundials_smalldense.h"        /* use generic DENSE linear solver for "small"          */
00080 #include "sundials_math.h"              /* contains UnitRoundoff, RSqrt, SQR functions          */
00081 #include "cvode.h"                      /* header file fpr CVODE                                */
00082 #include "cvode_spgmr.h"                /* the Krylov solver SPGMR in the context of CVODE      */
00083 #include "cvode_dense.h"                /* dense direct linear solver in the context of CVODE   */
00084 #include "nvector_serial.h"             /* defines the serial implementation NVECTOR-SERIAL     */
00085 
00086 
00087 
00088 /* PIHM Header Files */
00089 #include "pihm.h"                       /* Definations for all data Structure in PIHM           */
00090 //#include "et_is.h"
00091 
00092 /* Function declarations */
00093 void read_alloc(char *, Model_Data, Control_Data *);     /* read input from files :: read_alloc.c                */
00094 N_Vector N_VNew_Serial(int);                             
00095 void initialize(char *, Model_Data, Control_Data *, N_Vector);
00096                                                          /* Initialize model & Control Data :: initialize.c      */
00097 
00098 void* CVodeCreate(int , int);                                            
00099 int CVodeSetFdata(void *, void *);                       
00100 int CVodeSetInitStep(void *, realtype);                  
00101 int CVodeSetStabLimDet(void *, booleantype);             
00102 int CVodeSetMaxStep(void *, realtype);                   
00103 int CVodeMalloc(void *, CVRhsFn, realtype, N_Vector, int, realtype, void *); 
00104 int CVSpgmr(void *, int, int);                           
00105 int CVSpilsSetGSType(void *, int);                       
00107 void calET_IS(realtype, realtype, Model_Data, N_Vector); /* Calculates ET & IS    :: et_is.c                     */
00108 int CVode(void *, realtype, n_Vector, realtype *, int);  
00109 int  f(realtype, N_Vector, N_Vector, void *);            /* RHS of system of ODEs :: f.c                         */
00110 
00111 void setTSDiCounter(Model_Data mData, realtype t);       /* set the current position (iCounter) of TSD           */
00112 realtype Interpolation(TSD *Data, realtype t);           /* Data Value at time=t from a TimeSeries               */
00113 
00114 void FPrintInit(Model_Data);
00115 void FPrint(Model_Data, N_Vector, realtype);
00116 void FPrintInitFile(Model_Data, Control_Data, N_Vector, int);
00117 void FPrintCloseAll(void);
00118 
00119 
00120 
00121 
00122 
00123 
00124 int satEle, ovrEle;
00125 
00126 /* Main Function of PIHM */
00127 int main(int argc, char *argv[])
00128 {
00129     char *filename;                 /* File name prefix for input/output files    */
00130 
00131     Model_Data mData;               /* Model Data                                 */
00132     Control_Data cData;             /* Control Data                               */
00133     N_Vector CV_Y;                  /* State Variables Vector                     */
00134 
00135     void *cvode_mem;                /* pointer to the CVODE memory block          */
00136     int flag;                       /* return value of cvode function calls       */
00137 
00138     int N;                          /* Problem Size  (Numer of ODEs)              */
00139     int i,j,k;                        /* loop index variables                     */
00140     realtype t;                     /* simulation time (real time)                */
00141     realtype NextPtr, StepSize;     /* stress period & step size                  */
00142 
00143     clock_t start, end_r, end_s;    /* system clock at points                     */ //TODO: get rid of it
00144     realtype cputime_r, cputime_s;  /* for duration in realtype                   */ //TODO: get rid of it
00145 
00146     /***************************
00147     Next two lines of variable declarations are for printing flow to estuary/BC */
00148     realtype loc1_bcEle, loc_Avg_Y_Surf, loc_Avg_Y_Sub, loc_Distance, loc_Dif_Y_Sub, loc_Avg_Ksat, loc_Grad_Y_Sub, Sub_Bdd, loc_Dif_Y_Surf, loc_Grad_Y_Surf, Surf_Bdd;
00149     int loc_i, loc_j;
00150     /***************************/
00151 
00152     FILE *base2File, *over2File;
00153 
00154     char tmpFileName[100];
00155     setFileName(tmpFileName);        /* get File Name specified in calib.c        */
00156     filename = (char *)malloc(sizeof(char)*strlen(tmpFileName));
00157     strcpy(filename, tmpFileName);
00158 
00159     printf("\nBelt up!  PIHM 2.0 is starting ... \n");
00160     start = clock();
00161     /* allocate memory for model data structure */
00162     mData = (Model_Data)malloc(sizeof *mData);
00163 
00164 
00165 
00166 
00167     base2File=fopen("sc.base2","w");
00168     over2File=fopen("sc.over2","w");
00169 
00170 
00171     /* read the input files with "filename" as prefix */
00172     read_alloc(filename, mData, &cData);          /* function definition in read_alloc.c    */
00173 
00174     if(mData->UnsatMode ==1) //take off the option 1 from everywhere
00175     {
00176           N = 2*mData->NumEle + mData->NumRiv;    /* Set problem dimension                  */
00177       }
00178       if(mData->UnsatMode ==2)
00179     {
00180         N = 3*mData->NumEle + mData->NumRiv;      /* Set problem dimension                  */
00181       }
00182 
00183     CV_Y = N_VNew_Serial(N);                      /* Set Vector of initial values           */
00184 
00185 
00186     initialize(filename, mData, &cData, CV_Y);    /* initialize mode data structure         */
00187                                                   /* function definition in initialize.c    */
00188 
00189     FPrintInit(mData);
00190     //if(cData.Debug == 1) {PrintModelData(mData);}
00191 
00192     end_r = clock();
00193     cputime_r = (end_r - start)/(realtype)CLOCKS_PER_SEC;
00194 
00195     printf("\nSolving ODE system ... \n");
00196 
00197 
00198     /* Create the CVODE memory block and specify the Solution Method */
00199     cvode_mem = CVodeCreate(CV_BDF, CV_NEWTON);
00200     if(cvode_mem == NULL) { printf("CVodeCreate failed. \n"); return(1); }
00201 
00202 
00203     flag = CVodeSetFdata(cvode_mem, mData);                            /* Set Data for right-hand side function                  */
00204     flag = CVodeSetInitStep(cvode_mem,cData.InitStep);                 /* Set Initial step size                                  */
00205     flag = CVodeSetStabLimDet(cvode_mem,TRUE);                         /* ON/OFF the BDF stability limit detection algorithm     */
00206     flag = CVodeSetMaxStep(cvode_mem,cData.MaxStep);                   /* Specify the maximum absolute value of the step size    */
00207     flag = CVodeMalloc(cvode_mem, f, cData.StartTime, CV_Y, CV_SS, cData.reltol, &cData.abstol);
00208                                                                        /* provide required problem specifications,
00209                                                                          allocate internal memory for CVODE, and initialize CVODE*/
00210     flag = CVSpgmr(cvode_mem, PREC_NONE, 0);                           /* selects the CVSPGMR linear solver                      */
00211     flag = CVSpilsSetGSType(cvode_mem, MODIFIED_GS);                   /* specifies Gram-Schmidt orthogonalization to be used    */
00212 
00213 
00214     /*allocate and copy to get output file name */
00215 
00216 
00217 
00218     t = cData.StartTime;                                               /* set "t" to simulation start time                       */
00219 
00220     /* start CVODE solver in loops for NumStep number of times */
00221     for(i=0; i<cData.NumSteps; i++)
00222     {
00223         /*if (cData.Verbose != 1)
00224         {
00225             printf("  Running: %-4.1f%% ... ", (100*(i+1)/((realtype) cData.NumSteps)));
00226             fflush(stdout);
00227         }*/
00228 
00229         /* inner loops to next output points with ET/IS step size control */
00230         while(t < cData.Tout[i+1])
00231         {
00232             if (t + cData.ETStep >= cData.Tout[i+1])
00233             {
00234                 NextPtr = cData.Tout[i+1];
00235             }
00236             else
00237             {
00238                 NextPtr = t + cData.ETStep;
00239             }
00240             StepSize = NextPtr - t;
00241 
00242             calET_IS(t, StepSize, mData, CV_Y);                        /* Calculate Evaporation/Interception Rates             */
00243 
00244 
00245 /******************************************************************************************/
00246     //tempNetPrep = 0.0;
00247     //tempPrep    = 0.0;
00248     //tempET0     = 0.0;
00249     //tempTF      = 0.0;
00250     /*
00251     for(tempI=0; tempI<mData->NumEle; tempI++){
00252         tempPrep += mData->ElePrep[tempI];
00253         tempNetPrep += mData->EleNetPrep[tempI];
00254         tempET0 += mData->EleET[tempI][0]*mData->Ele[tempI].VegFrac;
00255         tempET1 += mData->EleET[tempI][1];
00256         tempET2 += mData->EleET[tempI][2];
00257         tempTF  += mData->EleTF[tempI]*mData->Ele[tempI].VegFrac;
00258     }
00259     fprintf(prepFile, "%lf\t%lf\t%lf\t%lf\t%lf\t%lf\t%lf\n", t, tempPrep/mData->NumEle, tempNetPrep/mData->NumEle, tempET0/mData->NumEle, tempET1/mData->NumEle, tempET2/mData->NumEle,tempTF/mData->NumEle);//fflush(prepFile);
00260 */
00261 /******************************************************************************************/
00262 
00263             Tsteps=t;
00264             printf("\n Tsteps = %f ",t);
00265 
00266             flag = CVode(cvode_mem, NextPtr, CV_Y, &t, CV_NORMAL);    /* Advance solution in time                                */
00267 
00268             setTSDiCounter(mData, t);
00269             FPrint(mData, CV_Y, t);
00270 
00271 
00272 /***************************************************************************************************************/
00273 /*  Specially for Rhode :: Estuary Discharge                                                                   */
00274 /***************************************************************************************************************/
00275     for(loc_i=0; loc_i<mData->NumEle; loc_i++){
00276         for(loc_j=0; loc_j<3; loc_j++){
00277             if(mData->Ele[loc_i].BC > 0){         // Dirichlet BC
00278                 loc1_bcEle = Interpolation(&mData->TSD_EleBC[(mData->Ele[loc_i].BC)-1], t);
00279                 if( loc1_bcEle < mData->Ele[loc_i].zmin ){
00280                     loc_Avg_Y_Surf = NV_Ith_S(CV_Y, loc_i)/2; //(DummyY[i])/2;
00281                     loc_Avg_Y_Sub  = NV_Ith_S(CV_Y, loc_i+2*mData->NumEle);//(DummyY[i+2*mData->NumEle])/2;
00282                 }
00283                 else if(loc1_bcEle < mData->Ele[loc_i].zmax){
00284                     loc_Avg_Y_Surf = NV_Ith_S(CV_Y, loc_i)/2;//(DummyY[i])/2;
00285                     loc_Avg_Y_Sub  = (loc1_bcEle - mData->Ele[loc_i].zmin + NV_Ith_S(CV_Y, loc_i+2*mData->NumEle))/2.0;
00286                 }
00287                 else{
00288                     loc_Avg_Y_Surf = (NV_Ith_S(CV_Y, loc_i) + loc1_bcEle - mData->Ele[loc_i].zmax)/2;
00289                     loc_Avg_Y_Sub  = mData->Ele[loc_i].zmax - mData->Ele[loc_i].zmin;
00290                 }
00291                 loc_Distance = sqrt(pow(mData->Ele[loc_i].edge[0]*mData->Ele[loc_i].edge[1]*mData->Ele[loc_i].edge[2]/(4*mData->Ele[loc_i].area), 2) - pow(mData->Ele[loc_i].edge[loc_j]/2, 2));
00292                 //Sub
00293                 loc_Dif_Y_Sub = NV_Ith_S(CV_Y, loc_i+2*mData->NumEle) + mData->Ele[loc_i].zmin - loc1_bcEle;
00294                 loc_Avg_Ksat = mData->Ele[loc_i].Ksat;
00295                 loc_Grad_Y_Sub = loc_Dif_Y_Sub/loc_Distance;
00296                 Sub_Bdd = loc_Avg_Ksat * loc_Grad_Y_Sub * loc_Avg_Y_Sub * mData->Ele[loc_i].edge[loc_j];
00297 
00298                 //Surf
00299                 loc_Dif_Y_Surf = NV_Ith_S(CV_Y, loc_i) + mData->Ele[loc_i].zmax - loc1_bcEle;
00300                 loc_Grad_Y_Surf = loc_Dif_Y_Surf / loc_Distance;
00301                 //Surf_Bdd = 0.1*(Grad_Y_Surf>0?1:-1)* (Avg_Y_Sub * mData->Ele[i].edge[loc_j] ) * (pow(pow(Avg_Y_Surf, 1.0/3.0),2)/(mData->Ele[i].Rough)) * sqrt((Grad_Y_Surf>0?1:-1)*Grad_Y_Surf);
00302 
00303                 Surf_Bdd = sqrt(NV_Ith_S(CV_Y, loc_i)/loc_Distance)*pow(NV_Ith_S(CV_Y, loc_i),2.0/3.0)*(NV_Ith_S(CV_Y, loc_i)/2)*mData->Ele[loc_i].edge[loc_j]/mData->Ele[loc_i].Rough;
00304 
00305                 fprintf(base2File, "%lf\t", Sub_Bdd);
00306                 fprintf(over2File, "%lf\t", Surf_Bdd);
00307             }
00308         }
00309     }
00310     fprintf(base2File, "\n");
00311     fprintf(over2File, "\n");
00312 
00313 /***************************************************************************************************************/
00314 
00315 
00316 
00317 
00318 
00319             satEle=0;
00320             ovrEle=0;
00321             for(k=0; k<mData->NumEle;k++)
00322             {
00323                 if(NV_Ith_S(CV_Y, k+2*mData->NumEle)/(mData->Ele[k].zmax-mData->Ele[k].zmin)>0.99){
00324                     satEle++;
00325                 }
00326                 if(NV_Ith_S(CV_Y, k)>1E-4){
00327                     ovrEle++;
00328                 }
00329             }
00330 
00331 
00332             //fprintf(res_flux_file,"\n"); //fflush(res_flux_file);
00333 
00334 
00335 /**********************************************MASS BALANCE*****************************************************/
00336 
00337         }
00338 
00339 
00340         /* print out results to files at every output time */
00341            //if (cData.res_out == 1) {FPrintY(mData, CV_Y, t, res_state_file);}
00342         //if (cData.flux_out == 1) {FPrintFlux(mData, t, res_flux_file);}
00343         //if (cData.etis_out == 1) {FPrintETIS(mData, t, res_etis_file);}
00344           /*  if (cData.q_out == 1) {FPrintQ(mData, t, res_q_file);} */
00345 
00346         /*    if (cData.Verbose != 1) {printf("\r");}
00347 
00348         if(flag != SUCCESS) {printf("CVode failed, flag = %d. \n", flag); return(flag);}
00349           */
00350         /* clear buffer */
00351            fflush(stdout);
00352        }
00353 
00354     FPrintInitFile(mData, cData, CV_Y, i);                    /* Routine for .init File : print.c     */
00355     FPrintCloseAll();
00356 
00357 
00358 
00359 
00360 
00361 
00362 
00363 
00364   /* capture time */
00365     end_s = clock();
00366     cputime_s = (end_s - end_r)/(realtype)CLOCKS_PER_SEC;
00367 
00368     /* print out simulation statistics */
00369     /*PrintFarewell(cData, iopt, ropt, cputime_r, cputime_s);
00370     if (cData.res_out == 1) {FPrintFarewell(cData, res_state_file, iopt, ropt, cputime_r, cputime_s);}
00371     */
00372     /* close output files */
00373     //if (cData.res_out == 1)  {fclose(res_state_file);}
00374     //if (cData.flux_out == 1) {fclose(res_flux_file);}
00375     //if (cData.etis_out == 1) {fclose(res_etis_file);}
00376     //if (cData.q_out == 1)    {fclose(res_q_file);}
00377 
00378     free(mData);
00379 
00380     return 0;
00381 }
00382 
00383 
00384 
00385 
00386 
00387 
00388 
00389 
00390 void setTSDiCounter(Model_Data mData, realtype t)
00392 
00395 {
00396 
00397     int k;
00398 
00399     for(k=0; k<mData->NumPrep; k++)
00400     {
00401         while(mData->TSD_Prep[k].iCounter < mData->TSD_Prep[k].length && t/(24.0*60.0) > mData->TSD_Prep[k].TS[mData->TSD_Prep[k].iCounter+1][0]){
00402             mData->TSD_Prep[k].iCounter++;
00403         }
00404     }
00405     for(k=0; k<mData->NumTemp; k++){
00406         while(mData->TSD_Temp[k].iCounter < mData->TSD_Temp[k].length && t/(24.0*60.0) > mData->TSD_Temp[k].TS[mData->TSD_Temp[k].iCounter+1][0]){
00407             mData->TSD_Temp[k].iCounter++;
00408         }
00409     }
00410     for(k=0; k<mData->NumHumidity; k++){
00411         while(mData->TSD_Humidity[k].iCounter < mData->TSD_Humidity[k].length && t/(24.0*60.0) > mData->TSD_Humidity[k].TS[mData->TSD_Humidity[k].iCounter+1][0]){
00412             mData->TSD_Humidity[k].iCounter++;
00413         }
00414     }
00415     for(k=0; k<mData->NumWindVel; k++){
00416         while(mData->TSD_WindVel[k].iCounter < mData->TSD_WindVel[k].length && t/(24.0*60.0) > mData->TSD_WindVel[k].TS[mData->TSD_WindVel[k].iCounter+1][0]){
00417             mData->TSD_WindVel[k].iCounter++;
00418         }
00419     }
00420     for(k=0; k<mData->NumRn; k++){
00421         while(mData->TSD_Rn[k].iCounter < mData->TSD_Rn[k].length && t/(24.0*60.0) > mData->TSD_Rn[k].TS[mData->TSD_Rn[k].iCounter+1][0]){
00422             mData->TSD_Rn[k].iCounter++;
00423         }
00424     }
00425     for(k=0; k<mData->NumG; k++){
00426         while(mData->TSD_G[k].iCounter < mData->TSD_G[k].length && t/(24.0*60.0) > mData->TSD_G[k].TS[mData->TSD_G[k].iCounter+1][0]){
00427             mData->TSD_G[k].iCounter++;
00428         }
00429     }
00430     for(k=0; k<mData->NumP; k++){
00431         while(mData->TSD_Pressure[k].iCounter < mData->TSD_Pressure[k].length && t/(24.0*60.0) > mData->TSD_Pressure[k].TS[mData->TSD_Pressure[k].iCounter+1][0]){
00432             mData->TSD_Pressure[k].iCounter++;
00433         }
00434     }
00435     for(k=0; k<mData->NumLC; k++){
00436         while(mData->TSD_LAI[k].iCounter < mData->TSD_LAI[k].length && t/(24.0*60.0) > mData->TSD_LAI[k].TS[mData->TSD_LAI[k].iCounter+1][0]){
00437             mData->TSD_LAI[k].iCounter++;
00438         }
00439     }
00440     for(k=0; k<mData->NumLC; k++){
00441         while(mData->TSD_DH[k].iCounter < mData->TSD_DH[k].length && t/(24.0*60.0) > mData->TSD_DH[k].TS[mData->TSD_DH[k].iCounter+1][0]){
00442             mData->TSD_DH[k].iCounter++;
00443         }
00444     }
00445     for(k=0; k<mData->NumMeltF; k++){
00446         while(mData->TSD_MeltF[k].iCounter < mData->TSD_MeltF[k].length && t/(24.0*60.0) > mData->TSD_MeltF[k].TS[mData->TSD_MeltF[k].iCounter+1][0]){
00447             mData->TSD_MeltF[k].iCounter++;
00448         }
00449     }
00450     for(k=0; k<mData->NumSource; k++){
00451         while(mData->TSD_Source[k].iCounter < mData->TSD_Source[k].length && t/(24.0*60.0) > mData->TSD_Source[k].TS[mData->TSD_Source[k].iCounter+1][0]){
00452             mData->TSD_Source[k].iCounter++;
00453         }
00454     }
00455 
00456     for(k=0; k<mData->Num1BC+mData->Num2BC; k++){
00457         while(mData->TSD_EleBC[k].iCounter < mData->TSD_EleBC[k].length && t/(24.0*60.0) > mData->TSD_EleBC[k].TS[mData->TSD_EleBC[k].iCounter+1][0]){
00458             mData->TSD_EleBC[k].iCounter++;
00459         }
00460     }
00461 
00462 }