[2696] | 1 | /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ |
---|
| 2 | Driver for the Adjoint (ADJ) model |
---|
| 3 | ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/ |
---|
| 4 | |
---|
| 5 | int InitSaveData(); |
---|
| 6 | void Initialize(); |
---|
| 7 | int SaveData(); |
---|
| 8 | int CloseSaveData(); |
---|
| 9 | int GenerateMatlab( char * prefix ); |
---|
| 10 | void GetMass( KPP_REAL CL[], KPP_REAL Mass[] ); |
---|
| 11 | void INTEGRATE_ADJ( int NADJ, KPP_REAL Y[], KPP_REAL Lambda[][NVAR], |
---|
| 12 | KPP_REAL TIN, KPP_REAL TOUT, KPP_REAL ATOL_adj[][NVAR], |
---|
| 13 | KPP_REAL RTOL_adj[][NVAR], int ICNTRL_U[], |
---|
| 14 | KPP_REAL RCNTRL_U[], int ISTATUS_U[], |
---|
| 15 | KPP_REAL RSTATUS_U[] ); |
---|
| 16 | |
---|
| 17 | int main() { |
---|
| 18 | |
---|
| 19 | KPP_REAL T, DVAL[NSPEC]; |
---|
| 20 | int i, j, ind_1 = ind_O1D, ind_2 = ind_O3; |
---|
| 21 | |
---|
| 22 | /*~~> NADJ = Number of functionals for which sensitivities are computed |
---|
| 23 | Note: the setting below is for sensitivities of all final concentrations; |
---|
| 24 | the setting may have to be changed for other applications */ |
---|
| 25 | int NADJ = NVAR; |
---|
| 26 | KPP_REAL Y_adj[NADJ][NVAR], ATOL_adj[NADJ][NVAR], RTOL_adj[NADJ][NVAR]; |
---|
| 27 | |
---|
| 28 | /*~~> Control (in) and status (out) arguments for the integration */ |
---|
| 29 | KPP_REAL RCNTRL[20], RSTATUS[20]; |
---|
| 30 | int ICNTRL[20], ISTATUS[20]; |
---|
| 31 | |
---|
| 32 | STEPMIN = (double)0.0; |
---|
| 33 | STEPMAX = (double)0.0; |
---|
| 34 | |
---|
| 35 | /*~~~> Tolerances for calculating concentrations */ |
---|
| 36 | for( i=0; i<NVAR; i++) { |
---|
| 37 | RTOL[i] = 1.0e-4; |
---|
| 38 | ATOL[i] = 1.0e-3; |
---|
| 39 | } |
---|
| 40 | |
---|
| 41 | /*~~~> Tolerances for calculating adjoints |
---|
| 42 | are used for controlling adjoint truncation error |
---|
| 43 | and for solving the linear adjoint equations by iterations |
---|
| 44 | Note: Adjoints typically span many orders of magnitude |
---|
| 45 | and a careful tuning of ATOL_adj may be necessary */ |
---|
| 46 | for(i=0; i<NADJ; i++) { |
---|
| 47 | for(j=0; j<NVAR; j++) { |
---|
| 48 | RTOL_adj[i][j] = 1.0e-4; |
---|
| 49 | ATOL_adj[i][j] = 1.0e-10; |
---|
| 50 | } |
---|
| 51 | } |
---|
| 52 | |
---|
| 53 | Initialize(); |
---|
| 54 | |
---|
| 55 | /*~~~> The adjoint values at the final time */ |
---|
| 56 | for(i=0; i<NADJ; i++) { |
---|
| 57 | for(j=0; j<NVAR; j++) |
---|
| 58 | Y_adj[i][j] = (double)0.0; |
---|
| 59 | } |
---|
| 60 | for(j=0; j<NADJ; j++) |
---|
| 61 | Y_adj[j][j] = (double)1.0; |
---|
| 62 | |
---|
| 63 | /*~~~> Default control options */ |
---|
| 64 | for(i=0; i<20; i++) { |
---|
| 65 | ICNTRL[i] = 0; |
---|
| 66 | RCNTRL[i] = (double)0.0; |
---|
| 67 | ISTATUS[i] = 0; |
---|
| 68 | RSTATUS[i] = (double)0.0; |
---|
| 69 | } |
---|
| 70 | |
---|
| 71 | /*~~~> Begin time loop */ |
---|
| 72 | |
---|
| 73 | TIME = TSTART; |
---|
| 74 | InitSaveData(); |
---|
| 75 | |
---|
| 76 | T = TSTART; |
---|
| 77 | |
---|
| 78 | GetMass( C, DVAL ); |
---|
| 79 | printf("\n%6.1f%% %7.2f ", (TIME-TSTART)/(TEND-TSTART)*100, TIME/3600 ); |
---|
| 80 | for( i = 0; i < NMONITOR; i++ ) |
---|
| 81 | printf( "%9.3e ", C[ MONITOR[i] ]/CFACTOR ); |
---|
| 82 | for( i = 0; i < NMASS; i++ ) |
---|
| 83 | printf( "%9.3e ", DVAL[i]/CFACTOR ); |
---|
| 84 | |
---|
| 85 | TIME = T; |
---|
| 86 | SaveData(); |
---|
| 87 | |
---|
| 88 | INTEGRATE_ADJ( NADJ, VAR, Y_adj, T, TEND, ATOL_adj, RTOL_adj, ICNTRL, |
---|
| 89 | RCNTRL, ISTATUS, RSTATUS ); |
---|
| 90 | |
---|
| 91 | GetMass( C, DVAL ); |
---|
| 92 | |
---|
| 93 | printf("\n%6.1f%% %7.2f ", (TEND-TSTART)/(TEND-TSTART)*100, TIME/3600 ); |
---|
| 94 | for( i = 0; i < NMONITOR; i++ ) |
---|
| 95 | printf( "%9.3e ", C[ MONITOR[i] ]/CFACTOR ); |
---|
| 96 | for( i = 0; i < NMASS; i++ ) |
---|
| 97 | printf( "%9.3e ", DVAL[i]/CFACTOR ); |
---|
| 98 | |
---|
| 99 | TIME = T; |
---|
| 100 | SaveData(); |
---|
| 101 | |
---|
| 102 | /*~~~> End time loop ~~~~~~~~~~*/ |
---|
| 103 | |
---|
| 104 | printf( "\n\n****************************************************\n" ); |
---|
| 105 | printf( " Concentrations and Sensitivities at final time\n" ); |
---|
| 106 | printf( " were written in the file KPP_ROOT_ADJ_results.m\n"); |
---|
| 107 | printf( "****************************************************\n"); |
---|
| 108 | |
---|
| 109 | FILE *out; |
---|
| 110 | out = fopen("KPP_ROOT_ADJ_results.m", "w"); |
---|
| 111 | if(out == NULL) { |
---|
| 112 | printf("Unable to open file KPP_ROOT_ADJ_results.m\n"); |
---|
| 113 | exit(1); |
---|
| 114 | } |
---|
| 115 | |
---|
| 116 | for(j=0; j<NADJ; j++) { |
---|
| 117 | for(i=0; i<NVAR; i++) |
---|
| 118 | fprintf( out, "%7.6e ", Y_adj[j][i] ); |
---|
| 119 | fprintf( out, "\n" ); |
---|
| 120 | } |
---|
| 121 | |
---|
| 122 | fclose(out); |
---|
| 123 | |
---|
| 124 | printf( "\nADJ: d[%s](tf)/d[%s](t0) = %7.6e\n", SPC_NAMES[ind_1], |
---|
| 125 | SPC_NAMES[ind_1], Y_adj[ind_1][ind_1] ); |
---|
| 126 | printf( "ADJ: d[%s](tf)/d[%s](t0) = %7.6e\n", SPC_NAMES[ind_2], |
---|
| 127 | SPC_NAMES[ind_2], Y_adj[ind_2][ind_2] ); |
---|
| 128 | printf( "ADJ: d[%s](tf)/d[%s](t0) = %7.6e\n", SPC_NAMES[ind_2], |
---|
| 129 | SPC_NAMES[ind_1], Y_adj[ind_2][ind_1] ); |
---|
| 130 | printf( "ADJ: d[%s](tf)/d[%s](t0) = %7.6e\n\n", SPC_NAMES[ind_1], |
---|
| 131 | SPC_NAMES[ind_2], Y_adj[ind_1][ind_2] ); |
---|
| 132 | |
---|
| 133 | CloseSaveData(); |
---|
| 134 | |
---|
| 135 | /*~~~> The entire matrix of sensitivities */ |
---|
| 136 | printf(" "); |
---|
| 137 | for(i=0; i<NVAR; i++) |
---|
| 138 | printf( "%15s", SPC_NAMES[i] ); |
---|
| 139 | printf("\n"); |
---|
| 140 | for(j=0; j<NVAR; j++) { |
---|
| 141 | printf("%5s = ", SPC_NAMES[j]); |
---|
| 142 | for(i=0; i<NADJ; i++) |
---|
| 143 | printf( " %7.6e ", Y_adj[i][j] ); |
---|
| 144 | printf("\n"); |
---|
| 145 | } |
---|
| 146 | |
---|
| 147 | return 0; |
---|
| 148 | } |
---|