1 | /*!\file: transient.cpp
|
---|
2 | * \brief: transient solution
|
---|
3 | */
|
---|
4 |
|
---|
5 | #include "../issm.h"
|
---|
6 | #include "./parallel.h"
|
---|
7 |
|
---|
8 | #ifdef HAVE_CONFIG_H
|
---|
9 | #include "config.h"
|
---|
10 | #else
|
---|
11 | #error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!"
|
---|
12 | #endif
|
---|
13 |
|
---|
14 |
|
---|
15 | int main(int argc,char* *argv){
|
---|
16 |
|
---|
17 | /*I/O: */
|
---|
18 | FILE* fid=NULL;
|
---|
19 | char* inputfilename=NULL;
|
---|
20 | char* outputfilename=NULL;
|
---|
21 | char* lockname=NULL;
|
---|
22 | int numberofnodes;
|
---|
23 | int qmu_analysis=0;
|
---|
24 | double waitonlock=0;
|
---|
25 |
|
---|
26 | /*Model: */
|
---|
27 | Model* model=NULL;
|
---|
28 | int dim=-1;
|
---|
29 |
|
---|
30 | /*Results: */
|
---|
31 | DataSet* results=NULL;
|
---|
32 | DataSet* processed_results=NULL;
|
---|
33 | Result* result=NULL;
|
---|
34 |
|
---|
35 | ParameterInputs* inputs=NULL;
|
---|
36 |
|
---|
37 | /*inputs: */
|
---|
38 | double* u_g=NULL;
|
---|
39 | double* m_g=NULL;
|
---|
40 | double* a_g=NULL;
|
---|
41 | double dt;
|
---|
42 | double yts;
|
---|
43 | Param* param=NULL;
|
---|
44 |
|
---|
45 | /*time*/
|
---|
46 | double start, finish;
|
---|
47 | double start_core, finish_core;
|
---|
48 | double start_init, finish_init;
|
---|
49 |
|
---|
50 | MODULEBOOT();
|
---|
51 |
|
---|
52 | #if !defined(_PARALLEL_) || (defined(_PARALLEL_) && !defined(_HAVE_PETSC_))
|
---|
53 | ISSMERROR(" parallel executable was compiled without support of parallel libraries!");
|
---|
54 | #endif
|
---|
55 |
|
---|
56 | /*Initialize Petsc and get start time*/
|
---|
57 | PetscInitialize(&argc,&argv,(char *)0,"");
|
---|
58 | MPI_Barrier(MPI_COMM_WORLD); start=MPI_Wtime();
|
---|
59 |
|
---|
60 | /*Size and rank: */
|
---|
61 | MPI_Comm_rank(MPI_COMM_WORLD,&my_rank);
|
---|
62 | MPI_Comm_size(MPI_COMM_WORLD,&num_procs);
|
---|
63 |
|
---|
64 | _printf_("recover , input file name and output file name:\n");
|
---|
65 | inputfilename=argv[2];
|
---|
66 | outputfilename=argv[3];
|
---|
67 | lockname=argv[4];
|
---|
68 |
|
---|
69 | /*Open handle to data on disk: */
|
---|
70 | fid=pfopen(inputfilename,"rb");
|
---|
71 |
|
---|
72 | /*Initialize model structure: */
|
---|
73 | MPI_Barrier(MPI_COMM_WORLD); start_init=MPI_Wtime();
|
---|
74 | model=new Model();
|
---|
75 |
|
---|
76 | _printf_("read and create finite element model:\n");
|
---|
77 | _printf_("\n reading diagnostic horiz model data:\n");
|
---|
78 | model->AddFormulation(fid,DiagnosticAnalysisEnum,HorizAnalysisEnum);
|
---|
79 |
|
---|
80 | _printf_("\n reading diagnostic vert model data:\n");
|
---|
81 | model->AddFormulation(fid,DiagnosticAnalysisEnum,VertAnalysisEnum);
|
---|
82 |
|
---|
83 | _printf_("\n reading diagnostic stokes model data:\n");
|
---|
84 | model->AddFormulation(fid,DiagnosticAnalysisEnum,StokesAnalysisEnum);
|
---|
85 |
|
---|
86 | _printf_("\n reading diagnostic hutter model data:\n");
|
---|
87 | model->AddFormulation(fid,DiagnosticAnalysisEnum,HutterAnalysisEnum);
|
---|
88 |
|
---|
89 | _printf_("\n reading surface and bed slope computation model data:\n");
|
---|
90 | model->AddFormulation(fid,SlopecomputeAnalysisEnum);
|
---|
91 |
|
---|
92 | _printf_("\n reading prognositc model data:\n");
|
---|
93 | model->AddFormulation(fid,PrognosticAnalysisEnum);
|
---|
94 |
|
---|
95 | /*Do we run in 3d?, in which case we need thermal and melting also:*/
|
---|
96 | model->FindParam(&dim,"dim");
|
---|
97 | if(dim==3){
|
---|
98 | _printf_("read and create thermal finite element model:\n");
|
---|
99 | model->AddFormulation(fid,ThermalAnalysisEnum,TransientAnalysisEnum);
|
---|
100 | _printf_("read and create melting finite element model:\n");
|
---|
101 | model->AddFormulation(fid,MeltingAnalysisEnum,TransientAnalysisEnum);
|
---|
102 | }
|
---|
103 |
|
---|
104 | /*recover parameters: */
|
---|
105 | model->FindParam(&waitonlock,"waitonlock");
|
---|
106 | model->FindParam(&qmu_analysis,"qmu_analysis");
|
---|
107 |
|
---|
108 | _printf_("initialize inputs:\n");
|
---|
109 | model->FindParam(&u_g,NULL,NULL,"u_g",PrognosticAnalysisEnum);
|
---|
110 | model->FindParam(&m_g,NULL,NULL,"m_g",PrognosticAnalysisEnum);
|
---|
111 | model->FindParam(&a_g,NULL,NULL,"a_g",PrognosticAnalysisEnum);
|
---|
112 | model->FindParam(&numberofnodes,"numberofnodes");
|
---|
113 | model->FindParam(&dt,"dt");
|
---|
114 | model->FindParam(&yts,"yts");
|
---|
115 |
|
---|
116 |
|
---|
117 | inputs=new ParameterInputs;
|
---|
118 | inputs->Add("velocity",u_g,3,numberofnodes);
|
---|
119 | inputs->Add("melting",m_g,1,numberofnodes);
|
---|
120 | inputs->Add("accumulation",a_g,1,numberofnodes);
|
---|
121 | inputs->Add("dt",dt*yts);
|
---|
122 |
|
---|
123 | _printf_("initialize results:\n");
|
---|
124 | results=new DataSet(ResultsEnum);
|
---|
125 | MPI_Barrier(MPI_COMM_WORLD); finish_init=MPI_Wtime();
|
---|
126 |
|
---|
127 | /*are we running the solution sequence, or a qmu wrapper around it? : */
|
---|
128 | if(!qmu_analysis){
|
---|
129 |
|
---|
130 | /*run diagnostic analysis: */
|
---|
131 | _printf_("call computational core:\n");
|
---|
132 | MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
|
---|
133 | transient_core(results,model,inputs);
|
---|
134 | MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
|
---|
135 |
|
---|
136 | /*Add analysis_type to results: */
|
---|
137 | result=new Result(results->Size()+1,0,1,"analysis_type","transient");
|
---|
138 | results->AddObject(result);
|
---|
139 |
|
---|
140 | _printf_("process results:\n");
|
---|
141 | ProcessResults(&processed_results,results,model,TransientAnalysisEnum);
|
---|
142 |
|
---|
143 | _printf_("write results to disk:\n");
|
---|
144 | OutputResults(processed_results,outputfilename);
|
---|
145 | }
|
---|
146 | else{
|
---|
147 | /*run qmu analysis: */
|
---|
148 | _printf_("calling qmu analysis on transient core:\n");
|
---|
149 |
|
---|
150 | #ifdef _HAVE_DAKOTA_
|
---|
151 | MPI_Barrier(MPI_COMM_WORLD); start_core=MPI_Wtime( );
|
---|
152 | Qmux(model,inputs,TransientAnalysisEnum,NoneAnalysisEnum);
|
---|
153 | MPI_Barrier(MPI_COMM_WORLD); finish_core=MPI_Wtime( );
|
---|
154 | #else
|
---|
155 | ISSMERROR(" Dakota not present, cannot do qmu!");
|
---|
156 | #endif
|
---|
157 | }
|
---|
158 |
|
---|
159 | if (waitonlock>0){
|
---|
160 | _printf_("write lock file:\n");
|
---|
161 | WriteLockFile(lockname);
|
---|
162 | }
|
---|
163 |
|
---|
164 | /*Free ressources:*/
|
---|
165 | delete results;
|
---|
166 | delete processed_results;
|
---|
167 | delete model;
|
---|
168 | delete inputs;
|
---|
169 | xfree((void**)&u_g);
|
---|
170 | xfree((void**)&m_g);
|
---|
171 | xfree((void**)&a_g);
|
---|
172 |
|
---|
173 | /*Get finish time and close*/
|
---|
174 | MPI_Barrier(MPI_COMM_WORLD); finish = MPI_Wtime( );
|
---|
175 | _printf_("\n %-34s %f seconds \n","Model initialization elapsed time:",finish_init-start_init);
|
---|
176 | _printf_(" %-34s %f seconds \n","Core solution elapsed time:",finish_core-start_core);
|
---|
177 | _printf_("\n %s %i hrs %i min %i sec\n\n","Total elapsed time:",int((finish-start)/3600),int(int(finish-start)%3600/60),int(finish-start)%60);
|
---|
178 | _printf_("closing MPI and Petsc\n");
|
---|
179 | PetscFinalize();
|
---|
180 |
|
---|
181 | /*end module: */
|
---|
182 | MODULEEND();
|
---|
183 |
|
---|
184 | return 0; //unix success return;
|
---|
185 | }
|
---|