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

Last change on this file since 1881 was 1881, checked in by Eric.Larour, 16 years ago

Model and FemModel are now classes in their own right.
This changes the cores quite a bit.

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