Ice Sheet System Model  4.18
Code documentation
solutionsequence_fct.cpp
Go to the documentation of this file.
1 
5 #include "../toolkits/toolkits.h"
6 #include "../classes/classes.h"
7 #include "../shared/shared.h"
8 #include "../modules/modules.h"
9 #include "../analyses/analyses.h"
10 #define USEPENALTYMETHOD false
11 
12 #ifdef _HAVE_PETSC_
13 void CreateDMatrix(Mat* pD,Mat K){/*{{{*/
14  /*Create D matrix such that:
15  *
16  * d_ij = max( -k_ij,0,-k_ji) off diagonal
17  *
18  * d_ii = - sum_{i!=j} d_ij for the diagonal
19  *
20  */
21 
22  /*Intermediaries*/
23  int ncols,ncols2,rstart,rend;
24  double d,diagD;
25  Mat D = NULL;
26  Mat K_transp = NULL;
27  int* cols = NULL;
28  int* cols2 = NULL;
29  double* vals = NULL;
30  double* vals2 = NULL;
31 
32  /*First, we need to transpose K so that we access both k_ij and k_ji*/
33  MatTranspose(K,MAT_INITIAL_MATRIX,&K_transp);
34 
35  /*Initialize output (D has the same non zero pattern as K)*/
36  MatDuplicate(K,MAT_SHARE_NONZERO_PATTERN,&D);
37 
38  /*Go through the rows of K an K' and build D*/
39  MatGetOwnershipRange(K,&rstart,&rend);
40  for(int row=rstart; row<rend; row++){
41  diagD = 0.;
42  MatGetRow(K ,row,&ncols, (const int**)&cols, (const double**)&vals);
43  MatGetRow(K_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
44  _assert_(ncols==ncols2);
45  for(int j=0; j<ncols; j++) {
46  _assert_(cols[j]==cols2[j]);
47  d = max(max(-vals[j],-vals2[j]),0.);
48  MatSetValue(D,row,cols[j],d,INSERT_VALUES);
49  if(cols[j]!=row) diagD -= d;
50  }
51  MatSetValue(D,row,row,diagD,INSERT_VALUES);
52  MatRestoreRow(K ,row,&ncols, (const int**)&cols, (const double**)&vals);
53  MatRestoreRow(K_transp,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
54  }
55  MatAssemblyBegin(D,MAT_FINAL_ASSEMBLY);
56  MatAssemblyEnd( D,MAT_FINAL_ASSEMBLY);
57 
58  /*Clean up and assign output pointer*/
59  MatFree(&K_transp);
60  *pD = D;
61 }/*}}}*/
62 void CreateLHS(Mat* pLHS,IssmDouble* pdmax,Mat K,Mat D,Vec Ml,IssmDouble theta,IssmDouble deltat,FemModel* femmodel,int configuration_type){/*{{{*/
63  /*Create Left Hand side of Lower order solution
64  *
65  * LHS = [ML - theta*detlat *(K+D)^n+1]
66  *
67  */
68 
69  /*Intermediaries*/
70  int dof,ncols,ncols2,rstart,rend;
71  double d,mi,dmax = 0.;
72  Mat LHS = NULL;
73  int* cols = NULL;
74  int* cols2 = NULL;
75  double* vals = NULL;
76  double* vals2 = NULL;
77 
78  MatDuplicate(K,MAT_SHARE_NONZERO_PATTERN,&LHS);
79  MatGetOwnershipRange(K,&rstart,&rend);
80  for(int row=rstart; row<rend; row++){
81  MatGetRow(K,row,&ncols, (const int**)&cols, (const double**)&vals);
82  MatGetRow(D,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
83  _assert_(ncols==ncols2);
84  for(int j=0; j<ncols; j++) {
85  _assert_(cols[j]==cols2[j]);
86  d = -theta*deltat*(vals[j] + vals2[j]);
87  if(cols[j]==row){
88  VecGetValues(Ml,1,(const int*)&cols[j],&mi);
89  d += mi;
90  }
91  if(fabs(d)>dmax) dmax = fabs(d);
92  MatSetValue(LHS,row,cols[j],d,INSERT_VALUES);
93  }
94  MatRestoreRow(K,row,&ncols, (const int**)&cols, (const double**)&vals);
95  MatRestoreRow(D,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
96  }
97  MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
98  MatAssemblyEnd( LHS,MAT_FINAL_ASSEMBLY);
99 
100  /*Deal with Dirichlet conditions: 2 options, penalties or zeros in K matrix*/
101  if(USEPENALTYMETHOD){
102  /*Option 1: Penalty method*/
103 
104  /*Broadcast max(dmax)*/
105  IssmDouble dmax_all;
106  ISSM_MPI_Reduce(&dmax,&dmax_all,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,0,IssmComm::GetComm() );
107  ISSM_MPI_Bcast(&dmax_all,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
108  dmax = dmax_all;
109 
110  dmax = dmax * 1.e+3;
111  for(int i=0;i<femmodel->constraints->Size();i++){
112  Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
113  constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
114  if(dof!=-1){
115  MatSetValue(LHS,dof,dof,dmax,INSERT_VALUES);
116  }
117  }
118  }
119  else{
120  /*Option 2: zero stiffness matrix, 1 one diagonal*/
121  int numrows = femmodel->constraints->Size();
122  int* rows = xNew<int>(numrows);
123  IssmDouble* rows_spc = xNew<IssmDouble>(numrows);
124  numrows = 0;
125 
126  dmax = dmax * 1.e+3;
127  for(int i=0;i<femmodel->constraints->Size();i++){
128  Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
129  constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
130  if(dof!=-1){
131  rows[numrows] = dof;
132  numrows++;
133  }
134  }
135  MatZeroRows(LHS,numrows,rows,1.,NULL,NULL);
136  xDelete<int>(rows);
137  xDelete<IssmDouble>(rows_spc);
138  }
139  MatAssemblyBegin(LHS,MAT_FINAL_ASSEMBLY);
140  MatAssemblyEnd( LHS,MAT_FINAL_ASSEMBLY);
141 
142  /*Clean up and assign output pointer*/
143  *pdmax = dmax;
144  *pLHS = LHS;
145 }/*}}}*/
146 void CreateRHS(Vec* pRHS,Mat K,Mat D,Vec Ml,Vec p,Vec u,IssmDouble theta,IssmDouble deltat,IssmDouble dmax,FemModel* femmodel,int configuration_type){/*{{{*/
147  /*Create Right Hand side of Lower order solution
148  *
149  * RHS = [ML + (1 - theta) deltaT L^n] u^n
150  *
151  * where L = K + D
152  *
153  */
154 
155  /*Intermediaries*/
156  Vec Ku = NULL;
157  Vec Du = NULL;
158  Vec RHS = NULL;
159  int dof;
160  IssmDouble d;
161 
162  /*Initialize vectors*/
163  VecDuplicate(u,&Ku);
164  VecDuplicate(u,&Du);
165  VecDuplicate(u,&RHS);
166 
167  /*Create RHS = M*u + (1-theta)*deltat*K*u + (1-theta)*deltat*D*u*/
168  MatMult(K,u,Ku);
169  MatMult(D,u,Du);
170  VecPointwiseMult(RHS,Ml,u);
171  VecAXPBYPCZ(RHS,(1-theta)*deltat,(1-theta)*deltat,1,Ku,Du);
172  VecAXPBY(RHS,deltat,1,p);
173  VecFree(&Ku);
174  VecFree(&Du);
175 
176  VecAssemblyBegin(RHS);
177  VecAssemblyEnd( RHS);
178 
179  /*Deal with Dirichlet conditions: 2 options, penalties or zeros in K matrix*/
180  if(USEPENALTYMETHOD){
181  /*Option 1: Penalty method*/
182  for(int i=0;i<femmodel->constraints->Size();i++){
183  Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
184  constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
185  d = d*dmax;
186  if(dof!=-1){
187  VecSetValues(RHS,1,&dof,(const double*)&d,INSERT_VALUES);
188  }
189  }
190  }
191  else{
192  /*Option 2: zero stiffness matrix, 1 one diagonal*/
193  int numrows = femmodel->constraints->Size();
194  int* rows = xNew<int>(numrows);
195  IssmDouble* rows_spc = xNew<IssmDouble>(numrows);
196  numrows = 0;
197  for(int i=0;i<femmodel->constraints->Size();i++){
198  Constraint* constraint=(Constraint*)femmodel->constraints->GetObjectByOffset(i);
199  constraint->PenaltyDofAndValue(&dof,&d,femmodel->nodes,femmodel->parameters);
200  if(dof!=-1){
201  rows[numrows] = dof;
202  rows_spc[numrows] = d;
203  numrows++;
204  }
205  }
206  VecSetValues(RHS,numrows,rows,rows_spc,INSERT_VALUES);
207  xDelete<int>(rows);
208  xDelete<IssmDouble>(rows_spc);
209  }
210  VecAssemblyBegin(RHS);
211  VecAssemblyEnd( RHS);
212 
213  /*Assign output pointer*/
214  *pRHS = RHS;
215 
216 }/*}}}*/
217 void RichardsonUdot(Vec* pudot,Vec u,Vec Ml,Mat K,Mat Mc){/*{{{*/
218  /*Use Richardson's formula to get udot using 5 steps and an initial guess of 0
219  *
220  * udot_new = udot_old + Ml^-1 (K^(n+1) u - Mc udot_old)
221  *
222  * */
223 
224  /*Intermediaries*/
225  Vec udot = NULL;
226  Vec temp1 = NULL;
227  Vec temp2 = NULL;
228 
229  /*Initialize vectors*/
230  VecDuplicate(u,&udot);
231  VecDuplicate(u,&temp1);
232  VecDuplicate(u,&temp2);
233 
234  /*Initial guess*/
235  VecZeroEntries(udot);
236 
237  /*Richardson's iterations*/
238  for(int i=0;i<5;i++){
239  MatMult(Mc,udot,temp1);
240  MatMult(K, u, temp2);
241  VecAXPBY(temp2,-1.,1.,temp1); // temp2 = (K^(n+1) u -- Mc udot_old)
242  VecPointwiseDivide(temp1,temp2,Ml); // temp1 = Ml^-1 temp2
243  VecAXPBY(udot,1.,1.,temp1);
244  }
245 
246  /*Clean up and assign output pointer*/
247  VecFree(&temp1);
248  VecFree(&temp2);
249  *pudot=udot;
250 
251 }/*}}}*/
252 void CreateRis(IssmDouble** pRi_plus_serial,IssmDouble** pRi_minus_serial,Mat Mc,Mat D,IssmDouble* ml_serial,Vec uvec,IssmDouble* u,IssmDouble* udot,IssmDouble* ulmin,IssmDouble* ulmax,IssmDouble deltat){/*{{{*/
253 
254  /*Intermediaries*/
255  Vec Ri_plus = NULL;
256  Vec Ri_minus = NULL;
257  int ncols,ncols2,rstart,rend;
258  double d;
259  int *cols = NULL;
260  int *cols2 = NULL;
261  double *vals = NULL;
262  double *vals2 = NULL;
263 
264  /*Initialize vectors*/
265  VecDuplicate(uvec,&Ri_plus);
266  VecDuplicate(uvec,&Ri_minus);
267 
268  /*Get global extremas*/
269  IssmDouble uLmin_global = ulmin[0];
270  IssmDouble uLmax_global = ulmax[0];
271  VecGetSize(uvec,&ncols);
272  for(int i=1;i<ncols;i++) if(ulmin[i]<uLmin_global) uLmin_global = ulmin[i];
273  for(int i=1;i<ncols;i++) if(ulmax[i]>uLmax_global) uLmax_global = ulmax[i];
274  //printf("%g %g\n",uLmin_global,uLmax_global);
275 
276  /*Go through D and calculate Ris*/
277  MatGetOwnershipRange(D,&rstart,&rend);
278  for(int row=rstart; row<rend; row++){
279  double Pi_plus = 0.;
280  double Pi_minus = 0.;
281  MatGetRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
282  MatGetRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
283  _assert_(ncols==ncols2);
284  for(int j=0; j<ncols; j++) {
285  _assert_(cols[j]==cols2[j]);
286  d = vals[j]*(udot[row] - udot[cols[j]]) + vals2[j]*(u[row] - u[cols[j]]);
287  if(row!=cols[j]){
288  if(d>0.){
289  Pi_plus += d;
290  }
291  else{
292  Pi_minus += d;
293  }
294  }
295  }
296 
297  /*Compute Qis and Ris*/
298  //double Qi_plus = ml_serial[row]/deltat*(3. - u[row]);
299  //double Qi_minus = ml_serial[row]/deltat*(2. - u[row]);
300  double Qi_plus = ml_serial[row]/deltat*(ulmax[row] - u[row]);
301  double Qi_minus = ml_serial[row]/deltat*(ulmin[row] - u[row]);
302  //double Qi_plus = ml_serial[row]/deltat*(uLmax_global - u[row]);
303  //double Qi_minus = ml_serial[row]/deltat*(uLmin_global - u[row]);
304  _assert_(Qi_plus >= 0.);
305  _assert_(Qi_minus <= 0.);
306  d = 1.;
307  if(Pi_plus!=0.) d = min(1.,Qi_plus/Pi_plus);
308  VecSetValue(Ri_plus,row,d,INSERT_VALUES);
309  d = 1.;
310  if(Pi_minus!=0.) d = min(1.,Qi_minus/Pi_minus);
311  VecSetValue(Ri_minus,row,d,INSERT_VALUES);
312 
313  MatRestoreRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
314  MatRestoreRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
315  }
316  VecAssemblyBegin(Ri_plus);
317  VecAssemblyEnd( Ri_plus);
318  VecAssemblyBegin(Ri_minus);
319  VecAssemblyEnd( Ri_minus);
320 
321  /*Serialize Ris*/
322  IssmDouble* Ri_plus_serial = NULL;
323  IssmDouble* Ri_minus_serial = NULL;
324  VecToMPISerial(&Ri_plus_serial,Ri_plus,IssmComm::GetComm());
325  VecToMPISerial(&Ri_minus_serial,Ri_minus,IssmComm::GetComm());
326  VecFree(&Ri_plus);
327  VecFree(&Ri_minus);
328 
329  /*Assign output pointer*/
330  *pRi_plus_serial = Ri_plus_serial;
331  *pRi_minus_serial = Ri_minus_serial;
332 }/*}}}*/
333 void CreateFbar(Vec* pFbar,IssmDouble* Ri_plus,IssmDouble* Ri_minus,Mat Mc,Mat D,IssmDouble* udot,IssmDouble* u,Vec uvec){/*{{{*/
334 
335  /*Intermediaries*/
336  Vec Fbar = NULL;
337  int ncols,ncols2,rstart,rend;
338  double d,f_ij;
339  int *cols = NULL;
340  int *cols2 = NULL;
341  double *vals = NULL;
342  double *vals2 = NULL;
343 
344  /*Build fbar*/
345  VecDuplicate(uvec,&Fbar);
346  MatGetOwnershipRange(D,&rstart,&rend);
347  for(int row=rstart; row<rend; row++){
348  MatGetRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
349  MatGetRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
350  _assert_(ncols==ncols2);
351  d = 0.;
352  for(int j=0; j<ncols; j++) {
353  _assert_(cols[j]==cols2[j]);
354  if(row==cols[j]) continue;
355  f_ij = vals[j]*(udot[row] - udot[cols[j]]) + vals2[j]*(u[row] - u[cols[j]]);
356  if(f_ij>0){
357  d += min(Ri_plus[row],Ri_minus[cols[j]]) * f_ij;
358  }
359  else{
360  d += min(Ri_minus[row],Ri_plus[cols[j]]) * f_ij;
361  }
362  }
363  VecSetValue(Fbar,row,d,INSERT_VALUES);
364  MatRestoreRow(Mc,row,&ncols, (const int**)&cols, (const double**)&vals);
365  MatRestoreRow(D ,row,&ncols2,(const int**)&cols2,(const double**)&vals2);
366  }
367  VecAssemblyBegin(Fbar);
368  VecAssemblyEnd( Fbar);
369 
370  /*Assign output pointer*/
371  *pFbar = Fbar;
372 }/*}}}*/
373 void UpdateSolution(Vec u,Vec udot,Vec Ml,Vec Fbar,IssmDouble deltat){/*{{{*/
374 
375  /*Intermediary*/
376  Vec temp1 = NULL;
377 
378  /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
379  VecDuplicate(u,&temp1);
380  VecPointwiseDivide(temp1,Fbar,Ml); //temp1 = Ml^-1 temp2
381  VecAXPBY(udot,1.,1.,temp1);
382  VecAXPY(u,deltat,temp1);
383 
384  /*CLean up and return*/
385  VecFree(&temp1);
386 }/*}}}*/
387 #endif
388 void solutionsequence_fct(FemModel* femmodel){/*{{{*/
389 
390  /*intermediary: */
391  IssmDouble theta,deltat,dmax;
392  int configuration_type,analysis_type;
393  Matrix<IssmDouble>* K = NULL;
394  Matrix<IssmDouble>* Mc = NULL;
395  Vector<IssmDouble>* p = NULL;
396  Vector<IssmDouble>* ug = NULL;
397  Vector<IssmDouble>* uf = NULL;
398  MasstransportAnalysis* manalysis = NULL;
399  DamageEvolutionAnalysis* danalysis = NULL;
400 
401  /*Recover parameters: */
402  femmodel->parameters->FindParam(&deltat,TimesteppingTimeStepEnum);
403  femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
404  femmodel->UpdateConstraintsx();
405  theta = 0.5; //0.5==Crank-Nicolson, 1==Backward Euler
406 
407  /*Create analysis*/
408  femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum);
409  switch(analysis_type){
410  case MasstransportAnalysisEnum:
411  manalysis = new MasstransportAnalysis();
412  manalysis->MassMatrix(&Mc,femmodel);
413  manalysis->FctKMatrix(&K,NULL,femmodel);
414  manalysis->FctPVector(&p,femmodel);
415  break;
416  case DamageEvolutionAnalysisEnum:
417  danalysis = new DamageEvolutionAnalysis();
418  danalysis->MassMatrix(&Mc,femmodel);
419  danalysis->FctKMatrix(&K,NULL,femmodel);
420  break;
421  default: _error_("analysis type " << EnumToStringx(analysis_type) << " not supported for FCT\n");
422  }
423  delete manalysis;
424  delete danalysis;
425 
426  #ifdef _HAVE_PETSC_
427 
428  /*Convert matrices to PETSc matrices*/
429  Mat D_petsc = NULL;
430  Mat LHS = NULL;
431  Vec RHS = NULL;
432  Vec u = NULL;
433  Vec udot = NULL;
434  Mat K_petsc = K->pmatrix->matrix;
435  Mat Mc_petsc = Mc->pmatrix->matrix;
436  Vec p_petsc = p->pvector->vector;
437 
438  /*Get previous solution u^n*/
439  GetSolutionFromInputsx(&ug,femmodel);
440  Reducevectorgtofx(&uf, ug, femmodel->nodes,femmodel->parameters);
441  delete ug;
442 
443  /*Compute lumped mass matrix*/
444  Vec Ml_petsc = NULL;
445  VecDuplicate(uf->pvector->vector,&Ml_petsc);
446  MatGetRowSum(Mc_petsc,Ml_petsc);
447 
448  /*Create D Matrix*/
449  CreateDMatrix(&D_petsc,K_petsc);
450 
451  /*Create LHS: [ML - theta*detlat *(K+D)^n+1]*/
452  CreateLHS(&LHS,&dmax,K_petsc,D_petsc,Ml_petsc,theta,deltat,femmodel,configuration_type);
453 
454  /*Create RHS: [ML + (1 - theta) deltaT L^n] u^n */
455  CreateRHS(&RHS,K_petsc,D_petsc,Ml_petsc,p_petsc,uf->pvector->vector,theta,deltat,dmax,femmodel,configuration_type);
456  delete uf;
457  delete p;
458 
459  /*Go solve lower order solution*/
460  femmodel->profiler->Start(SOLVER);
461  SolverxPetsc(&u,LHS,RHS,NULL,NULL, femmodel->parameters);
462  femmodel->profiler->Stop(SOLVER);
463  MatFree(&LHS);
464  VecFree(&RHS);
465 
466  /*Richardson to calculate udot*/
467  RichardsonUdot(&udot,u,Ml_petsc,K_petsc,Mc_petsc);
468  delete K;
469 
470  /*Serialize u and udot*/
471  IssmDouble* udot_serial = NULL;
472  IssmDouble* u_serial = NULL;
473  IssmDouble* ml_serial = NULL;
474  VecToMPISerial(&udot_serial,udot ,IssmComm::GetComm());
475  VecToMPISerial(&u_serial ,u ,IssmComm::GetComm());
476  VecToMPISerial(&ml_serial ,Ml_petsc,IssmComm::GetComm());
477 
478  /*Anti diffusive fluxes*/
479  Vec Fbar = NULL;
480  IssmDouble *Ri_minus_serial = NULL;
481  IssmDouble *Ri_plus_serial = NULL;
482  IssmDouble *ulmin = NULL;
483  IssmDouble *ulmax = NULL;
484  femmodel->GetInputLocalMinMaxOnNodesx(&ulmin,&ulmax,u_serial);
485  CreateRis(&Ri_plus_serial,&Ri_minus_serial,Mc_petsc,D_petsc,ml_serial,u,u_serial,udot_serial,ulmin,ulmax,deltat);
486  CreateFbar(&Fbar,Ri_plus_serial,Ri_minus_serial,Mc_petsc,D_petsc,udot_serial,u_serial,u);
487  xDelete<IssmDouble>(Ri_plus_serial);
488  xDelete<IssmDouble>(Ri_minus_serial);
489  xDelete<IssmDouble>(ulmin);
490  xDelete<IssmDouble>(ulmax);
491 
492  /*Clean up*/
493  MatFree(&D_petsc);
494  delete Mc;
495  xDelete<IssmDouble>(udot_serial);
496  xDelete<IssmDouble>(u_serial);
497  xDelete<IssmDouble>(ml_serial);
498 
499  /*Compute solution u^n+1 = u_L + deltat Ml^-1 fbar*/
500  UpdateSolution(u,udot,Ml_petsc,Fbar,deltat);
501  uf =new Vector<IssmDouble>(u);
502  VecFree(&u);
503  VecFree(&Fbar);
504  VecFree(&udot);
505  VecFree(&Ml_petsc);
506 
507  /*Update Element inputs*/
508  InputUpdateFromSolutionx(femmodel,uf);
509  delete uf;
510 
511  #else
512  _error_("PETSc needs to be installed");
513  #endif
514 }/*}}}*/
DataSet::Size
int Size()
Definition: DataSet.cpp:399
Matrix< IssmDouble >
SOLVER
#define SOLVER
Definition: Profiler.h:16
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
IssmDouble
double IssmDouble
Definition: types.h:37
MasstransportAnalysisEnum
@ MasstransportAnalysisEnum
Definition: EnumDefinitions.h:1163
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
DamageEvolutionAnalysisEnum
@ DamageEvolutionAnalysisEnum
Definition: EnumDefinitions.h:1026
ISSM_MPI_MAX
#define ISSM_MPI_MAX
Definition: issmmpi.h:131
VecFree
void VecFree(Vec *pvec)
Definition: VecFree.cpp:16
SolverxPetsc
void SolverxPetsc(Vec *puf, Mat Kff, Vec pf, Vec uf0, Vec df, Parameters *parameters)
Definition: PetscSolver.cpp:35
ISSM_MPI_DOUBLE
#define ISSM_MPI_DOUBLE
Definition: issmmpi.h:125
MasstransportAnalysis::MassMatrix
void MassMatrix(Matrix< IssmDouble > **pMff, FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:1176
FemModel::UpdateConstraintsx
void UpdateConstraintsx(void)
Definition: FemModel.cpp:3027
FemModel::nodes
Nodes * nodes
Definition: FemModel.h:56
VecToMPISerial
int VecToMPISerial(double **pgathered_vector, Vec vector, ISSM_MPI_Comm comm, bool broadcast=true)
Definition: VecToMPISerial.cpp:14
FemModel::GetInputLocalMinMaxOnNodesx
void GetInputLocalMinMaxOnNodesx(IssmDouble **pmin, IssmDouble **pmax, IssmDouble *ug)
Definition: FemModel.cpp:1308
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
ConfigurationTypeEnum
@ ConfigurationTypeEnum
Definition: EnumDefinitions.h:101
FemModel::constraints
Constraints * constraints
Definition: FemModel.h:52
DamageEvolutionAnalysis::FctKMatrix
void FctKMatrix(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, FemModel *femmodel)
Definition: DamageEvolutionAnalysis.cpp:868
Profiler::Stop
void Stop(int tagenum, bool dontmpisync=true)
Definition: Profiler.cpp:179
ISSM_MPI_Bcast
int ISSM_MPI_Bcast(void *buffer, int count, ISSM_MPI_Datatype datatype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:162
Reducevectorgtofx
void Reducevectorgtofx(Vector< IssmDouble > **puf, Vector< IssmDouble > *ug, Nodes *nodes, Parameters *parameters)
Definition: Reducevectorgtofx.cpp:8
FemModel
Definition: FemModel.h:31
MasstransportAnalysis::FctKMatrix
void FctKMatrix(Matrix< IssmDouble > **pKff, Matrix< IssmDouble > **pKfs, FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:1128
USEPENALTYMETHOD
#define USEPENALTYMETHOD
Definition: solutionsequence_fct.cpp:10
DamageEvolutionAnalysis
Definition: DamageEvolutionAnalysis.h:11
MatFree
void MatFree(Mat *pmat)
Definition: MatFree.cpp:16
Constraint::PenaltyDofAndValue
virtual void PenaltyDofAndValue(int *dof, IssmDouble *value, Nodes *nodes, Parameters *parameters)=0
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
Profiler::Start
void Start(int tagenum, bool dontmpisync=true)
Definition: Profiler.cpp:139
InputUpdateFromSolutionx
void InputUpdateFromSolutionx(FemModel *femmodel, Vector< IssmDouble > *solution)
Definition: InputUpdateFromSolutionx.cpp:9
GetSolutionFromInputsx
void GetSolutionFromInputsx(Vector< IssmDouble > **psolution, FemModel *femmodel)
Definition: GetSolutionFromInputsx.cpp:9
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
min
IssmDouble min(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:14
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
AnalysisTypeEnum
@ AnalysisTypeEnum
Definition: EnumDefinitions.h:36
Constraint
Definition: Constraint.h:17
DamageEvolutionAnalysis::MassMatrix
void MassMatrix(Matrix< IssmDouble > **pMff, FemModel *femmodel)
Definition: DamageEvolutionAnalysis.cpp:896
solutionsequence_fct
void solutionsequence_fct(FemModel *femmodel)
Definition: solutionsequence_fct.cpp:388
ISSM_MPI_Reduce
int ISSM_MPI_Reduce(void *sendbuf, void *recvbuf, int count, ISSM_MPI_Datatype datatype, ISSM_MPI_Op op, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:373
max
IssmDouble max(IssmDouble a, IssmDouble b)
Definition: extrema.cpp:24
MasstransportAnalysis::FctPVector
void FctPVector(Vector< IssmDouble > **ppf, FemModel *femmodel)
Definition: MasstransportAnalysis.cpp:1156
FemModel::profiler
Profiler * profiler
Definition: FemModel.h:42
Vector< IssmDouble >
MasstransportAnalysis
Definition: MasstransportAnalysis.h:11
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:22 for Ice Sheet System Model by   doxygen 1.8.19