From 0377ffb6dc578f7b71e7d20478599198238cc314 Mon Sep 17 00:00:00 2001 From: vyaas Date: Tue, 18 Feb 2020 10:43:17 -0600 Subject: [PATCH] compatible with sundials-5.1.0; added command line option for sorting --- src/Makefile | 4 ++-- src/sensBVP.cpp | 49 ++++++++++++++++++++++++++++++++++------------- src/sensBrute.cpp | 3 ++- 3 files changed, 40 insertions(+), 16 deletions(-) diff --git a/src/Makefile b/src/Makefile index 37f6ff0..8f4afdc 100644 --- a/src/Makefile +++ b/src/Makefile @@ -1,7 +1,7 @@ compiler =g++ CANTERA_DIR =/opt/scientific/cantera-2.4_gnu_blas -CVODE_DIR =/opt/scientific/sundials-3.1.1_gnu_blas -KINSOL_DIR =/opt/scientific/sundials-3.1.1_gnu_blas +CVODE_DIR =/opt/scientific/sundials-5.1.0 +KINSOL_DIR =/opt/scientific/sundials-5.1.0 BVPEXE =sensBVP BRUTEEXE =sensBrute DESTDIR =~/bin diff --git a/src/sensBVP.cpp b/src/sensBVP.cpp index 35ae0a6..8ecc367 100644 --- a/src/sensBVP.cpp +++ b/src/sensBVP.cpp @@ -117,6 +117,10 @@ typedef struct { realtype T0; //Initial Temperature (in K) N_Vector Y0; //Initial Mass Fractions realtype rho0; //initial density + realtype TIgnFac; //factor that determines the ignition + //temperature: should be anywhere between 0 and 1 where 0 corresponds to the + //initial temperature and 1 corresponds to adiabatic flame + //temperature realtype TIgn; //Ignition Temperature (in K) realtype tauIgn; //Ignition delay time (in s) N_Vector x; //The grid @@ -158,6 +162,7 @@ typedef struct { int pert_index; //index of reaction whose rate is perturbed bool sensitivityAnalysis;// boolean to activate perturbations for //sensitivity analysis + bool sort; //boolean to activate sorting bool sens; //if true, perform sensitivity analysis FILE *output; //output file for temperature and species profiles FILE *ignSensOutput; //output file for ignition sensitivities @@ -293,7 +298,8 @@ int main(int argc, char *argv[]) printOutput(0.0e0,y,data); /* Create a CVode solver object for solution of the IVP: */ - mem = CVodeCreate(CV_BDF,CV_NEWTON); + //mem = CVodeCreate(CV_BDF,CV_NEWTON); + mem = CVodeCreate(CV_BDF); ier=check_flag((void *)mem, "CVodeCreate", 0); /*Associate the user data with the solver object: */ @@ -481,7 +487,8 @@ int main(int argc, char *argv[]) SUNMatrix J; int mu,ml; mu = data->nvar+2; ml = data->nvar+2; - J = SUNBandMatrix(data->KSneq, mu, ml, mu+ml); + //J = SUNBandMatrix(data->KSneq, mu, ml, mu+ml); + J = SUNBandMatrix(data->KSneq, mu, ml); ier = check_flag((void *)J, "SUNBandMatrix", 0); /* Create dense SUNLinearSolver object for use by KINSOL: */ @@ -549,7 +556,8 @@ int main(int argc, char *argv[]) SUNMatrix J; int mu,ml; mu = data->nvar+1; ml = data->nvar+1; - J = SUNBandMatrix(data->KSneq, mu, ml, mu+ml); + //J = SUNBandMatrix(data->KSneq, mu, ml, mu+ml); + J = SUNBandMatrix(data->KSneq, mu, ml); ier = check_flag((void *)J, "SUNBandMatrix", 0); /*Create a residue vector and two temporary vectors:*/ @@ -658,14 +666,17 @@ int main(int argc, char *argv[]) } if(data->sensSuccess){ - /*Sort the sensitivities in ascending order. Note the advancing - * of the beginning indices of the sensCoeffs and indices - * arrays. This is due to the Numerical recipes convention for - * array indexing used in sort subroutine:*/ - //sort(data->nreac,sensCoeffs-1,indices-1); - - /*Print out the sensitivities:*/ - printIgnSensitivities(sensCoeffs,indices,data); + + /*Sort the sensitivities in ascending order. Note the advancing + * of the beginning indices of the sensCoeffs and indices + * arrays. This is due to the Numerical recipes convention for + * array indexing used in sort subroutine:*/ + + if(data->sort){ + sort(data->nreac,sensCoeffs-1,indices-1); + } + /*Print out the sensitivities:*/ + printIgnSensitivities(sensCoeffs,indices,data); } N_VDestroy_Serial(scale); @@ -762,6 +773,10 @@ static int parseInput(UserData data, int argc, char *argv[]){ data->IVPSuccess=false; data->BVPSuccess=false; data->sensSuccess=false; + /*TIgnFac:*/ + data->TIgnFac=0.20; + /*Sorting*/ + data->sort=sort; /*****************************************************/ int ier; @@ -770,7 +785,7 @@ static int parseInput(UserData data, int argc, char *argv[]){ char comp[BUFSIZE+1]; bool enteredT0, enteredP, enteredMech, enteredComp; enteredT0=enteredP=enteredMech=enteredComp=false; - while((opt=getopt(argc,argv,"a:r:f:T:P:m:c:t:vsd")) != -1){ + while((opt=getopt(argc,argv,"a:r:f:T:P:m:c:t:q:vsod")) != -1){ switch(opt){ case 'a': data->atol=RCONST(atof(optarg)); @@ -810,9 +825,15 @@ static int parseInput(UserData data, int argc, char *argv[]){ case 's': data->sens=true; break; + case 'o': + data->sort=true; + break; case 'd': data->imposedPdt=true; break; + case 'q': + data->TIgnFac=RCONST(atof(optarg)); + break; default: printInstructions(); return(-1); @@ -887,7 +908,7 @@ static int parseInput(UserData data, int argc, char *argv[]){ data->gas->equilibrate("HP"); } data->TIgn=data->T0+(data->gas->temperature() - -data->T0)*RCONST(0.20); + -data->T0)*data->TIgnFac; printf("Ignition Temperature: %15.6e K\n",data->TIgn); data->gas->setState_TPX(data->T0, data->P, @@ -1461,6 +1482,8 @@ static void printInstructions(){ printf("\n-s :Enables ignition delay sensitivity output\n"); printf("\n-S :Enables species mass fraction sensitivity output (works only with -s)\n"); printf("\n-d :Enables manual dPdt entryspecies sensitivity output\n"); + printf("\n-o :Enables sorting of ignition delay sensitivity output\n"); + printf("\n-q :factor that determines ignition temperature (default: 0.2)\n"); printf("\nexample: "); printf("-T 1200.0 -P 1.0 -m gri30.cti"); printf(" -c H2:1.0,N2:3.76,O2:1.0"); diff --git a/src/sensBrute.cpp b/src/sensBrute.cpp index 12f3414..cda0241 100644 --- a/src/sensBrute.cpp +++ b/src/sensBrute.cpp @@ -242,7 +242,8 @@ int main(int argc, char *argv[]) printOutput(0.0e0,y,data); /* Create a CVode solver object for solution of the IVP: */ - mem = CVodeCreate(CV_BDF,CV_NEWTON); + //mem = CVodeCreate(CV_BDF,CV_NEWTON); + mem = CVodeCreate(CV_BDF); ier=check_flag((void *)mem, "CVodeCreate", 0); /*Associate the user data with the solver object: */