source: issm/trunk/src/c/parallel/transient.cpp@ 3567

Last change on this file since 3567 was 3567, checked in by Mathieu Morlighem, 15 years ago

all enums are now C Enums, no need to provide a number: the numbering is done automatically
File size: 5.3 KB
Line 
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
15int 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}
Note: See TracBrowser for help on using the repository browser.