source: palm/trunk/UTIL/chemistry/gasphase_preproc/kpp/drv/general_adj.c @ 4180

Last change on this file since 4180 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

File size: 4.2 KB
Line 
1/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2   Driver for the Adjoint (ADJ) model
3~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
4
5int  InitSaveData();
6void Initialize();
7int  SaveData();
8int  CloseSaveData();
9int  GenerateMatlab( char * prefix );
10void GetMass( KPP_REAL CL[], KPP_REAL Mass[] );
11void 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
17int 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}
Note: See TracBrowser for help on using the repository browser.