Ice Sheet System Model  4.18
Code documentation
Functions
cores.h File Reference
#include "../shared/io/Comm/IssmComm.h"
#include "../shared/Numerics/types.h"

Go to the source code of this file.

Functions

void adjointstressbalance_core (FemModel *femmodel)
 
void adjointbalancethickness_core (FemModel *femmodel)
 
void adjointbalancethickness2_core (FemModel *femmodel)
 
void stressbalance_core (FemModel *femmodel)
 
void hydrology_core (FemModel *femmodel)
 
void thermal_core (FemModel *femmodel)
 
void surfaceslope_core (FemModel *femmodel)
 
void levelsetfunctionslope_core (FemModel *femmodel)
 
void movingfront_core (FemModel *femmodel)
 
void bedslope_core (FemModel *femmodel)
 
void meshdeformation_core (FemModel *femmodel)
 
void control_core (FemModel *femmodel)
 
void controltao_core (FemModel *femmodel)
 
void controlm1qn3_core (FemModel *femmodel)
 
void controladm1qn3_core (FemModel *femmodel)
 
void controlvalidation_core (FemModel *femmodel)
 
void masstransport_core (FemModel *femmodel)
 
void depthaverage_core (FemModel *femmodel)
 
void extrudefrombase_core (FemModel *femmodel)
 
void extrudefromtop_core (FemModel *femmodel)
 
void balancethickness_core (FemModel *femmodel)
 
void balancethickness2_core (FemModel *femmodel)
 
void balancevelocity_core (FemModel *femmodel)
 
void slopecompute_core (FemModel *femmodel)
 
void steadystate_core (FemModel *femmodel)
 
void transient_core (FemModel *femmodel)
 
void dakota_core (FemModel *femmodel)
 
void ad_core (FemModel *femmodel)
 
void adgradient_core (FemModel *femmodel)
 
void dummy_core (FemModel *femmodel)
 
void gia_core (FemModel *femmodel)
 
void love_core (FemModel *femmodel)
 
void esa_core (FemModel *femmodel)
 
void smb_core (FemModel *femmodel)
 
void bmb_core (FemModel *femmodel)
 
void damage_core (FemModel *femmodel)
 
void sealevelchange_core (FemModel *femmodel)
 
void grd_core (FemModel *femmodel)
 
void dynstr_core (FemModel *femmodel)
 
void sealevelrise_core_geometry (FemModel *femmodel)
 
SealevelMasks * sealevelrise_core_masks (FemModel *femmodel)
 
Vector< IssmDouble > * sealevelrise_core_eustatic (FemModel *femmodel, SealevelMasks *mask, IssmDouble *poceanarea)
 
Vector< IssmDouble > * sealevelrise_core_noneustatic (FemModel *femmodel, SealevelMasks *masks, Vector< IssmDouble > *RSLg_eustatic, IssmDouble oceanarea)
 
void sealevelrise_core_elastic (Vector< IssmDouble > **pU_radial, Vector< IssmDouble > **pU_north, Vector< IssmDouble > **pU_east, FemModel *femmodel, Vector< IssmDouble > *RSLg, SealevelMasks *masks)
 
void sealevelrise_core_viscous (Vector< IssmDouble > **pU_gia, Vector< IssmDouble > **pN_gia, FemModel *femmodel, Vector< IssmDouble > *RSLg)
 
void sealevelrise_diagnostics (FemModel *femmodel, Vector< IssmDouble > *RSLg)
 
IssmDouble objectivefunction (IssmDouble search_scalar, FemModel *femmodel)
 
void GetStericRate (Vector< IssmDouble > **psteric_rate_g, FemModel *femmodel)
 
void GetDynamicRate (Vector< IssmDouble > **pdynamic_rate_g, FemModel *femmodel)
 
int GradJSearch (IssmDouble *search_vector, FemModel *femmodel, int step)
 
void ProcessArguments (int *solution, char **pbinname, char **poutbinname, char **ptoolkitsname, char **plockname, char **prestartname, char **prootpath, int argc, char **argv)
 
void WriteLockFile (char *filename)
 
void ResetBoundaryConditions (FemModel *femmodel, int analysis_type)
 
void PrintBanner (void)
 
void TransferForcing (FemModel *femmodel, int forcingenum)
 
void TransferSealevel (FemModel *femmodel, int forcingenum)
 
void EarthMassTransport (FemModel *femmodel)
 
void slrconvergence (bool *pconverged, Vector< IssmDouble > *RSLg, Vector< IssmDouble > *RSLg_old, IssmDouble eps_rel, IssmDouble eps_abs)
 
void CorePointerFromSolutionEnum (void(**psolutioncore)(FemModel *), Parameters *parameters, int solutiontype)
 
void WrapperCorePointerFromSolutionEnum (void(**psolutioncore)(FemModel *), Parameters *parameters, int solutiontype, bool nodakotacore=false)
 
void AdjointCorePointerFromSolutionEnum (void(**padjointcore)(FemModel *), int solutiontype)
 

Function Documentation

◆ adjointstressbalance_core()

void adjointstressbalance_core ( FemModel *  femmodel)

Definition at line 12 of file adjointstressbalance_core.cpp.

12  {
13 
14  /*parameters: */
15  bool isFS,isSSA,isHO;
16  bool save_results;
17  bool conserve_loads = true;
18  int fe_FS;
19 
20  /*retrieve parameters:*/
21  femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
22  femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum);
23  femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum);
24  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
25  femmodel->parameters->FindParam(&fe_FS,FlowequationFeFSEnum);
26 
27  /*Compute velocities*/
28  if(VerboseSolution()) _printf0_(" computing velocities\n");
29  femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
30 
31  bool is_schur_cg_solver = false;
32  #ifdef _HAVE_PETSC_
33  int solver_type;
34  PetscOptionsDetermineSolverType(&solver_type);
35 
36  if(solver_type==FSSolverEnum) is_schur_cg_solver = true;
37  #endif
38 
39  if(isFS){
40  if (fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum)
41  solutionsequence_la(femmodel);
42  else if(is_schur_cg_solver)
43  solutionsequence_schurcg(femmodel);
44  else
45  solutionsequence_nonlinear(femmodel,conserve_loads);
46  }
47  else{
48  solutionsequence_nonlinear(femmodel,conserve_loads);
49  }
50 
51  /*Call SurfaceAreax, because some it might be needed by PVector*/
52  SurfaceAreax(NULL,femmodel);
53 
54  /*Compute adjoint*/
55  if(VerboseSolution()) _printf0_(" computing adjoint\n");
56  femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum,AdjointHorizAnalysisEnum);
57  solutionsequence_adjoint_linear(femmodel);
58 
59  /*Save results*/
60  if(save_results || true){
61  if(VerboseSolution()) _printf0_(" saving results\n");
62  if(isFS){
63  if(fe_FS==LATaylorHoodEnum || fe_FS==LACrouzeixRaviartEnum){
64  int outputs[2] = {AdjointxEnum,AdjointyEnum};
65  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
66  }
67  else{
68  if(!isSSA && !isHO){
69  int outputs[3] = {AdjointxEnum,AdjointyEnum,AdjointpEnum};
70  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
71  }
72  else{
73  int outputs[2] = {AdjointxEnum,AdjointyEnum};
74  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
75  }
76  }
77  }
78  else{
79  int outputs[2] = {AdjointxEnum,AdjointyEnum};
80  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
81  }
82  }
83 }

◆ adjointbalancethickness_core()

void adjointbalancethickness_core ( FemModel *  femmodel)

Definition at line 12 of file adjointbalancethickness_core.cpp.

12  {
13 
14  /*parameters: */
15  bool save_results;
16 
17  /*retrieve parameters:*/
18  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
19 
20  /*compute thickness */
21  if(VerboseSolution()) _printf0_(" computing thickness\n");
22  femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
23  solutionsequence_linear(femmodel);
24 
25  /*Call SurfaceAreax, because some it might be needed by PVector*/
26  SurfaceAreax(NULL,femmodel);
27 
28  /*compute adjoint*/
29  if(VerboseSolution()) _printf0_(" computing adjoint\n");
30  femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum,AdjointBalancethicknessAnalysisEnum);
31  solutionsequence_adjoint_linear(femmodel);
32 
33  /*Save results*/
34  if(save_results){
35  if(VerboseSolution()) _printf0_(" saving results\n");
36  int outputs[1] = {AdjointEnum};
37  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
38  }
39 }

◆ adjointbalancethickness2_core()

void adjointbalancethickness2_core ( FemModel *  femmodel)

Definition at line 12 of file adjointbalancethickness2_core.cpp.

12  {
13 
14  /*parameters: */
15  bool save_results;
16 
17  /*retrieve parameters:*/
18  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
19 
20  /*compute thickness2 */
21  if(VerboseSolution()) _printf0_(" computing thickness2\n");
22  femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
23  solutionsequence_linear(femmodel);
24 
25  /*Call SurfaceAreax, because some it might be needed by PVector*/
26  //SurfaceAreax(NULL,femmodel);
27 
28  /*compute adjoint*/
29  if(VerboseSolution()) _printf0_(" computing adjoint\n");
30  femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum,AdjointBalancethickness2AnalysisEnum);
31  solutionsequence_adjoint_linear(femmodel);
32 
33  /*Save results*/
34  if(save_results || true){
35  if(VerboseSolution()) _printf0_(" saving results\n");
36  int outputs[1] = {AdjointEnum};
37  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
38  }
39 }

◆ stressbalance_core()

void stressbalance_core ( FemModel *  femmodel)

Definition at line 13 of file stressbalance_core.cpp.

13  {
14 
15  /*Start profiler*/
16  femmodel->profiler->Start(STRESSBALANCECORE);
17 
18  /*parameters: */
19  bool dakota_analysis,control_analysis;
20  int domaintype;
21  bool isSIA,isSSA,isL1L2,isHO,isFS;
22  bool save_results;
23  int solution_type;
24  int numoutputs = 0;
25  char **requested_outputs = NULL;
26  Analysis *analysis = NULL;
27 
28  /* recover parameters:*/
29  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
30  femmodel->parameters->FindParam(&isSIA,FlowequationIsSIAEnum);
31  femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum);
32  femmodel->parameters->FindParam(&isL1L2,FlowequationIsL1L2Enum);
33  femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum);
34  femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
35  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
36  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
37  femmodel->parameters->FindParam(&numoutputs,StressbalanceNumRequestedOutputsEnum);
38  femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum);
39  if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,StressbalanceRequestedOutputsEnum);
40 
41  if(VerboseSolution()) _printf0_(" computing new velocity\n");
42  /*Compute slopes if necessary */
43  if(isSIA || (isFS && domaintype==Domain2DverticalEnum)) surfaceslope_core(femmodel);
44  if(isFS){
45  /*We need bed slopoes for the non-penetration condition*/
46  bedslope_core(femmodel);
47  femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
48  ResetFSBasalBoundaryConditionx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
49 
50  /*We need basal melt rates for the shelf dampening*/
51  bmb_core(femmodel);
52  }
53 
54  /*Compute SIA velocities*/
55  if(isSIA){
56 
57  /*Take the last velocity into account so that the velocity on the SSA domain is not zero*/
58  if(isSSA || isL1L2 || isHO ){
59  femmodel->SetCurrentConfiguration(StressbalanceSIAAnalysisEnum);
60  ResetBoundaryConditions(femmodel,StressbalanceSIAAnalysisEnum);
61  }
62 
63  analysis = new StressbalanceSIAAnalysis();
64  analysis->Core(femmodel);
65  delete analysis;
66 
67  /*Reset velocities for other ice flow models*/
68  if(isSSA || isL1L2 || isHO){
69  femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum);
70  ResetBoundaryConditions(femmodel,StressbalanceAnalysisEnum);
71  }
72  }
73 
74  /*Compute stressbalance for SSA L1L2 HO and FS*/
75  if(isSSA || isL1L2 || isHO || isFS){
76  analysis = new StressbalanceAnalysis();
77  analysis->Core(femmodel);
78  delete analysis;
79  }
80 
81  /*Compute vertical velocities*/
82  if (domaintype==Domain3DEnum && (isSIA || isSSA || isL1L2 || isHO)){
83 
84  /*We need basal melt rates for vertical velocity*/
85  bmb_core(femmodel);
86 
87  analysis = new StressbalanceVerticalAnalysis();
88  analysis->Core(femmodel);
89  delete analysis;
90  }
91 
92  if(save_results){
93  if(VerboseSolution()) _printf0_(" saving stressbalance results\n");
94  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
95  if(VerboseSolution()) _printf0_(" results saved\n");
96  }
97 
98  if(solution_type==StressbalanceSolutionEnum && !control_analysis)femmodel->RequestedDependentsx();
99 
100  /*Free ressources:*/
101  if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
102 
103  /*End profiler*/
104  femmodel->profiler->Stop(STRESSBALANCECORE);
105 }

◆ hydrology_core()

void hydrology_core ( FemModel *  femmodel)

Definition at line 12 of file hydrology_core.cpp.

12  {
13 
14  /*Start profiler*/
15  femmodel->profiler->Start(HYDROLOGYCORE);
16 
17  /*intermediary*/
18  int hydrology_model;
19  int solution_type;
20  int numoutputs = 0;
21  bool save_results;
22  bool modify_loads = true;
23  char **requested_outputs = NULL;
24  IssmDouble ThawedNodes;
25 
26  /*first recover parameters common to all solutions*/
27  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
28  femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
29  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
30  femmodel->parameters->FindParam(&numoutputs,HydrologyNumRequestedOutputsEnum);
31  if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,HydrologyRequestedOutputsEnum);
32 
33  /*Using the Shreve based Model*/
34  if (hydrology_model==HydrologyshreveEnum){
35  if(VerboseSolution()) _printf0_(" computing water heads\n");
36  /*first compute slopes: */
37  surfaceslope_core(femmodel);
38  bedslope_core(femmodel);
39  /*and then go to water column*/
40  if(VerboseSolution()) _printf0_(" computing water column\n");
41  femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum);
42  solutionsequence_nonlinear(femmodel,modify_loads);
43  /*transfer water column thickness to old water column thickness: */
44  InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
45 
46  }
47 
48  /*Using the double continuum model*/
49  else if (hydrology_model==HydrologydcEnum){
50  /*intermediary: */
51  bool isefficientlayer;
52  /*recover parameters: */
53  femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum);
54 
55  /*first we exclude frozen nodes of the solved nodes*/
56  femmodel->SetCurrentConfiguration(HydrologyDCInefficientAnalysisEnum);
57  femmodel->HydrologyIDSupdateDomainx(&ThawedNodes);
58 
59  if(ThawedNodes>0){
60  /*check if we need sub steps*/
61  int dtslices;
62  femmodel->parameters->FindParam(&dtslices,HydrologyStepsPerStepEnum);
63 
64  if(dtslices>1){
65  int substep, numaveragedinput, hydro_averaging;
66  IssmDouble global_time, subtime, yts;
67  IssmDouble dt, subdt;
68 
69  femmodel->parameters->FindParam(&global_time,TimeEnum);
70  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
71  femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
72  femmodel->parameters->FindParam(&hydro_averaging,HydrologyAveragingEnum);
73 
74  subtime=global_time-dt; //getting the time back to the start of the timestep
75  subdt=dt/dtslices; //computing hydro dt from dt and a divider
76  substep=0;
77  femmodel->parameters->SetParam(subdt,TimesteppingTimeStepEnum);
78  femmodel->parameters->SetParam(subtime,TimeEnum);
79 
80  /*intermiedaries to deal with averaging*/
81  static const int substeplist[4] = {EffectivePressureSubstepEnum,SedimentHeadSubstepEnum,EplHeadSubstepEnum,HydrologydcEplThicknessSubstepEnum};
82  static const int transientlist[4] = {EffectivePressureTransientEnum,SedimentHeadTransientEnum,EplHeadTransientEnum,HydrologydcEplThicknessTransientEnum};
83  static const int averagelist[4] = {EffectivePressureEnum,SedimentHeadEnum,EplHeadEnum,HydrologydcEplThicknessEnum};
84  std::vector<int> substepinput;
85  std::vector<int> transientinput;
86  std::vector<int> averagedinput;
87 
88  if (isefficientlayer){
89  /*define which variable needs to be averaged on the sub-timestep and initialize as needed*/
90  numaveragedinput = 4;
91  substepinput.assign(substeplist,substeplist+4);
92  transientinput.assign(transientlist,transientlist+4);
93  averagedinput.assign(averagelist,averagelist+4);
94  }
95  else{
96  numaveragedinput = 2;
97  substepinput.assign(substeplist,substeplist+2);
98  transientinput.assign(transientlist,transientlist+2);
99  averagedinput.assign(averagelist,averagelist+2);
100  }
101  femmodel->InitTransientInputx(&transientinput[0],numaveragedinput);
102  while(substep<dtslices){ //loop on hydro dts
103  substep+=1;
104  subtime+=subdt;
105  /*Setting substep time as global time*/
106  femmodel->parameters->SetParam(subtime,TimeEnum);
107  if(VerboseSolution()) _printf0_("sub iteration " << substep << "/" << dtslices << " time [yr]: " << setprecision(4) << subtime/yts << " (time step: " << subdt/yts << ")\n");
108  if(VerboseSolution()) _printf0_(" computing water heads\n");
109  /*save preceding timestep*/
110  InputDuplicatex(femmodel,SedimentHeadSubstepEnum,SedimentHeadOldEnum);
111  if (isefficientlayer){
112  InputDuplicatex(femmodel,EplHeadSubstepEnum,EplHeadOldEnum);
113  InputDuplicatex(femmodel,HydrologydcEplThicknessSubstepEnum,HydrologydcEplThicknessOldEnum);
114  }
115  /*Proceed now to heads computations*/
116  solutionsequence_hydro_nonlinear(femmodel);
117  /*If we have a sub-timestep we store the substep inputs in a transient input here*/
118  femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
119  }
120  /*averaging the stack*/
121  femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,hydro_averaging);
122  /*And reseting to global time*/
123  femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
124  femmodel->parameters->SetParam(global_time,TimeEnum);
125  }
126  else{
127  InputDuplicatex(femmodel,SedimentHeadSubstepEnum,SedimentHeadOldEnum);
128  if (isefficientlayer){
129  InputDuplicatex(femmodel,EplHeadSubstepEnum,EplHeadOldEnum);
130  InputDuplicatex(femmodel,HydrologydcEplThicknessSubstepEnum,HydrologydcEplThicknessOldEnum);
131  }
132  /*Proceed now to heads computations*/
133  if(VerboseSolution()) _printf0_(" computing water heads\n");
134  solutionsequence_hydro_nonlinear(femmodel);
135  /*If no substeps are present we want to duplicate the results for coupling purposes*/
136  InputDuplicatex(femmodel,SedimentHeadSubstepEnum,SedimentHeadEnum);
137  InputDuplicatex(femmodel,EffectivePressureSubstepEnum,EffectivePressureEnum);
138  if (isefficientlayer){
139  InputDuplicatex(femmodel,EplHeadSubstepEnum,EplHeadEnum);
140  InputDuplicatex(femmodel,HydrologydcEplThicknessSubstepEnum,HydrologydcEplThicknessEnum);
141  }
142  }
143  }
144  if(VerboseSolution())_printf0_(" hydroDC done\n");
145  }
146 
147  /*Using the SHAKTI model*/
148  else if (hydrology_model==HydrologyshaktiEnum){
149  femmodel->SetCurrentConfiguration(HydrologyShaktiAnalysisEnum);
150  InputDuplicatex(femmodel,HydrologyHeadEnum,HydrologyHeadOldEnum);
151  solutionsequence_shakti_nonlinear(femmodel);
152  if(VerboseSolution()) _printf0_(" updating gap height\n");
153  HydrologyShaktiAnalysis* analysis = new HydrologyShaktiAnalysis();
154  analysis->UpdateGapHeight(femmodel);
155  delete analysis;
156  }
157 
158  /*Using the GlaDS model*/
159  else if (hydrology_model==HydrologyGlaDSEnum){
160  HydrologyGlaDSAnalysis* analysis = new HydrologyGlaDSAnalysis();
161  femmodel->SetCurrentConfiguration(HydrologyGlaDSAnalysisEnum);
162 
163  /*Set fields as old*/
164  InputDuplicatex(femmodel,HydraulicPotentialEnum,HydraulicPotentialOldEnum);
165  InputDuplicatex(femmodel,HydrologySheetThicknessEnum,HydrologySheetThicknessOldEnum);
166  analysis->SetChannelCrossSectionOld(femmodel);
167 
168  /*Solve for new potential*/
169  solutionsequence_glads_nonlinear(femmodel);
170 
171  if(VerboseSolution()) _printf0_(" updating effective pressure\n");
172  analysis->UpdateEffectivePressure(femmodel);
173  delete analysis;
174  }
175 
176  /*Using the PISM hydrology model*/
177  else if (hydrology_model==HydrologypismEnum){
178  femmodel->SetCurrentConfiguration(HydrologyPismAnalysisEnum);
179  if(VerboseSolution()) _printf0_(" updating water column\n");
180  HydrologyPismAnalysis* analysis = new HydrologyPismAnalysis();
181  InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum);
182  analysis->UpdateWaterColumn(femmodel);
183  delete analysis;
184  }
185  else{
186  _error_("Hydrology model "<< EnumToStringx(hydrology_model) <<" not supported yet");
187  }
188  if(save_results){
189  if(VerboseSolution()) _printf0_(" saving hydrology results \n");
190  if(hydrology_model==HydrologydcEnum && ThawedNodes==0){
191  if(VerboseSolution()) _printf0_(" No thawed node hydro is skiped \n");}
192  else{
193  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
194  }
195  }
196  /*Free ressources:*/
197  if(numoutputs){
198  for(int i=0;i<numoutputs;i++){
199  xDelete<char>(requested_outputs[i]);
200  }
201  xDelete<char*>(requested_outputs);
202  }
203  /*End profiler*/
204  femmodel->profiler->Stop(HYDROLOGYCORE);
205 }

◆ thermal_core()

void thermal_core ( FemModel *  femmodel)

Definition at line 13 of file thermal_core.cpp.

13  {
14 
15  /*Start profiler*/
16  femmodel->profiler->Start(THERMALCORE);
17 
18  /*intermediary*/
19  bool save_results,isenthalpy;
20  bool dakota_analysis;
21  int solution_type,numoutputs;
22  char** requested_outputs = NULL;
23  EnthalpyAnalysis * enthalpy_analysis = NULL;
24 
25  /*first recover parameters common to all solutions*/
26  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
27  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
28  femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
29  femmodel->parameters->FindParam(&numoutputs,ThermalNumRequestedOutputsEnum);
30  if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,ThermalRequestedOutputsEnum);
31 
32  /*Calculate geothermalflux*/
33  GeothermalFluxx(femmodel);
34 
35  if(isenthalpy){
36  femmodel->InputMakeDiscontinuous(BasalforcingsGroundediceMeltingRateEnum);
37  enthalpy_analysis = new EnthalpyAnalysis();
38  enthalpy_analysis->Core(femmodel);
39  delete enthalpy_analysis;
40  }
41  else{
42  if(VerboseSolution()) _printf0_(" computing temperatures\n");
43  femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
44  solutionsequence_thermal_nonlinear(femmodel);
45 
46  if(VerboseSolution()) _printf0_(" computing melting\n");
47  femmodel->SetCurrentConfiguration(MeltingAnalysisEnum);
48  solutionsequence_linear(femmodel);
49  }
50 
51  if(save_results){
52  if(VerboseSolution()) _printf0_(" saving thermal results\n");
53  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
54  }
55 
56  /*Free ressources:*/
57  if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
58 
59  /*End profiler*/
60  femmodel->profiler->Stop(THERMALCORE);
61 }

◆ surfaceslope_core()

void surfaceslope_core ( FemModel *  femmodel)

Definition at line 12 of file surfaceslope_core.cpp.

12  {
13 
14  /*parameters: */
15  bool save_results;
16  int domaintype;
17 
18  /*Recover some parameters: */
19  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
20  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
21 
22  if(VerboseSolution()) _printf0_("computing slope...\n");
23 
24  /*Call on core computations: */
25  femmodel->SetCurrentConfiguration(L2ProjectionBaseAnalysisEnum);
26 
27  femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToL2ProjectEnum);
28  solutionsequence_linear(femmodel);
29 
30  if(domaintype!=Domain2DverticalEnum){
31  femmodel->parameters->SetParam(SurfaceSlopeYEnum,InputToL2ProjectEnum);
32  solutionsequence_linear(femmodel);
33  }
34  if(domaintype==Domain2DverticalEnum){
35  femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToExtrudeEnum);
36  extrudefrombase_core(femmodel);
37  }
38 
39  if(save_results){
40  if(VerboseSolution()) _printf0_("saving surface slopes results:\n");
41  if(domaintype!=Domain2DverticalEnum){
42  int outputs[2] = {SurfaceSlopeXEnum,SurfaceSlopeYEnum};
43  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
44 
45  }
46  else{
47  int outputs = SurfaceSlopeXEnum;
48  femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
49  }
50  }
51 
52 }

◆ levelsetfunctionslope_core()

void levelsetfunctionslope_core ( FemModel *  femmodel)

Definition at line 12 of file levelsetfunctionslope_core.cpp.

12  {
13 
14  /*parameters: */
15  bool save_results;
16  int domaintype;
17 
18  /*Recover some parameters: */
19  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
20  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
21 
22  if(VerboseSolution()) _printf0_(" computing slope of levelset function...\n");
23 
24  /*Call on core computations: */
25  femmodel->SetCurrentConfiguration(L2ProjectionBaseAnalysisEnum);
26 
27  femmodel->parameters->SetParam(LevelsetfunctionSlopeXEnum,InputToL2ProjectEnum);
28  solutionsequence_linear(femmodel);
29 
30  if(domaintype!=Domain2DverticalEnum){
31  femmodel->parameters->SetParam(LevelsetfunctionSlopeYEnum,InputToL2ProjectEnum);
32  solutionsequence_linear(femmodel);
33  }
34  if(domaintype==Domain2DverticalEnum){
35  femmodel->parameters->SetParam(LevelsetfunctionSlopeXEnum,InputToExtrudeEnum);
36  extrudefrombase_core(femmodel);
37  }
38 }

◆ movingfront_core()

void movingfront_core ( FemModel *  femmodel)

Definition at line 12 of file movingfront_core.cpp.

12  {
13 
14  /*Start profiler*/
15  femmodel->profiler->Start(MOVINGFRONTCORE);
16 
17  /* intermediaries */
18  bool save_results,isstressbalance,ismasstransport,isthermal,isenthalpy,islevelset,ismovingfront,killicebergs;
19  int domaintype, num_extrapol_vars, index,reinit_frequency,step;
20  int* extrapol_vars=NULL;
21  Analysis *analysis=NULL;
22 
23  /* recover parameters */
24  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
25  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
26  femmodel->parameters->FindParam(&isstressbalance,TransientIsstressbalanceEnum);
27  femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
28  femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
29  femmodel->parameters->FindParam(&ismovingfront,TransientIsmovingfrontEnum);
30  femmodel->parameters->FindParam(&reinit_frequency,LevelsetReinitFrequencyEnum);
31  femmodel->parameters->FindParam(&killicebergs,LevelsetKillIcebergsEnum);
32  femmodel->parameters->FindParam(&step,StepEnum);
33  if(isthermal && domaintype==Domain3DEnum) femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
34 
35  if(!ismovingfront) return;
36 
37  /* start the work from here */
38  Calvingx(femmodel);
39  FrontalForcingsx(femmodel);
40  if(VerboseSolution()) _printf0_(" computing level set transport\n");
41 
42  /* smoothen slope of lsf for computation of normal on ice domain*/
43  levelsetfunctionslope_core(femmodel);
44 
45  /* determine variables for extrapolation */
46  num_extrapol_vars=0;
47  if(isstressbalance){
48  if(domaintype==Domain3DEnum)
49  num_extrapol_vars+=3;
50  else
51  num_extrapol_vars+=2;
52  }
53  if(ismasstransport) num_extrapol_vars+=1;
54  if(isthermal && domaintype==Domain3DEnum) num_extrapol_vars+=1;
55  extrapol_vars=xNew<int>(num_extrapol_vars);
56  index=0;
57  if(isstressbalance){
58  extrapol_vars[index]=VxEnum; index++;
59  extrapol_vars[index]=VyEnum; index++;
60  if(domaintype==Domain3DEnum){
61  extrapol_vars[index]=VzEnum; index++;
62  }
63  }
64  if(ismasstransport){
65  extrapol_vars[index]=ThicknessEnum; index++;
66  }
67  if(isthermal && domaintype==Domain3DEnum){
68  if(isenthalpy){
69  extrapol_vars[index]=EnthalpyEnum;
70  }
71  else{
72  extrapol_vars[index]=TemperatureEnum;
73  }
74  index++;
75  }
76 
77  /* extrapolate */
78  analysis = new ExtrapolationAnalysis();
79  for(int iv=0;iv<num_extrapol_vars;iv++){
80  femmodel->parameters->SetParam(extrapol_vars[iv],ExtrapolationVariableEnum);
81  analysis->Core(femmodel);
82  }
83  xDelete<int>(extrapol_vars);
84  delete analysis;
85 
86  /* solve level set equation */
87  analysis = new LevelsetAnalysis();
88  analysis->Core(femmodel);
89  delete analysis;
90 
91  /*Kill ice berg to avoid free body motion*/
92  if(killicebergs){
93  if(VerboseSolution()) _printf0_(" looking for icebergs to kill\n");
94  KillIcebergsx(femmodel);
95  }
96 
97  /*Reset levelset if needed*/
98  if(reinit_frequency && (step%reinit_frequency==0)){
99  if(VerboseSolution()) _printf0_(" reinitializing level set\n");
100  femmodel->ResetLevelset();
101  }
102 
103  /* update vertices included for next calculation */
104  GetMaskOfIceVerticesLSMx(femmodel);
105 
106  /* add computation domain mask to outputs */
107  if(save_results){
108  int outputs[1] = {IceMaskNodeActivationEnum};
109  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
110  }
111 
112  /*End profiler*/
113  femmodel->profiler->Stop(MOVINGFRONTCORE);
114 }

◆ bedslope_core()

void bedslope_core ( FemModel *  femmodel)

Definition at line 12 of file bedslope_core.cpp.

12  {
13 
14  /*parameters: */
15  bool save_results;
16  int domaintype;
17 
18  /*Recover some parameters: */
19  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
20  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
21 
22  if(VerboseSolution()) _printf0_(" computing slope\n");
23 
24  /*Call on core computations: */
25  femmodel->SetCurrentConfiguration(L2ProjectionBaseAnalysisEnum);
26 
27  femmodel->parameters->SetParam(BedSlopeXEnum,InputToL2ProjectEnum);
28  solutionsequence_linear(femmodel);
29 
30  if(domaintype!=Domain2DverticalEnum){
31  femmodel->parameters->SetParam(BedSlopeYEnum,InputToL2ProjectEnum);
32  solutionsequence_linear(femmodel);
33  }
34 
35  if(save_results){
36  if(VerboseSolution()) _printf0_(" saving bedslopes results\n");
37  if(domaintype!=Domain2DverticalEnum){
38  int outputs[2] = {BedSlopeXEnum,BedSlopeYEnum};
39  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
40  }
41  else{
42  int outputs[1] = {BedSlopeXEnum};
43  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1);
44  }
45  }
46 
47 }

◆ meshdeformation_core()

void meshdeformation_core ( FemModel *  femmodel)

◆ control_core()

void control_core ( FemModel *  femmodel)

Definition at line 22 of file control_core.cpp.

22  {/*{{{*/
23 
24  /*parameters: */
25  int num_controls,nsize,nsteps;
26  int solution_type;
27  bool isFS,dakota_analysis;
28  int *control_type = NULL;
29  int *maxiter = NULL;
30  IssmDouble *cm_jump = NULL;
31  IssmDouble *J = NULL;
32 
33  /*Solution and Adjoint core pointer*/
34  void (*solutioncore)(FemModel*) = NULL;
35  void (*adjointcore)(FemModel*) = NULL;
36 
37  /*Recover parameters used throughout the solution*/
38  femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum);
39  femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum);
40  femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum);
41  femmodel->parameters->FindParam(&maxiter,NULL,InversionMaxiterPerStepEnum);
42  femmodel->parameters->FindParam(&cm_jump,NULL,InversionStepThresholdEnum);
43  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
44  femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
45  femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
46  femmodel->parameters->SetParam(false,SaveResultsEnum);
47 
48  /*out of solution_type, figure out solution core and adjoint function pointer*/
49  CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
50  AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
51 
52  /*Launch once a complete solution to set up all inputs*/
53  if(VerboseControl()) _printf0_(" preparing initial solution\n");
54  if(isFS) solutioncore(femmodel);
55 
56  /*Get initial guess*/
57  Vector<IssmDouble> *Xpetsc = NULL;
58  GetVectorFromControlInputsx(&Xpetsc,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
59  IssmDouble* X0 = Xpetsc->ToMPISerial();
60  Xpetsc->GetSize(&nsize);
61  delete Xpetsc;
62 
63  /*Initialize some of the BrentSearch arguments: */
64  OptPars optpars;
65  optpars.xmin = 0;
66  optpars.xmax = 1;
67  optpars.nsteps = nsteps;
68  optpars.nsize = nsize;
69  optpars.maxiter = maxiter;
70  optpars.cm_jump = cm_jump;
71 
72  /*Initialize function argument*/
73  AppCtx usr;
74  usr.femmodel = femmodel;
75  usr.nsize = nsize;
76 
77  /*Call Brent optimization*/
78  BrentSearch(&J,optpars,X0,&FormFunction,&FormFunctionGradient,(void*)&usr);
79 
80  if(VerboseControl()) _printf0_(" preparing final solution\n");
81  IssmDouble *XL = NULL;
82  IssmDouble *XU = NULL;
83  GetVectorFromControlInputsx(&XL,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound");
84  GetVectorFromControlInputsx(&XU,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound");
85  for(long i=0;i<nsize;i++){
86  if(X0[i]>XU[i]) X0[i]=XU[i];
87  if(X0[i]<XL[i]) X0[i]=XL[i];
88  }
89  xDelete<IssmDouble>(XU);
90  xDelete<IssmDouble>(XL);
91  SetControlInputsFromVectorx(femmodel,X0);
92  femmodel->parameters->SetParam(true,SaveResultsEnum);
93  solutioncore(femmodel);
94 
95  /*some results not computed by steadystate_core or stressbalance_core: */
96  if(!dakota_analysis){ //do not save this if we are running the control core from a qmu run!
97  femmodel->OutputControlsx(&femmodel->results);
98 
99  #ifdef _HAVE_AD_
100  IssmPDouble* J_passive=xNew<IssmPDouble>(nsteps);
101  for(int i=0;i<nsteps;i++) J_passive[i]=reCast<IssmPDouble>(J[i]);
102  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,nsteps,1,0,0));
103  xDelete<IssmPDouble>(J_passive);
104  #else
105  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J,nsteps,1,0,0));
106  #endif
107  }
108 
109  /*Free ressources: */
110  xDelete<int>(control_type);
111  xDelete<int>(maxiter);
112  xDelete<IssmDouble>(cm_jump);
113  xDelete<IssmDouble>(J);
114  xDelete<IssmDouble>(X0);
115 }/*}}}*/

◆ controltao_core()

void controltao_core ( FemModel *  femmodel)

Definition at line 213 of file controltao_core.cpp.

213  {
214  _error_("TAO not installed or PETSc version not supported");
215 }

◆ controlm1qn3_core()

void controlm1qn3_core ( FemModel *  femmodel)

Definition at line 304 of file controlm1qn3_core.cpp.

304 {_error_("M1QN3 not installed");}

◆ controladm1qn3_core()

void controladm1qn3_core ( FemModel *  femmodel)

Definition at line 657 of file controladm1qn3_core.cpp.

657 {_error_("M1QN3 or ADOLC/CoDiPack not installed");}

◆ controlvalidation_core()

void controlvalidation_core ( FemModel *  femmodel)

Definition at line 155 of file controlvalidation_core.cpp.

155  {
156 
157  int solution_type,n;
158  int num_responses;
159  IssmDouble j0,j;
160  IssmDouble Ialpha,exponent,alpha;
161  IssmDouble* scaling_factors = NULL;
162  IssmDouble* jlist = NULL;
163  int my_rank=IssmComm::GetRank();
164 
165  /*Recover parameters used throughout the solution*/
166  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
167  femmodel->parameters->SetParam(false,SaveResultsEnum);
168  femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum);
169  femmodel->parameters->FindParam(&scaling_factors,NULL,InversionControlScalingFactorsEnum);
170 
171  /*Get initial guess*/
172  IssmPDouble* X0 = NULL;
173  GetPassiveVectorFromControlInputsx(&X0,&n,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value");
174 
175  /*Allocate vectors*/
176  IssmDouble* X = xNew<IssmDouble>(n);
177  IssmPDouble* G = xNew<IssmPDouble>(n);
178 
179  /*out of solution_type, figure out solution core and adjoint function pointer*/
180  void (*solutioncore)(FemModel*)=NULL;
181  CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type);
182 
183  #if defined(_HAVE_ADOLC_)
184  /*{{{*/
185  IssmDouble* aX=xNew<IssmDouble>(n);
186  if(my_rank==0){
187  for(int i=0;i<n;i++){
188  aX[i]<<=X0[i];
189  }
190  }
191  _error_("not implemented yet...");
192  /*}}}*/
193  #elif defined(_HAVE_CODIPACK_)
194  simul_starttrace2(femmodel);
195  IssmDouble* aX=xNew<IssmDouble>(n);
196  auto& tape_codi = IssmDouble::getGlobalTape();
197  codi_global.input_indices.clear();
198  if(my_rank==0){
199  for (int i=0;i<n;i++) {
200  aX[i]=X0[i];
201  tape_codi.registerInput(aX[i]);
202  codi_global.input_indices.push_back(aX[i].getGradientData());
203  }
204  }
205  SetControlInputsFromVectorx(femmodel,aX);
206  xDelete(aX);
207 
208  if(VerboseControl()) _printf0_(" Compute Initial cost function\n");
209  solutioncore(femmodel);
210 
211  /*Get Dependents*/
212  IssmDouble output_value;
213  int num_dependents;
214  IssmPDouble *dependents;
215  DataSet* dependent_objects=NULL;
216  IssmDouble J=0.;
217  femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
218  femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
219 
220  /*Go through our dependent variables, and compute the response:*/
221  dependents=xNew<IssmPDouble>(num_dependents);
222  codi_global.output_indices.clear();
223  for(int i=0;i<dependent_objects->Size();i++){
224  DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
225  if(solution_type==TransientSolutionEnum){
226  output_value = dep->GetValue();
227  }
228  else{
229  dep->Responsex(&output_value,femmodel);
230  }
231  _printf0_("=== output ="<<output_value<<" \n");
232  if(my_rank==0) {
233  tape_codi.registerOutput(output_value);
234  dependents[i] = output_value.getValue();
235  codi_global.output_indices.push_back(output_value.getGradientData());
236  J+=output_value;
237  }
238  }
239  j0 = J;
240  _printf0_("Initial cost function J(x) = "<<setw(12)<<setprecision(7)<<j0<<"\n");
241  _assert_(j0>0.);
242  simul_stoptrace2();
243  /*initialize direction index in the weights vector: */
244  if(my_rank==0){
245  tape_codi.setGradient(codi_global.output_indices[0],1.0);
246  }
247  tape_codi.evaluate();
248 
249  /*Get gradient for this dependent */
250  auto in_size = codi_global.input_indices.size();
251  for(size_t i = 0; i < in_size; ++i) {
252  G[i] = tape_codi.getGradient(codi_global.input_indices[i]);
253  }
254  #else
255  /*{{{*/
256  void (*adjointcore)(FemModel*) = NULL;
257  AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type);
258 
259  if(VerboseControl()) _printf0_(" Compute Initial solution\n");
260  solutioncore(femmodel);
261  if(VerboseControl()) _printf0_(" Compute Adjoint\n");
262  adjointcore(femmodel);
263 
264  if(VerboseControl()) _printf0_(" Compute Initial cost function\n");
265  femmodel->CostFunctionx(&j0,&jlist,NULL);
266  _printf0_("Initial cost function J(x) = "<<setw(12)<<setprecision(7)<<j0<<"\n");
267  xDelete<IssmDouble>(jlist);
268 
269  if(VerboseControl()) _printf0_(" Compute Gradient\n");
270  Gradjx(&G,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
271  for(int i=0;i<n;i++) G[i] = -G[i];
272  /*}}}*/
273  #endif
274 
275  /*Allocate output*/
276  int num = 26;
277  IssmDouble* output = xNew<IssmDouble>(2*num);
278 
279  /*Start loop*/
280  _printf0_(" alpha Ialpha \n");
281  _printf0_("_________________________\n");
282  for(int m=0;m<num;m++){
283 
284  /*Create new vector*/
285  alpha = pow(2.,-m);
286  for(int i=0;i<n;i++) X[i] = X0[i] + alpha*scaling_factors[0];
287 
288  /*Calculate j(k+alpha delta k) */
289  SetControlInputsFromVectorx(femmodel,X);
290  solutioncore(femmodel);
291 
292  #if defined(_HAVE_CODIPACK_)
293  j=0.;
294  for(int i=0;i<dependent_objects->Size();i++){
295  DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
296  if(solution_type==TransientSolutionEnum){
297  output_value = dep->GetValue();
298  }
299  else{
300  dep->Responsex(&output_value,femmodel);
301  }
302  j+=output_value;
303  }
304  #else
305  femmodel->CostFunctionx(&j,NULL,NULL);
306  #endif
307 
308  IssmDouble Den = 0.;
309  for(int i=0;i<n;i++) Den += alpha* G[i] * scaling_factors[0];
310  Ialpha = fabs((j - j0)/Den - 1.);
311  _assert_(fabs(Den)>0.);
312 
313  _printf0_(" " << setw(11) << setprecision (5)<<alpha<<" " << setw(11) << setprecision (5)<<Ialpha<<"\n");
314  output[m*2+0] = alpha;
315  output[m*2+1] = Ialpha;
316  }
317 
318  /*output*/
319  #ifdef _HAVE_AD_
320  IssmPDouble* J_passive=xNew<IssmPDouble>(2*num);
321  for(int i=0;i<2*num;i++) J_passive[i]=reCast<IssmPDouble>(output[i]);
322  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,J_passive,num,2,0,0));
323  xDelete<IssmPDouble>(J_passive);
324  IssmDouble* aG=xNew<IssmDouble>(n);
325  for(int i=0;i<n;i++) aG[i] = G[i];
326  ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,aG);
327  xDelete<IssmDouble>(aG);
328  #else
329  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,JEnum,output,num,2,0,0));
330  ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G);
331  #endif
332  femmodel->OutputControlsx(&femmodel->results);
333 
334  /*Clean up and return*/
335  xDelete<IssmDouble>(output);
336  xDelete<IssmPDouble>(G);
337  xDelete<IssmDouble>(X);
338  xDelete<double>(X0);
339  xDelete<IssmDouble>(scaling_factors);
340 }

◆ masstransport_core()

void masstransport_core ( FemModel *  femmodel)

Definition at line 12 of file masstransport_core.cpp.

12  {
13 
14  /*Start profiler*/
15  femmodel->profiler->Start(MASSTRANSPORTCORE);
16 
17  /*parameters: */
18  int numoutputs,domaintype;
19  bool save_results;
20  bool isFS,isfreesurface,dakota_analysis;
21  int solution_type,stabilization;
22  char** requested_outputs = NULL;
23 
24  /*activate configuration*/
25  femmodel->SetCurrentConfiguration(MasstransportAnalysisEnum);
26 
27  /*recover parameters: */
28  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
29  femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
30  femmodel->parameters->FindParam(&isfreesurface,MasstransportIsfreesurfaceEnum);
31  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
32  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
33  femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum);
34  femmodel->parameters->FindParam(&stabilization,MasstransportStabilizationEnum);
35  if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum);
36 
37  if(VerboseSolution()) _printf0_(" computing mass transport\n");
38 
39  /*Transport mass or free surface*/
40  if(isFS && isfreesurface){
41  if(VerboseSolution()) _printf0_(" call free surface computational core\n");
42  femmodel->SetCurrentConfiguration(FreeSurfaceBaseAnalysisEnum);
43  solutionsequence_linear(femmodel);
44  femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
45  extrudefrombase_core(femmodel);
46 
47  femmodel->SetCurrentConfiguration(FreeSurfaceTopAnalysisEnum);
48  solutionsequence_linear(femmodel);
49  femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
50  extrudefromtop_core(femmodel);
51  }
52  else{
53  if(VerboseSolution()) _printf0_(" call computational core\n");
54  femmodel->parameters->SetParam(VxEnum,InputToDepthaverageInEnum);
55  femmodel->parameters->SetParam(VxAverageEnum,InputToDepthaverageOutEnum);
56  depthaverage_core(femmodel);
57  if(domaintype!=Domain2DverticalEnum){
58  femmodel->parameters->SetParam(VyEnum,InputToDepthaverageInEnum);
59  femmodel->parameters->SetParam(VyAverageEnum,InputToDepthaverageOutEnum);
60  depthaverage_core(femmodel);
61  }
62  femmodel->SetCurrentConfiguration(MasstransportAnalysisEnum);
63  InputDuplicatex(femmodel,ThicknessEnum,ThicknessOldEnum);
64  InputDuplicatex(femmodel,BaseEnum,BaseOldEnum);
65  InputDuplicatex(femmodel,SurfaceEnum,SurfaceOldEnum);
66  InputDuplicatex(femmodel,SealevelriseCumDeltathicknessEnum,SealevelriseCumDeltathicknessOldEnum);
67  if(stabilization==4){
68  solutionsequence_fct(femmodel);
69  }
70  else{
71  solutionsequence_linear(femmodel);
72  // ThicknessAverage: method not totally tested
73  //if(stabilization==3){
74  // if(VerboseSolution()) _printf0_(" call thickness average core\n");
75  // femmodel->ThicknessAverage();
76  //}
77  }
78  femmodel->parameters->SetParam(ThicknessEnum,InputToExtrudeEnum);
79  extrudefrombase_core(femmodel);
80  femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
81  extrudefrombase_core(femmodel);
82  femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
83  extrudefrombase_core(femmodel);
84  }
85 
86  if(save_results){
87  if(VerboseSolution()) _printf0_(" saving results\n");
88  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
89  }
90 
91  if(solution_type==MasstransportSolutionEnum)femmodel->RequestedDependentsx();
92 
93  /*Free ressources:*/
94  if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
95 
96  /*profiler*/
97  femmodel->profiler->Stop(MASSTRANSPORTCORE);
98 }

◆ depthaverage_core()

void depthaverage_core ( FemModel *  femmodel)

Definition at line 12 of file depthaverage_core.cpp.

12  {
13 
14  /*Intermediaries*/
15  int domaintype,elementtype;
16  int inputenum,input_average_enum;
17 
18  if(VerboseSolution()) _printf0_(" depth averaging solution...\n");
19 
20  /*Get parameters*/
21  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
22  femmodel->parameters->FindParam(&elementtype,MeshElementtypeEnum);
23  femmodel->parameters->FindParam(&inputenum,InputToDepthaverageInEnum);
24  femmodel->parameters->FindParam(&input_average_enum,InputToDepthaverageOutEnum);
25 
26  /*If this is a 2D horizontal domain: no need to do anything, just copy input*/
27  if(domaintype==Domain2DhorizontalEnum){
28  InputDuplicatex(femmodel,inputenum,input_average_enum);
29  return;
30  }
31 
32  /*Special method for Penta, otherwise call solution sequence*/
33  if(elementtype==PentaEnum){
34  InputDepthAverageAtBasex(femmodel,inputenum,input_average_enum);
35  }
36  else{
37  /*Call on core computations: */
38  femmodel->SetCurrentConfiguration(DepthAverageAnalysisEnum);
39  solutionsequence_linear(femmodel);
40  }
41 }

◆ extrudefrombase_core()

void extrudefrombase_core ( FemModel *  femmodel)

Definition at line 12 of file extrudefrombase_core.cpp.

12  {
13 
14  /*Intermediaries*/
15  int elementtype,domaintype;
16 
17  if(VerboseSolution()) _printf0_(" extruding solution from base...\n");
18 
19  /*Get parameters*/
20  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
21  femmodel->parameters->FindParam(&elementtype,MeshElementtypeEnum);
22 
23  /*If this is a 2D horizontal domain, return (no need to extrude)*/
24  if(domaintype==Domain2DhorizontalEnum) return;
25 
26  /*Special method for Penta, otherwise call solution sequence*/
27  if(elementtype==PentaEnum){
28  int inputenum; femmodel->parameters->FindParam(&inputenum,InputToExtrudeEnum);
29  InputExtrudex(femmodel,inputenum,-1);
30  }
31  else{
32  /*Call on core computations: */
33  femmodel->SetCurrentConfiguration(ExtrudeFromBaseAnalysisEnum);
34  femmodel->UpdateConstraintsExtrudeFromBasex();
35  solutionsequence_linear(femmodel);
36  }
37 }

◆ extrudefromtop_core()

void extrudefromtop_core ( FemModel *  femmodel)

Definition at line 12 of file extrudefromtop_core.cpp.

12  {
13 
14  /*Intermediaries*/
15  int elementtype,domaintype;
16 
17  if(VerboseSolution()) _printf0_(" extruding solution from top...\n");
18 
19  /*Get parameters*/
20  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
21  femmodel->parameters->FindParam(&elementtype,MeshElementtypeEnum);
22 
23  /*If this is a 2D horizontal domain, return (no need to extrude)*/
24  if(domaintype==Domain2DhorizontalEnum) return;
25 
26  /*Special method for Penta, otherwise call solution sequence*/
27  if(elementtype==PentaEnum){
28  int inputenum; femmodel->parameters->FindParam(&inputenum,InputToExtrudeEnum);
29  InputExtrudex(femmodel,inputenum,+1);
30  }
31  else{
32  /*Call on core computations: */
33  femmodel->SetCurrentConfiguration(ExtrudeFromTopAnalysisEnum);
34  femmodel->UpdateConstraintsExtrudeFromTopx();
35  solutionsequence_linear(femmodel);
36  }
37 }

◆ balancethickness_core()

void balancethickness_core ( FemModel *  femmodel)

Definition at line 12 of file balancethickness_core.cpp.

12  {
13 
14  /*parameters: */
15  bool save_results;
16 
17  /*activate formulation: */
18  femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum);
19 
20  /*recover parameters: */
21  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
22 
23  if(VerboseSolution()) _printf0_(" call computational core:\n");
24  solutionsequence_linear(femmodel);
25 
26  if(save_results){
27  if(VerboseSolution()) _printf0_(" saving results\n");
28  int outputs = ThicknessEnum;
29  femmodel->RequestedOutputsx(&femmodel->results,&outputs,1);
30  }
31 
32 }

◆ balancethickness2_core()

void balancethickness2_core ( FemModel *  femmodel)

Definition at line 12 of file balancethickness2_core.cpp.

12  {
13 
14  /*parameters: */
15  bool save_results;
16  //IssmDouble l = 3.;
17 
18  /*recover parameters: */
19  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
20 
21  //if(VerboseSolution()) _printf0_("computing smooth surface slopes:\n");
22  //femmodel->parameters->SetParam(l,SmoothThicknessMultiplierEnum);
23  //femmodel->SetCurrentConfiguration(SmoothAnalysisEnum);
24  //femmodel->parameters->SetParam(SurfaceSlopeXEnum,InputToSmoothEnum);
25  //solutionsequence_linear(femmodel);
26  //femmodel->parameters->SetParam(SurfaceSlopeYEnum,InputToSmoothEnum);
27  //solutionsequence_linear(femmodel);
28  //surfaceslope_core(femmodel);
29 
30  femmodel->SetCurrentConfiguration(Balancethickness2AnalysisEnum);
31  solutionsequence_linear(femmodel);
32  //solutionsequence_nonlinear(femmodel,false);
33 
34  if(save_results){
35  if(VerboseSolution()) _printf0_(" saving results\n");
36  const int numoutputs = 1;
37  int outputs[numoutputs] = {ThicknessEnum};
38  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],numoutputs);
39  }
40 
41 }

◆ balancevelocity_core()

void balancevelocity_core ( FemModel *  femmodel)

Definition at line 12 of file balancevelocity_core.cpp.

12  {
13 
14  /*parameters: */
15  bool save_results;
16  IssmDouble l = 8.;
17 
18  /*recover parameters: */
19  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
20 
21  if(VerboseSolution()) _printf0_("computing smooth driving stress:\n");
22  femmodel->parameters->SetParam(l,SmoothThicknessMultiplierEnum);
23  femmodel->SetCurrentConfiguration(SmoothAnalysisEnum);
24  femmodel->parameters->SetParam(DrivingStressXEnum,InputToSmoothEnum);
25  solutionsequence_linear(femmodel);
26  femmodel->parameters->SetParam(DrivingStressYEnum,InputToSmoothEnum);
27  solutionsequence_linear(femmodel);
28 
29  if(VerboseSolution()) _printf0_(" call computational core:\n");
30  femmodel->SetCurrentConfiguration(BalancevelocityAnalysisEnum);
31  solutionsequence_linear(femmodel);
32 
33  if(save_results){
34  if(VerboseSolution()) _printf0_(" saving results\n");
35  int outputs[3] = {DrivingStressXEnum,DrivingStressYEnum,VelEnum};
36  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
37  }
38 
39 }

◆ slopecompute_core()

void slopecompute_core ( FemModel *  femmodel)

◆ steadystate_core()

void steadystate_core ( FemModel *  femmodel)

Definition at line 20 of file steadystate_core.cpp.

20  { //{{{
21 
22  /*intermediary: */
23  int step;
24  Vector<IssmDouble>* ug = NULL;
25  Vector<IssmDouble>* ug_old = NULL;
26  Vector<IssmDouble>* tg = NULL;
27  Vector<IssmDouble>* tg_old = NULL;
28 
29  /*parameters: */
30  bool save_results,isenthalpy;
31  int maxiter;
32  IssmDouble reltol;
33  int numoutputs = 0;
34  char** requested_outputs = NULL;
35 
36  /* recover parameters:*/
37  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
38  femmodel->parameters->FindParam(&maxiter,SteadystateMaxiterEnum);
39  femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
40  femmodel->parameters->FindParam(&reltol,SteadystateReltolEnum);
41  femmodel->parameters->SetParam(false,SaveResultsEnum);
42  femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum);
43  if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SteadystateRequestedOutputsEnum);
44 
45  /*intialize counters: */
46  step=1;
47 
48  for(;;){
49 
50  /* Compute first velocity, then temperature due to high sensitivity of temperature to velocity. */
51  if(VerboseSolution()) _printf0_("\n======================================================\n");
52  if(VerboseSolution()) _printf0_(" computing velocity and temperature for step: " << step << "\n");
53  if(VerboseSolution()) _printf0_("====================================================\n");
54 
55  if(VerboseSolution()) _printf0_("\n -- computing new velocity -- \n\n");
56  stressbalance_core(femmodel);
57  GetSolutionFromInputsx(&ug,femmodel);
58 
59  if(VerboseSolution()) _printf0_("\n -- computing new temperature --\n\n");
60  thermal_core(femmodel);
61  if(!isenthalpy)femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);/*Could be MeltingAnalysis...*/
62  GetSolutionFromInputsx(&tg,femmodel);
63 
64  if(step>1){
65  if(VerboseSolution()) _printf0_(" checking steadystate convergence\n");
66  if(steadystateconvergence(tg,tg_old,ug,ug_old,reltol)) break;
67  }
68  if(step>maxiter){
69  if(VerboseSolution()) _printf0_(" maximum number steadystate iterations " << maxiter << " reached\n");
70  break;
71  }
72 
73  /*update results and increase counter*/
74  delete tg_old;tg_old=tg;
75  delete ug_old;ug_old=ug;
76  step++;
77  }
78 
79  if(save_results){
80  if(VerboseSolution()) _printf0_(" saving results\n");
81  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
82  }
83 
84  /*Free ressources:*/
85  delete tg_old;
86  delete ug_old;
87  delete tg;
88  delete ug;
89  if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
90 }//}}}

◆ transient_core()

void transient_core ( FemModel *  femmodel)

Definition at line 19 of file transient_core.cpp.

19  {
20 
21  /*parameters: */
22  IssmDouble finaltime,dt,yts,starttime;
23  bool isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isesa,isslr,iscoupler,ismovingfront,isdamageevolution,ishydrology,isoceancoupling,iscontrol,isautodiff;
24  bool save_results,dakota_analysis;
25  int timestepping;
26  int output_frequency;
27  int sb_coupling_frequency;
28  int recording_frequency;
29  int domaintype,groundingline_migration,smb_model,amr_frequency,amr_restart;
30  int numoutputs;
31  Analysis *analysis = NULL;
32  char **requested_outputs = NULL;
33 
34  /*intermediary: */
35  int step;
36  IssmDouble time;
37 
38  /*first, figure out if there was a check point, if so, do a reset of the FemModel* femmodel structure. */
39  femmodel->parameters->FindParam(&recording_frequency,SettingsRecordingFrequencyEnum);
40  if(recording_frequency) femmodel->Restart();
41 
42  /*then recover parameters common to all solutions*/
43  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
44  femmodel->parameters->FindParam(&step,StepEnum);
45  femmodel->parameters->FindParam(&time,TimeEnum);
46  femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum);
47  femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum);
48  femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
49  femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
50  femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum);
51  femmodel->parameters->FindParam(&sb_coupling_frequency,SettingsSbCouplingFrequencyEnum);
52  femmodel->parameters->FindParam(&timestepping,TimesteppingTypeEnum);
53  femmodel->parameters->FindParam(&isstressbalance,TransientIsstressbalanceEnum);
54  femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum);
55  femmodel->parameters->FindParam(&issmb,TransientIssmbEnum);
56  femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
57  femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
58  femmodel->parameters->FindParam(&isesa,TransientIsesaEnum);
59  femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
60  femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
61  femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
62  femmodel->parameters->FindParam(&ismovingfront,TransientIsmovingfrontEnum);
63  femmodel->parameters->FindParam(&isoceancoupling,TransientIsoceancouplingEnum);
64  femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
65  femmodel->parameters->FindParam(&ishydrology,TransientIshydrologyEnum);
66  femmodel->parameters->FindParam(&amr_frequency,TransientAmrFrequencyEnum);
67  femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum);
68  femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum);
69  femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
70  if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
71  femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
72  if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum);
73 
74  #if defined(_HAVE_BAMG_) && !defined(_HAVE_AD_)
75  if(amr_frequency){
76  femmodel->parameters->FindParam(&amr_restart,AmrRestartEnum);
77  if(amr_restart) femmodel->ReMesh();
78  }
79  #endif
80 
81  #if defined(_HAVE_OCEAN_ )
82  if(isoceancoupling) OceanExchangeDatax(femmodel,true);
83  #endif
84 
85  IssmDouble output_value;
86  int num_dependents;
87  IssmPDouble *dependents;
88  DataSet* dependent_objects=NULL;
89  IssmDouble J=0.;
90 
91  if(iscontrol && isautodiff){
92 
93  femmodel->parameters->SetParam(starttime,TimeEnum);
94  femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
95  femmodel->parameters->FindParam(&dependent_objects,AutodiffDependentObjectsEnum);
96 
97  /*Go through our dependent variables, and compute the response:*/
98  dependents=xNew<IssmPDouble>(num_dependents);
99  }
100 
101  if(isslr) sealevelrise_core_geometry(femmodel);
102 
103  while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime.
104 
105  /*Time Increment*/
106  switch(timestepping){
107  case AdaptiveTimesteppingEnum:
108  femmodel->TimeAdaptx(&dt);
109  if(time+dt>finaltime) dt=finaltime-time;
110  femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
111  break;
112  case FixedTimesteppingEnum:
113  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
114  break;
115  default:
116  _error_("Time stepping \""<<EnumToStringx(timestepping)<<"\" not supported yet");
117  }
118  step+=1;
119  time+=dt;
120  femmodel->parameters->SetParam(time,TimeEnum);
121  femmodel->parameters->SetParam(step,StepEnum);
122 
123  if(VerboseSolution()) _printf0_("iteration " << step << "/" << ceil((finaltime-time)/dt)+step << " time [yr]: " << setprecision(4) << time/yts << " (time step: " << dt/yts << ")\n");
124  if(step%output_frequency==0 || (time >= finaltime - (yts*DBL_EPSILON)) || step==1)
125  save_results=true;
126  else
127  save_results=false;
128  femmodel->parameters->SetParam(save_results,SaveResultsEnum);
129 
130  #if defined(_HAVE_OCEAN_ )
131  if(isoceancoupling) OceanExchangeDatax(femmodel,false);
132  #endif
133 
134  if(isthermal && domaintype==Domain3DEnum){
135  if(issmb){
136  bool isenthalpy;
137  femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
138  femmodel->parameters->FindParam(&smb_model,SmbEnum);
139  if(isenthalpy){
140  if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
141  femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum);
142  ResetBoundaryConditions(femmodel,EnthalpyAnalysisEnum);
143  }
144  }
145  else{
146  if(smb_model==SMBpddEnum || smb_model==SMBd18opddEnum || smb_model==SMBpddSicopolisEnum){
147  femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);
148  ResetBoundaryConditions(femmodel,ThermalAnalysisEnum);
149  }
150  }
151  }
152  thermal_core(femmodel);
153  }
154  /* Using Hydrology dc coupled we need to compute smb in the hydrology inner time loop*/
155  if(issmb) {
156  if(VerboseSolution()) _printf0_(" computing smb\n");
157  smb_core(femmodel);
158  }
159 
160  if(ishydrology){
161  if(VerboseSolution()) _printf0_(" computing hydrology\n");
162  int hydrology_model;
163  hydrology_core(femmodel);
164  femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum);
165  if(hydrology_model!=HydrologydcEnum && issmb)smb_core(femmodel);
166  }
167 
168  if(isstressbalance && (step%sb_coupling_frequency==0 || step==1) ) {
169  if(VerboseSolution()) _printf0_(" computing stress balance\n");
170  stressbalance_core(femmodel);
171  }
172 
173  if(isdamageevolution) {
174  if(VerboseSolution()) _printf0_(" computing damage\n");
175  damage_core(femmodel);
176  }
177 
178  if(ismovingfront) {
179  if(VerboseSolution()) _printf0_(" computing moving front\n");
180  movingfront_core(femmodel);
181  }
182 
183  /* from here on, prepare geometry for next time step*/
184  //if(issmb) smb_core(femmodel);
185 
186  if(ismasstransport){
187  if(VerboseSolution()) _printf0_(" computing mass transport\n");
188  bmb_core(femmodel);
189  masstransport_core(femmodel);
190  femmodel->UpdateVertexPositionsx();
191  }
192 
193  if(isgroundingline){
194 
195  /*Start profiler*/
196  femmodel->profiler->Start(GROUNDINGLINECORE);
197 
198  if(VerboseSolution()) _printf0_(" computing new grounding line position\n");
199  GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);
200 
201  femmodel->parameters->SetParam(MaskOceanLevelsetEnum,InputToExtrudeEnum);
202  extrudefrombase_core(femmodel);
203  femmodel->parameters->SetParam(BaseEnum,InputToExtrudeEnum);
204  extrudefrombase_core(femmodel);
205  femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum);
206  extrudefrombase_core(femmodel);
207 
208  /*Stop profiler*/
209  femmodel->profiler->Stop(GROUNDINGLINECORE);
210 
211  if(save_results){
212  int outputs[3] = {SurfaceEnum,BaseEnum,MaskOceanLevelsetEnum};
213  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3);
214  }
215  }
216 
217  if(isgia){
218  if(VerboseSolution()) _printf0_(" computing glacial isostatic adjustment\n");
219  #ifdef _HAVE_GIA_
220  gia_core(femmodel);
221  #else
222  _error_("ISSM was not compiled with gia capabilities. Exiting");
223  #endif
224  }
225 
226  /*esa: */
227  if(isesa) esa_core(femmodel);
228 
229  /*Sea level rise: */
230  if(isslr) sealevelchange_core(femmodel);
231 
232  /*unload results*/
233  if(VerboseSolution()) _printf0_(" computing requested outputs\n");
234  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs,save_results);
235  if(isgroundingline && (groundingline_migration==AggressiveMigrationEnum || groundingline_migration==ContactEnum)){
236  int outputs[1] = {MaskOceanLevelsetEnum};
237  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1,save_results);
238  }
239 
240  if(save_results){
241  if(VerboseSolution()) _printf0_(" saving temporary results\n");
242  OutputResultsx(femmodel);
243  }
244 
245  if(recording_frequency && (step%recording_frequency==0)){
246  if(VerboseSolution()) _printf0_(" checkpointing model \n");
247  femmodel->CheckPoint();
248  }
249 
250  /*Adaptive mesh refinement*/
251  if(amr_frequency){
252 
253  #if !defined(_HAVE_AD_)
254  if(save_results) femmodel->WriteMeshInResults();
255  if(step%amr_frequency==0 && time<finaltime){
256  if(VerboseSolution()) _printf0_(" refining mesh\n");
257  femmodel->ReMesh();//Do not refine the last step
258  }
259 
260  #else
261  _error_("AMR not suppored with AD");
262  #endif
263  }
264 
265  if (iscontrol && isautodiff) {
266  /*Go through our dependent variables, and compute the response:*/
267  for(int i=0;i<dependent_objects->Size();i++){
268  DependentObject* dep=(DependentObject*)dependent_objects->GetObjectByOffset(i);
269  dep->Responsex(&output_value,femmodel);
270  dep->AddValue(output_value);
271  }
272  }
273  }
274 
275  if(!iscontrol || !isautodiff) femmodel->RequestedDependentsx();
276  if(iscontrol && isautodiff) femmodel->parameters->SetParam(dependent_objects,AutodiffDependentObjectsEnum);
277 
278  /*Free ressources:*/
279  if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
280 }

◆ dakota_core()

void dakota_core ( FemModel *  femmodel)

Definition at line 313 of file dakota_core.cpp.

313  {
314  _error_("dakota_core for versions of Dakota >=6 should not be used anymore! Use instead the issm_dakota executable!");
315 }

◆ ad_core()

void ad_core ( FemModel *  femmodel)

Definition at line 24 of file ad_core.cpp.

24  {
25 
26  /*diverse: */
27  int i;
28  int dummy;
29  int num_dependents=0;
30  int num_independents=0;
31  bool isautodiff,iscontrol;
32  char *driver = NULL;
33  size_t tape_stats[15];
34 
35  /*state variables: */
36  IssmDouble *axp = NULL;
37  double *xp = NULL;
38  int my_rank=IssmComm::GetRank();
39 
40  /*AD mode on?: */
41  femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
42  femmodel->parameters->FindParam(&iscontrol,InversionIscontrolEnum);
43 
44  if(isautodiff && !iscontrol){
45 
46  #if defined(_HAVE_ADOLC_)
47  if(VerboseAutodiff())_printf0_(" start ad core\n");
48 
49  /*First, stop tracing: */
50  trace_off();
51 
52  /*Print tape statistics so that user can kill this run if something is off already:*/
53  if(VerboseAutodiff()){ /*{{{*/
54  tapestats(my_rank,tape_stats); //reading of tape statistics
55  int commSize=IssmComm::GetSize();
56  int *sstats=new int[7];
57  sstats[0]=tape_stats[NUM_OPERATIONS];
58  sstats[1]=tape_stats[OP_FILE_ACCESS];
59  sstats[2]=tape_stats[NUM_LOCATIONS];
60  sstats[3]=tape_stats[LOC_FILE_ACCESS];
61  sstats[4]=tape_stats[NUM_VALUES];
62  sstats[5]=tape_stats[VAL_FILE_ACCESS];
63  sstats[6]=tape_stats[TAY_STACK_SIZE];
64  int *rstats=NULL;
65  if (my_rank==0) rstats=new int[commSize*7];
66  ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
67  if (my_rank==0) {
68  int offset=50;
69  int rOffset=(commSize/10)+1;
70  _printf_(" ADOLC statistics: \n");
71  _printf_(" "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
72  _printf_(" "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
73  _printf_(" "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
74  _printf_(" operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
75  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
76  for (int r=0;r<commSize;++r)
77  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+0] << (rstats[r*7+1]?" ->file":"") << "\n");
78  _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n");
79  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
80  for (int r=0;r<commSize;++r)
81  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+2] << (rstats[r*7+3]?" ->file":"") << "\n");
82  _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n");
83  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
84  for (int r=0;r<commSize;++r)
85  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+4] << (rstats[r*7+5]?" ->file":"") << "\n");
86  _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
87  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
88  for (int r=0;r<commSize;++r)
89  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+6] << (rstats[r*7+6]>tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n");
90  delete []rstats;
91  }
92  delete [] sstats;
93  } /*}}}*/
94 
95  /*retrieve parameters: */
96  femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
97  femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
98 
99  /*if no dependents, no point in running a driver: */
100  if(!(num_dependents*num_independents)) return;
101 
102  /*for adolc to run in parallel, we 0 out on rank~=0:*/
103  if (my_rank!=0){
104  num_dependents=0; num_independents=0;
105  }
106 
107  /*retrieve state variable: */
108  femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
109 
110  /* driver argument */
111  xp=xNew<double>(num_independents);
112  for(i=0;i<num_independents;i++){
113  xp[i]=reCast<double,IssmDouble>(axp[i]);
114  }
115 
116  /*get the EDF pointer:*/
117  ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
118 
119  /*Branch according to AD driver: */
120  femmodel->parameters->FindParam(&driver,AutodiffDriverEnum);
121 
122  if (strcmp(driver,"fos_forward")==0){ /*{{{*/
123 
124  int anIndepIndex;
125  double *tangentDir = NULL;
126  double *jacTimesTangentDir = NULL;
127  double *theOutput = NULL;
128 
129  /*retrieve direction index: */
130  femmodel->parameters->FindParam(&anIndepIndex,AutodiffFosForwardIndexEnum);
131 
132  if (anIndepIndex<0 || anIndepIndex>=num_independents) _error_("index value for AutodiffFosForwardIndexEnum should be in [0,num_independents-1]");
133 
134  tangentDir=xNewZeroInit<double>(num_independents);
135  tangentDir[anIndepIndex]=1.0;
136 
137  jacTimesTangentDir=xNew<double>(num_dependents);
138  theOutput=xNew<double>(num_dependents);
139 
140  /*set the forward method function pointer: */
141 #ifdef _HAVE_GSL_
142  anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx;
143 #endif
144 
145  /*call driver: */
146  fos_forward(my_rank,num_dependents,num_independents, 0, xp, tangentDir, theOutput, jacTimesTangentDir );
147 
148  /*add to results*/
149  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,jacTimesTangentDir,num_dependents,1,0,0.0));
150 
151  /*free resources :*/
152  xDelete(theOutput);
153  xDelete(jacTimesTangentDir);
154  xDelete(tangentDir);
155  } /*}}}*/
156  else if ((strcmp(driver,"fov_forward")==0) || (strcmp(driver,"fov_forward_all")==0)){ /*{{{*/
157 
158  int tangentDirNum;
159  int dummy;
160  int *indepIndices = NULL;
161  double **jacTimesSeed = NULL;
162  double **seed = NULL;
163  double *theOutput = NULL;
164  std::set<unsigned int> anIndexSet;
165 
166  /*retrieve directions:*/
167  if (strcmp(driver,"fov_forward_all")==0){
168  tangentDirNum=num_independents;
169  indepIndices=xNewZeroInit<int>(tangentDirNum);
170  for(i=0;i<num_independents;i++)indepIndices[i]=1;
171  }
172  else{
173  femmodel->parameters->FindParam(&indepIndices,&tangentDirNum,&dummy,AutodiffFovForwardIndicesEnum);
174  }
175 
176  /*Some checks: */
177  if (tangentDirNum<1 || tangentDirNum>num_independents) _error_("tangentDirNum should be in [1,num_independents]");
178 
179  /* full Jacobian or Jacobian projection:*/
180  jacTimesSeed=xNew<double>(num_dependents,tangentDirNum);
181 
182  /*set the forward method function pointers: */
183 #ifdef _HAVE_GSL_
184  anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx;
185 #endif
186  // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
187 
188  /*seed matrix: */
189  seed=xNewZeroInit<double>(num_independents,tangentDirNum);
190 
191  /*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/
192  for (int i=0; i<tangentDirNum; ++i) {
193  /* make sure the index is in range*/
194  if (indepIndices[i]>num_independents) {
195  _error_("indepIndices values must be in [0,num_independents-1]");
196  }
197  if (anIndexSet.find(indepIndices[i])!=anIndexSet.end()) {
198  _error_("duplicate indepIndices values are not allowed until we implement Jacobian decompression");
199  }
200  anIndexSet.insert(indepIndices[i]);
201  /* now populate the seed matrix from the set of independent indices;
202  * simple setup with a single 1.0 per column and at most a single 1.0 per row*/
203  seed[indepIndices[i]][i]=1.0;
204  }
205 
206  /*allocate output: */
207  theOutput=xNew<double>(num_dependents);
208 
209  /*call driver: */
210  fov_forward(my_rank,num_dependents,num_independents, tangentDirNum, xp, seed, theOutput, jacTimesSeed );
211  /*Free resources: */
212  xDelete(theOutput);
213  xDelete(indepIndices);
214  xDelete(seed);
215 
216  /*add to results: */
217  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,*jacTimesSeed,num_dependents*tangentDirNum,1,0,0.0));
218 
219  /*Free resources: */
220  xDelete(jacTimesSeed);
221  xDelete(indepIndices);
222  } /*}}}*/
223  else if (strcmp(driver,"fos_reverse")==0) { /*{{{*/
224 
225  int aDepIndex=0;
226  double *aWeightVector=NULL;
227  double *weightVectorTimesJac=NULL;
228 
229  /*retrieve direction index: */
230  femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum);
231  aWeightVector=xNewZeroInit<double>(num_dependents);
232  if (my_rank==0) {
233  if (aDepIndex<0 || aDepIndex>=num_dependents) _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]");
234  aWeightVector[aDepIndex]=1.0;
235  }
236  weightVectorTimesJac=xNew<double>(num_independents);
237 
238  /*set the forward method function pointer: */
239 #ifdef _HAVE_GSL_
240  anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
241 #endif
242 #ifdef _HAVE_MUMPS_
243  anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
244 #endif
245 
246  /*call driver: */
247  fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
248 
249  /*add to results*/
250  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,weightVectorTimesJac,num_independents,1,0,0.0));
251 
252  /*free resources :*/
253  xDelete(weightVectorTimesJac);
254  xDelete(aWeightVector);
255  } /*}}}*/
256  else if ((strcmp(driver,"fov_reverse")==0) || (strcmp(driver,"fov_reverse_all")==0)){ /*{{{*/
257 
258  int* depIndices=NULL;
259  int weightNum;
260  int dummy;
261  double **weightsTimesJac=NULL;
262  double **weights=NULL;
263  std::set<unsigned int> anIndexSet;
264 
265  /*retrieve directions:*/
266  if (strcmp(driver,"fov_reverse_all")==0){
267  weightNum=num_dependents;
268  depIndices=xNewZeroInit<int>(weightNum);
269  for(i=0;i<num_dependents;i++)depIndices[i]=1;
270  }
271  else{
272  femmodel->parameters->FindParam(&depIndices,&weightNum,&dummy,AutodiffFovForwardIndicesEnum);
273  }
274 
275  /*Some checks: */
276  if (weightNum<1 || weightNum>num_dependents) _error_("tangentDirNum should be in [1,num_dependents]");
277 
278  /* full Jacobian or Jacobian projection:*/
279  weightsTimesJac=xNew<double>(weightNum,num_independents);
280 
281  /*set the forward method function pointers: */
282  #ifdef _HAVE_GSL_
283  anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx;
284  #endif
285 
286  /*seed matrix: */
287  weights=xNewZeroInit<double>(weightNum,num_dependents);
288 
289  /*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/
290  for (int i=0; i<weightNum; ++i) {
291  /* make sure the index is in range*/
292  if (depIndices[i]>num_dependents) {
293  _error_("depIndices values must be in [0,num_dependents-1]");
294  }
295  if (anIndexSet.find(depIndices[i])!=anIndexSet.end()) {
296  _error_("duplicate depIndices values are not allowed until we implement Jacobian decompression");
297  }
298  anIndexSet.insert(depIndices[i]);
299  /* now populate the seed matrix from the set of independent indices;
300  * simple setup with a single 1.0 per column and at most a single 1.0 per row*/
301  weights[depIndices[i]][i]=1.0;
302  }
303 
304  /*call driver: */
305  fov_reverse(my_rank,num_dependents,num_independents, weightNum, weights, weightsTimesJac );
306 
307  /*add to results: */
308  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,*weightsTimesJac,weightNum*num_independents,1,0,0.0));
309 
310  /*Free resources: */
311  xDelete(weights);
312  xDelete(weightsTimesJac);
313  xDelete(depIndices);
314  } /*}}}*/
315  else _error_("driver: " << driver << " not yet supported!");
316 
317  if(VerboseAutodiff())_printf0_(" end AD core\n");
318 
319  /*Free resources: */
320  xDelete(xp);
321  xDelete(axp);
322  xDelete(driver);
323 
324  #elif defined(_HAVE_CODIPACK_)
325  if(VerboseAutodiff())_printf0_(" start CoDiPack ad core\n");
326 
327  /*First, stop tracing: */
328  auto& tape_codi = IssmDouble::getGlobalTape();
329  tape_codi.setPassive();
330 
331  if(VerboseAutodiff()){ /*{{{*/
332  if(my_rank == 0) {
333  // FIXME codi "just because" for now
334  tape_codi.printStatistics(std::cout);
335  codi_global.print(std::cout);
336  }
337  }
338 
339  /*retrieve parameters: */
340  femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
341  femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
342 
343  /*if no dependents, no point in running a driver: */
344  if(!(num_dependents*num_independents)) return;
345 
346  /*Branch according to AD driver: */
347  femmodel->parameters->FindParam(&driver,AutodiffDriverEnum);
348  if(VerboseAutodiff())_printf0_(" driver: " << driver << "\n");
349 
350  if (strcmp(driver,"fos_reverse")==0) { /*{{{*/
351  if(VerboseAutodiff())_printf0_(" CoDiPack fos_reverse\n");
352  int aDepIndex=0;
353  double *weightVectorTimesJac=NULL;
354 
355  /*retrieve direction index: */
356  femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum);
357  if (my_rank==0) {
358  if (aDepIndex<0 || aDepIndex>=num_dependents
359  || codi_global.output_indices.size() <= aDepIndex){
360  _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]");
361  }
362  tape_codi.setGradient(codi_global.output_indices[aDepIndex], 1.0);
363  }
364 
365  tape_codi.evaluate();
366 
367  weightVectorTimesJac=xNew<double>(num_independents);
368  /*call driver: */
369  auto in_size = codi_global.input_indices.size();
370  for(size_t i = 0; i < in_size; ++i) {
371  weightVectorTimesJac[i] = tape_codi.getGradient(codi_global.input_indices[i]);
372  }
373 
374  /*add to results*/
375  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,weightVectorTimesJac,num_independents,1,0,0.0));
376 
377  /*free resources :*/
378  xDelete(weightVectorTimesJac);
379  } /*}}}*/
380  else _error_("driver: " << driver << " not yet supported!");
381 
382  if(VerboseAutodiff())_printf0_(" end CoDiPack ad core\n");
383  xDelete(driver);
384  #else
385  _error_("Should not be requesting AD drivers when an AD library is not available!");
386  #endif
387  }
388 }

◆ adgradient_core()

void adgradient_core ( FemModel *  femmodel)

Definition at line 22 of file adgradient_core.cpp.

22  {
23 
24  /*diverse: */
25  int i;
26  int dummy;
27  int num_dependents=0;
28  int num_dependents_old=0;
29  int num_independents=0;
30  bool isautodiff = false;
31  int aDepIndex=0;
32  int my_rank=IssmComm::GetRank();
33 
34  /*state variables: */
35  IssmDouble *axp = NULL;
36  IssmPDouble *xp = NULL;
37 
38  /*intermediary: */
39  IssmPDouble *aWeightVector=NULL;
40  IssmPDouble *weightVectorTimesJac=NULL;
41 
42  /*output: */
43  IssmPDouble *totalgradient=NULL;
44 
45  /*AD mode on?: */
46  femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum);
47 
48  if(isautodiff){
49 
50  #if defined(_HAVE_ADOLC_)
51  if(VerboseAutodiff())_printf0_(" start ad core\n");
52 
53  /*First, stop tracing: */
54  trace_off();
55 
56  if(VerboseAutodiff()){ /*{{{*/
57  size_t tape_stats[15];
58  tapestats(my_rank,tape_stats); //reading of tape statistics
59  int commSize=IssmComm::GetSize();
60  int *sstats=new int[7];
61  sstats[0]=tape_stats[NUM_OPERATIONS];
62  sstats[1]=tape_stats[OP_FILE_ACCESS];
63  sstats[2]=tape_stats[NUM_LOCATIONS];
64  sstats[3]=tape_stats[LOC_FILE_ACCESS];
65  sstats[4]=tape_stats[NUM_VALUES];
66  sstats[5]=tape_stats[VAL_FILE_ACCESS];
67  sstats[6]=tape_stats[TAY_STACK_SIZE];
68  int *rstats=NULL;
69  if (my_rank==0) rstats=new int[commSize*7];
70  ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm());
71  if (my_rank==0) {
72  int offset=50;
73  int rOffset=(commSize/10)+1;
74  _printf_(" ADOLC statistics: \n");
75  _printf_(" "<<setw(offset)<<left<<"#independents: " <<setw(12)<<right<<tape_stats[NUM_INDEPENDENTS] << "\n");
76  _printf_(" "<<setw(offset)<<left<<"#dependents: " <<setw(12)<<right<<tape_stats[NUM_DEPENDENTS] << "\n");
77  _printf_(" "<<setw(offset)<<left<<"max #live active variables: " <<setw(12)<<right<<tape_stats[NUM_MAX_LIVES] << "\n");
78  _printf_(" operations: entry size "<< sizeof(unsigned char) << " Bytes \n");
79  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffObufsizeEnum) " <<setw(12)<<right<<tape_stats[OP_BUFFER_SIZE] << "\n");
80  for (int r=0;r<commSize;++r)
81  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+0] << (rstats[r*7+1]?" ->file":"") << "\n");
82  _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n");
83  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffLbufsizeEnum) " <<setw(12)<<right<<tape_stats[LOC_BUFFER_SIZE] << "\n");
84  for (int r=0;r<commSize;++r)
85  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+2] << (rstats[r*7+3]?" ->file":"") << "\n");
86  _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n");
87  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffCbufsizeEnum) " <<setw(12)<<right<<tape_stats[VAL_BUFFER_SIZE] << "\n");
88  for (int r=0;r<commSize;++r)
89  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+4] << (rstats[r*7+5]?" ->file":"") << "\n");
90  _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n");
91  _printf_(" "<<setw(offset)<<left<<" #entries in buffer (AutodiffTbufsizeEnum) " <<setw(12)<<right<<tape_stats[TAY_BUFFER_SIZE] << "\n");
92  for (int r=0;r<commSize;++r)
93  _printf_(" ["<<setw(rOffset)<<right<<r<<"]"<<setw(offset-rOffset-4)<<left<<" #entries total" <<setw(12)<<right<<rstats[r*7+6] << (rstats[r*7+6]>tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n");
94  delete []rstats;
95  }
96  delete [] sstats;
97  } /*}}}*/
98 
99  /*retrieve parameters: */
100  femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum);
101  femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum);
102 
103  /*if no dependents, no point in running a driver: */
104  if(!(num_dependents*num_independents)) return;
105 
106  /*for adolc to run in parallel, we 0 out on rank~=0. But we still keep track of num_dependents:*/
107  num_dependents_old=num_dependents;
108  if (my_rank!=0){
109  num_dependents=0; num_independents=0;
110  }
111 
112  /*retrieve state variable: */
113  femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum);
114 
115  /* driver argument */
116  xp=xNew<double>(num_independents);
117  for(i=0;i<num_independents;i++){
118  xp[i]=reCast<double,IssmDouble>(axp[i]);
119  }
120 
121  /*get the EDF pointer:*/
122  ext_diff_fct *anEDF_for_solverx_p=xDynamicCast<GenericParam<Adolc_edf> * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p;
123 
124  /* these are always needed regardless of the interpreter */
125  anEDF_for_solverx_p->dp_x=xNew<double>(anEDF_for_solverx_p->max_n);
126  anEDF_for_solverx_p->dp_y=xNew<double>(anEDF_for_solverx_p->max_m);
127 
128  /* Ok, now we are going to call the fos_reverse in a loop on the index, from 0 to num_dependents, so
129  * as to generate num_dependents gradients: */
130 
131  /*Initialize outputs: */
132  totalgradient=xNewZeroInit<IssmPDouble>(num_independents);
133 
134  for(aDepIndex=0;aDepIndex<num_dependents_old;aDepIndex++){
135 
136  /*initialize direction index in the weights vector: */
137  aWeightVector=xNewZeroInit<IssmPDouble>(num_dependents);
138  if (my_rank==0) aWeightVector[aDepIndex]=1.0;
139 
140  /*initialize output gradient: */
141  weightVectorTimesJac=xNew<IssmPDouble>(num_independents);
142 
143  /*set the forward method function pointer: */
144  #ifdef _HAVE_GSL_
145  anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx;
146  #endif
147  #ifdef _HAVE_MUMPS_
148  anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF;
149  #endif
150 
151  anEDF_for_solverx_p->dp_U=xNew<IssmPDouble>(anEDF_for_solverx_p->max_m);
152  anEDF_for_solverx_p->dp_Z=xNew<IssmPDouble>(anEDF_for_solverx_p->max_n);
153 
154  /*call driver: */
155  fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac );
156 
157  /*Add to totalgradient: */
158  if(my_rank==0)for(i=0;i<num_independents;i++)totalgradient[i]+=weightVectorTimesJac[i];
159 
160  /*free resources :*/
161  xDelete(weightVectorTimesJac);
162  xDelete(aWeightVector);
163  }
164 
165  /*add totalgradient to results*/
166  femmodel->results->AddObject(new GenericExternalResult<IssmPDouble*>(femmodel->results->Size()+1,AutodiffJacobianEnum,totalgradient,num_independents,1,1,0.0));
167 
168  if(VerboseAutodiff())_printf0_(" end ad core\n");
169 
170  /* delete the allocated space for the parameters and free ressources:{{{*/
171  xDelete(anEDF_for_solverx_p->dp_x);
172  xDelete(anEDF_for_solverx_p->dp_X);
173  xDelete(anEDF_for_solverx_p->dpp_X);
174  xDelete(anEDF_for_solverx_p->dp_y);
175  xDelete(anEDF_for_solverx_p->dp_Y);
176  xDelete(anEDF_for_solverx_p->dpp_Y);
177  xDelete(anEDF_for_solverx_p->dp_U);
178  xDelete(anEDF_for_solverx_p->dpp_U);
179  xDelete(anEDF_for_solverx_p->dp_Z);
180  xDelete(anEDF_for_solverx_p->dpp_Z);
181  xDelete(xp);
182  xDelete(totalgradient);
183  xDelete(axp); /*}}}*/
184 
185  #elif defined(_HAVE_CODIPACK_)
186  fprintf(stderr, "*** Codipack adgradient_core()\n");
187  #else
188  _error_("Should not be requesting AD drivers when an AD library is not available!");
189  #endif
190  }
191 }

◆ dummy_core()

void dummy_core ( FemModel *  femmodel)

Definition at line 7 of file dummy_core.cpp.

7  {
8 
9  /*We do not do anything*/
10 
11 }

◆ gia_core()

void gia_core ( FemModel *  femmodel)

Definition at line 13 of file gia_core.cpp.

13  {
14 
15  Vector<IssmDouble> *wg = NULL;
16  Vector<IssmDouble> *dwdtg = NULL;
17  IssmDouble *x = NULL;
18  IssmDouble *y = NULL;
19  int nummodels,giamodel;
20  IssmDouble modelid;
21 
22  /*parameters: */
23  bool save_results;
24  int gsize;
25 
26  /*Start profiler*/
27  femmodel->profiler->Start(GIACORE);
28 
29  /*Recover some parameters: */
30  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
31 
32  if(VerboseSolution()) _printf0_(" computing GIA\n");
33 
34  /*Call on core computations: */
35  femmodel->SetCurrentConfiguration(GiaAnalysisEnum);
36 
37  /*Figure out type of gia model:*/
38  femmodel->parameters->FindParam(&giamodel,GiaModelEnum);
39 
40  if(giamodel==1){ //GiaIvins
41  /*Figure out size of g-set deflection vector and allocate solution vector: */
42  gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
43  wg = new Vector<IssmDouble>(gsize);
44  dwdtg = new Vector<IssmDouble>(gsize);
45 
46  /*first, recover x and y vectors from vertices: */
47  VertexCoordinatesx(&x,&y,NULL,femmodel->vertices); //no need for z coordinate
48 
49  /*call the main module: */
50  femmodel->Deflection(wg,dwdtg,x,y);
51 
52  /*assemble vector: */
53  wg->Assemble();
54  dwdtg->Assemble();
55 
56  InputUpdateFromVectorx(femmodel,wg,UGiaEnum,VertexSIdEnum);
57  InputUpdateFromVectorx(femmodel,dwdtg,UGiaRateEnum,VertexSIdEnum);
58 
59  if(save_results){
60  if(VerboseSolution()) _printf0_(" saving results\n");
61  int outputs[2] = {UGiaEnum,UGiaRateEnum};
62  femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2);
63  }
64 
65  xDelete<IssmDouble>(x);
66  xDelete<IssmDouble>(y);
67  }
68  else if(giamodel==2){ //GiaCaron
69  /*not implemneted yet*/
70  }
71  else if(giamodel==3){ //GiaMme
72 
73  /*Go grab the values for Ngia and Ugia in the offline GiaMme class:*/
74  femmodel->parameters->FindParam(&modelid,GiaModelidEnum);
75  modelid--; //from matlab
76 
77  /*find the Ngia and Ugia dataset:*/
78  DatasetInput2* dataset_input_ngia=femmodel->inputs2->GetDatasetInput2(GiaMmeNgiaEnum);
79  DatasetInput2* dataset_input_ugia=femmodel->inputs2->GetDatasetInput2(GiaMmeUgiaEnum);
80  /*Go find the modelid'th input:*/
81  TriaInput2* tria_input_ngia=dataset_input_ngia->GetTriaInputByOffset(modelid);
82  TriaInput2* tria_input_ugia=dataset_input_ugia->GetTriaInputByOffset(modelid);
83 
84  /*Plug into SealevelUGiaRate and SealevelNGiaRate inputs:*/
85  Input2* tria_input_copy_ngia=tria_input_ngia->copy();
86  Input2* tria_input_copy_ugia=tria_input_ugia->copy();
87  tria_input_copy_ngia->ChangeEnum(NGiaRateEnum);
88  tria_input_copy_ugia->ChangeEnum(UGiaRateEnum);
89  femmodel->inputs2->AddInput(tria_input_copy_ngia);
90  femmodel->inputs2->AddInput(tria_input_copy_ugia);
91  }
92 
93  /*End profiler*/
94  femmodel->profiler->Stop(GIACORE);
95 }

◆ love_core()

void love_core ( FemModel *  femmodel)

Definition at line 11 of file love_core.cpp.

11  {
12 
13  Vector<IssmDouble> *wg = NULL;
14  Vector<IssmDouble> *dwdtg = NULL;
15  IssmDouble *x = NULL;
16  IssmDouble *y = NULL;
17 
18  /*love parameters: */
19  IssmDouble *frequencies = NULL;
20  int nfreq,dummy;
21  int sh_nmin,sh_nmax;
22  IssmDouble g0,r0,mu0;
23  bool allow_layer_deletion;
24  bool love_kernels;
25  int forcing_type;
26  bool verbosemod = (int)VerboseModule();
27 
28  /*parameters: */
29  bool save_results;
30 
31  if(VerboseSolution()) _printf0_(" computing LOVE numbers\n");
32 
33  /*Recover some parameters: */
34  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
35 
36  /*recover love number parameters: */
37  femmodel->parameters->FindParam(&nfreq,LoveNfreqEnum);
38  femmodel->parameters->FindParam(&frequencies,&dummy,LoveFrequenciesEnum); _assert_(nfreq==dummy);
39  femmodel->parameters->FindParam(&sh_nmax,LoveShNmaxEnum);
40  femmodel->parameters->FindParam(&sh_nmin,LoveShNminEnum);
41  femmodel->parameters->FindParam(&g0,LoveG0Enum);
42  femmodel->parameters->FindParam(&r0,LoveR0Enum);
43  femmodel->parameters->FindParam(&mu0,LoveMu0Enum);
44  femmodel->parameters->FindParam(&allow_layer_deletion,LoveAllowLayerDeletionEnum);
45  femmodel->parameters->FindParam(&love_kernels,LoveKernelsEnum);
46  femmodel->parameters->FindParam(&forcing_type,LoveForcingTypeEnum);
47 
48  /*recover materials parameters: there is only one Matlitho, chase it down the hard way:*/
49  Matlitho* matlitho=NULL;
50  for (int i=femmodel->materials->Size()-1;i>=0;i--){
51  Material* material=xDynamicCast<Material*>(femmodel->materials->GetObjectByOffset(i));
52  if (material->ObjectEnum()==MatlithoEnum)matlitho=xDynamicCast<Matlitho*>(material);
53  }
54  _assert_(matlitho);
55 
56  /*Initialize three love matrices: geoid, vertical and horizontal displacement (real and imaginary parts) */
57  IssmDouble* LoveKr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
58  IssmDouble* LoveHr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
59  IssmDouble* LoveLr = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
60  IssmDouble* LoveKi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
61  IssmDouble* LoveHi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
62  IssmDouble* LoveLi = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1));
63 
64  /*Initialize love kernels (real and imaginary parts): */
65  IssmDouble* LoveKernelsReal = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
66  IssmDouble* LoveKernelsImag = xNewZeroInit<IssmDouble>(nfreq*(sh_nmax+1)*(matlitho->numlayers+1)*6);
67 
68  /*call the main module: */
69  FourierLoveCorex(LoveKr,LoveKi,LoveHr,LoveHi,LoveLr,LoveLi,LoveKernelsReal,LoveKernelsImag, //output
70  nfreq,frequencies,sh_nmax,sh_nmin,g0,r0,mu0,allow_layer_deletion,forcing_type,verbosemod, //parameter inputs
71  matlitho->numlayers, matlitho->radius, matlitho->viscosity, matlitho->lame_lambda, matlitho->lame_mu,
72  matlitho->burgers_viscosity, matlitho->burgers_mu, matlitho->density, matlitho->isburgers, matlitho->issolid //matlitho inputs
73  );
74 
75  /*Add love matrices to results:*/
76  femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKrEnum,LoveKr,sh_nmax+1,nfreq,0,0));
77  femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHrEnum,LoveHr,sh_nmax+1,nfreq,0,0));
78  femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLrEnum,LoveLr,sh_nmax+1,nfreq,0,0));
79  femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKiEnum,LoveKi,sh_nmax+1,nfreq,0,0));
80  femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveHiEnum,LoveHi,sh_nmax+1,nfreq,0,0));
81  femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveLiEnum,LoveLi,sh_nmax+1,nfreq,0,0));
82  /*Only when love_kernels is set unity*/
83  if (love_kernels==1) {
84  femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsRealEnum,LoveKernelsReal,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
85  femmodel->results->AddObject(new GenericExternalResult<IssmDouble*>(femmodel->results->Size()+1,LoveKernelsImagEnum,LoveKernelsImag,(sh_nmax+1)*(matlitho->numlayers+1)*6,nfreq,0,0));
86  }
87 
88  /*Free ressources:*/
89  xDelete<IssmDouble>(frequencies);
90  xDelete<IssmDouble>(LoveKr);
91  xDelete<IssmDouble>(LoveHr);
92  xDelete<IssmDouble>(LoveLr);
93  xDelete<IssmDouble>(LoveKi);
94  xDelete<IssmDouble>(LoveHi);
95  xDelete<IssmDouble>(LoveLi);
96  xDelete<IssmDouble>(LoveKernelsReal);
97  xDelete<IssmDouble>(LoveKernelsImag);
98 }

◆ esa_core()

void esa_core ( FemModel *  femmodel)

Definition at line 12 of file esa_core.cpp.

12  { /*{{{*/
13 
14  /*Start profiler*/
15  femmodel->profiler->Start(ESACORE);
16 
17  Vector<IssmDouble> *U_radial = NULL;
18  Vector<IssmDouble> *U_north = NULL;
19  Vector<IssmDouble> *U_east = NULL;
20  Vector<IssmDouble> *U_x = NULL;
21  Vector<IssmDouble> *U_y = NULL;
22  bool save_results,isesa,iscoupler;
23  int domaintype;
24  int solution_type;
25  int numoutputs = 0;
26  char **requested_outputs = NULL;
27 
28  /*additional parameters: */
29  int gsize;
30  bool spherical=true;
31  IssmDouble *latitude = NULL;
32  IssmDouble *longitude = NULL;
33  IssmDouble *radius = NULL;
34  IssmDouble *xx = NULL;
35  IssmDouble *yy = NULL;
36  IssmDouble *zz = NULL;
37 
38  /*Recover some parameters: */
39  femmodel->parameters->FindParam(&domaintype,DomainTypeEnum);
40  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
41  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
42  femmodel->parameters->FindParam(&isesa,TransientIsesaEnum);
43  femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
44 
45  /* recover coordinates of vertices: */
46  VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
47  VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
48 
49  /*Figure out size of g-set deflection vector and allocate solution vector: */
50  gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
51 
52  /*several cases here, depending on value of iscoupler and isesa:
53  solution_type == EsaSolutionEnum) we are running elastic adjustment core (no coupler)
54  ( !iscoupler & !isesa) we are not interested in being here :)
55  ( !iscoupler & isesa) we are running in uncoupled mode
56  ( iscoupler & isesa) we are running in coupled mode, and better be earth
57  ( iscoupler & !isesa) we are running in coupled mode, and better be an ice cap
58  */
59 
60  if(solution_type==EsaSolutionEnum){
61  isesa=1;
62  iscoupler=0;
63  }
64 
65  /*early return: */
66  if( !iscoupler & !isesa) return; //we are not interested in being here :)
67 
68  /*In what follows we assume we are all running esa, either in coupled, or uncoupled mode:*/
69  if(VerboseSolution()) _printf0_(" computing elastic adjustment\n");
70 
71  /*set configuration: */
72  if(isesa)femmodel->SetCurrentConfiguration(EsaAnalysisEnum);
73 
74  if(VerboseSolution()) _printf0_(" computing elastic geodetic core\n");
75  if(isesa){
76 
77  /*compute components of 3-D crustal motion: */
78  /*Initialize:*/
79  U_radial = new Vector<IssmDouble>(gsize);
80  U_north = new Vector<IssmDouble>(gsize);
81  U_east = new Vector<IssmDouble>(gsize);
82  U_x = new Vector<IssmDouble>(gsize);
83  U_y = new Vector<IssmDouble>(gsize);
84 
85  /*call the geodetic main modlule:*/
86  if(domaintype==Domain3DsurfaceEnum){
87  femmodel->EsaGeodetic3D(U_radial,U_north,U_east,latitude,longitude,radius,xx,yy,zz);
88  }
89  if(domaintype==Domain2DhorizontalEnum){
90  femmodel->EsaGeodetic2D(U_radial,U_north,U_east,U_x,U_y,xx,yy);
91  InputUpdateFromVectorx(femmodel,U_x,EsaXmotionEnum,VertexSIdEnum);
92  InputUpdateFromVectorx(femmodel,U_y,EsaYmotionEnum,VertexSIdEnum);
93  }
94 
95  /*get results into elements:*/
96  InputUpdateFromVectorx(femmodel,U_radial,EsaUmotionEnum,VertexSIdEnum); // radial displacement
97  InputUpdateFromVectorx(femmodel,U_north,EsaNmotionEnum,VertexSIdEnum); // north motion
98  InputUpdateFromVectorx(femmodel,U_east,EsaEmotionEnum,VertexSIdEnum); // east motion
99 
100  if(save_results){
101  if(VerboseSolution()) _printf0_(" saving results\n");
102  femmodel->parameters->FindParam(&requested_outputs,&numoutputs,EsaRequestedOutputsEnum);
103  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
104  }
105 
106  if(solution_type==EsaSolutionEnum)femmodel->RequestedDependentsx();
107 
108  /*Free ressources:*/
109  delete U_radial;
110  delete U_north;
111  delete U_east;
112  delete U_x;
113  delete U_y;
114  if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
115  }
116 
117  /*End profiler*/
118  femmodel->profiler->Stop(ESACORE);
119 
120 }

◆ smb_core()

void smb_core ( FemModel *  femmodel)

Definition at line 12 of file smb_core.cpp.

12  {
13 
14  /*Start profiler*/
15  femmodel->profiler->Start(SMBCORE);
16 
17  /*parameters: */
18  Analysis* analysis=NULL;
19  int smb_model;
20  int numoutputs;
21  bool save_results;
22  int solution_type;
23  char** requested_outputs = NULL;
24 
25  /*activate configuration*/
26  femmodel->SetCurrentConfiguration(SmbAnalysisEnum);
27 
28  /*recover parameters: */
29  femmodel->parameters->FindParam(&smb_model,SmbEnum);
30  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
31  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
32  femmodel->parameters->FindParam(&numoutputs,SmbNumRequestedOutputsEnum);
33  if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SmbRequestedOutputsEnum);
34 
35  /*sub steping specifics*/
36  int dtslices;
37  int numaveragedinput;
38  femmodel->parameters->FindParam(&dtslices,SmbStepsPerStepEnum);
39  /*intermediaries to deal with averaging*/
40  static const int substeplist[2] = {SmbMassBalanceSubstepEnum,SmbRunoffSubstepEnum};
41  static const int transientlist[2] = {SmbMassBalanceTransientEnum,SmbRunoffTransientEnum};
42  static const int averagelist[2] = {SmbMassBalanceEnum,SmbRunoffEnum};
43  std::vector<int> substepinput;
44  std::vector<int> transientinput;
45  std::vector<int> averagedinput;
46 
47  /*define which variable needs to be averaged on the sub-timestep and initialize as needed*/
48  if(smb_model==SMBgradientscomponentsEnum){
49  numaveragedinput = 2;
50  substepinput.assign(substeplist,substeplist+2);
51  transientinput.assign(transientlist,transientlist+2);
52  averagedinput.assign(averagelist,averagelist+2);
53  }
54 
55  /*if yes compute necessary intermediaries and start looping*/
56  if (dtslices>1){
57  int substep,smb_averaging;
58  IssmDouble global_time,subtime,yts;
59  IssmDouble dt,subdt;
60 
61  femmodel->parameters->FindParam(&global_time,TimeEnum);
62  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
63  femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
64  femmodel->parameters->FindParam(&smb_averaging,SmbAveragingEnum);
65 
66  subtime=global_time-dt; //getting the time back to the start of the timestep
67  subdt=dt/dtslices; //computing substep from dt and a divider
68  substep=0;
69  femmodel->parameters->SetParam(subdt,TimesteppingTimeStepEnum);
70 
71  femmodel->InitTransientInputx(&transientinput[0],numaveragedinput);
72  analysis = new SmbAnalysis();
73  while(substep<dtslices){ //loop on sub dts
74  substep+=1;
75  subtime+=subdt;
76  femmodel->parameters->SetParam(subtime,TimeEnum);
77  if(VerboseSolution()) _printf0_("sub iteration " << substep << "/" << dtslices << " time [yr]: " << setprecision(4) << subtime/yts << " (time step: " << subdt/yts << ")\n");
78  if(VerboseSolution()) _printf0_(" computing smb\n");
79  if(VerboseSolution()) _printf0_(" Calling core\n");
80  analysis->Core(femmodel);
81  /*If we have a sub-timestep we store the substep inputs in a transient input here*/
82  femmodel->StackTransientInputx(&substepinput[0],&transientinput[0],subtime,numaveragedinput);
83  }
84  delete analysis;
85  /*averaging the transient input*/
86  femmodel->AverageTransientInputx(&transientinput[0],&averagedinput[0],global_time-dt,subtime,numaveragedinput,smb_averaging);
87  /*and reset timesteping variables to original*/
88  femmodel->parameters->SetParam(global_time,TimeEnum);
89  femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum);
90  }
91  else{
92  if(VerboseSolution()) _printf0_(" computing smb \n");
93  analysis = new SmbAnalysis();
94  analysis->Core(femmodel);
95  /*If no substeps are present we want to duplicate the computed substep enum for coupling purposes*/
96  if(smb_model==SMBgradientscomponentsEnum){
97  for(int i=0;i<numaveragedinput;i++){
98  InputDuplicatex(femmodel,substepinput[i],averagedinput[i]);
99  }
100  }
101  delete analysis;
102  }
103 
104  if(save_results){
105  if(VerboseSolution()) _printf0_(" saving smb results\n");
106  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
107  }
108 
109  if(solution_type==SmbSolutionEnum)femmodel->RequestedDependentsx();
110 
111  /*Free ressources:*/
112  if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
113 
114  /*End profiler*/
115  femmodel->profiler->Stop(SMBCORE);
116 }

◆ bmb_core()

void bmb_core ( FemModel *  femmodel)

Definition at line 12 of file bmb_core.cpp.

12  {
13 
14  /*First, get BMB model from parameters*/
15  int basalforcing_model;
16  bool isplume = false;
17  femmodel->parameters->FindParam(&basalforcing_model,BasalforcingsEnum);
18 
19  if(VerboseSolution()) _printf0_(" computing basal mass balance\n");
20 
21  /*In some cases we need to run additional analyses to get the required input data*/
22  if(basalforcing_model==BasalforcingsPicoEnum){
23  femmodel->parameters->FindParam(&isplume,BasalforcingsPicoIsplumeEnum);
24  if(isplume){
25  femmodel->SetCurrentConfiguration(L2ProjectionBaseAnalysisEnum);
26  femmodel->parameters->SetParam(BaseSlopeXEnum,InputToL2ProjectEnum);
27  solutionsequence_linear(femmodel);
28  femmodel->parameters->SetParam(BaseSlopeYEnum,InputToL2ProjectEnum);
29  solutionsequence_linear(femmodel);
30  femmodel->SetCurrentConfiguration(GLheightadvectionAnalysisEnum);
31  solutionsequence_linear(femmodel);
32  }
33  }
34 
35  /*Call module now*/
36  FloatingiceMeltingRatex(femmodel);
37 
38  /*Extrude basal melt*/
39  femmodel->parameters->SetParam(BasalforcingsFloatingiceMeltingRateEnum,InputToExtrudeEnum);
40  extrudefrombase_core(femmodel);
41 }

◆ damage_core()

void damage_core ( FemModel *  femmodel)

Definition at line 12 of file damage_core.cpp.

12  {
13 
14  /*Start profiler*/
15  femmodel->profiler->Start(DAMAGECORE);
16 
17  /*intermediary*/
18  bool save_results;
19  bool dakota_analysis = false;
20  int solution_type,stabilization;
21  int numoutputs = 0;
22  char **requested_outputs = NULL;
23 
24  //first recover parameters common to all solutions
25  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
26  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
27  femmodel->parameters->FindParam(&numoutputs,DamageEvolutionNumRequestedOutputsEnum);
28  if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,DamageEvolutionRequestedOutputsEnum);
29  femmodel->parameters->FindParam(&stabilization,DamageStabilizationEnum);
30 
31  if(VerboseSolution()) _printf0_(" computing damage\n");
32  Damagex(femmodel); /* optionally calculate damage analytically first */
33  femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum);
34  if(stabilization==4){
35  solutionsequence_fct(femmodel);
36  }
37  else{
38  solutionsequence_linear(femmodel);
39  }
40 
41  if(save_results){
42  if(VerboseSolution()) _printf0_(" saving damage results\n");
43  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
44  }
45 
46  /*Free resources:*/
47  if(numoutputs){
48  for(int i=0;i<numoutputs;i++){
49  xDelete<char>(requested_outputs[i]);
50  }
51  xDelete<char*>(requested_outputs);
52  }
53 
54  /*End profiler*/
55  femmodel->profiler->Stop(DAMAGECORE);
56 }

◆ sealevelchange_core()

void sealevelchange_core ( FemModel *  femmodel)

Definition at line 17 of file sealevelchange_core.cpp.

17  { /*{{{*/
18 
19  /*Start profiler*/
20  femmodel->profiler->Start(SLRCORE);
21 
22  /*Parameters, variables:*/
23  bool save_results;
24  bool isslr=0;
25  bool isgia=0;
26  int solution_type;
27 
28  /*Retrieve parameters:*/
29  femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
30  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
31  femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
32 
33  /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
34  if(solution_type==SealevelriseSolutionEnum){
35  isslr=1;
36  isgia=1;
37  }
38 
39  /*Should we be here?:*/
40  if(!isslr)return;
41 
42  /*Verbose: */
43  if(VerboseSolution()) _printf0_(" computing sea level rise\n");
44 
45  /*Run gia core: */
46  if(isgia){
47  #ifdef _HAVE_GIA_
48  gia_core(femmodel);
49  #else
50  _error_("ISSM was not compiled with gia capabilities. Exiting");
51  #endif
52  }
53 
54  /*set SLR configuration: */
55  femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
56 
57  /*run geometry core: */
58  sealevelrise_core_geometry(femmodel);
59 
60  /*Run geodetic:*/
61  grd_core(femmodel);
62 
63  /*Run steric core for sure:*/
64  dynstr_core(femmodel);
65 
66  /*Save results: */
67  if(save_results){
68  int numoutputs = 0;
69  char **requested_outputs = NULL;
70 
71  if(VerboseSolution()) _printf0_(" saving results\n");
72  femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
73  femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
74  if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
75  }
76 
77  /*requested dependents: */
78  if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
79 
80  /*End profiler*/
81  femmodel->profiler->Stop(SLRCORE);
82 }

◆ grd_core()

void grd_core ( FemModel *  femmodel)

Definition at line 84 of file sealevelchange_core.cpp.

84  { /*{{{*/
85 
86  /*Gravity rotation deformation core GRD: */
87 
88  /*variables:*/
89  Vector<IssmDouble> *RSLg = NULL;
90  Vector<IssmDouble> *BPg = NULL;
91  Vector<IssmDouble> *RSLg_rate = NULL;
92  Vector<IssmDouble> *RSLg_eustatic = NULL;
93  Vector<IssmDouble> *U_esa = NULL;
94  Vector<IssmDouble> *U_esa_rate = NULL;
95  Vector<IssmDouble> *N_esa = NULL;
96  Vector<IssmDouble> *N_esa_rate = NULL;
97  Vector<IssmDouble> *U_north_esa = NULL;
98  Vector<IssmDouble> *U_east_esa = NULL;
99  Vector<IssmDouble> *N_gia= NULL;
100  Vector<IssmDouble> *U_gia= NULL;
101  Vector<IssmDouble> *N_gia_rate= NULL;
102  Vector<IssmDouble> *U_gia_rate= NULL;
103  SealevelMasks* masks=NULL;
104 
105  /*parameters:*/
106  bool iscoupler;
107  int solution_type;
108  int modelid,earthid;
109  bool istransientmasstransport;
110  int frequency,count;
111  int horiz;
112  int geodetic=0;
113  IssmDouble dt;
114  IssmDouble oceanarea;
115  int bp_compute_fingerprints=0;
116 
117  /*Should we even be here?:*/
118  femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum); if(!geodetic)return;
119 
120  /*Verbose: */
121  if(VerboseSolution()) _printf0_(" computing geodetic sea level rise\n");
122 
123  /*retrieve more parameters:*/
124  femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
125  femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
126  femmodel->parameters->FindParam(&count,SealevelriseRunCountEnum);
127  femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
128  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
129  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
130 
131  if(iscoupler){
132  femmodel->parameters->FindParam(&modelid,ModelIdEnum);
133  femmodel->parameters->FindParam(&earthid,EarthIdEnum);
134  femmodel->parameters->FindParam(&istransientmasstransport,TransientIsmasstransportEnum);
135  }
136  else{
137  /* we are here, we are not running in a coupler, so we will indeed compute SLR,
138  * so make sure we are identified as being the Earth.:*/
139  modelid=1; earthid=1;
140  /* in addition, if we are running solution_type SealevelriseSolutionEnum, make sure we
141  * run, irresepective of the time settings:*/
142  count=frequency;
143  }
144 
145  /*If we are running in coupled mode, the Earth model needs to run its own mass transport (if
146  * not already done by the mass trasnport module. For ice caps, they rely on the transient mass
147  * transport module exclusively:*/
148  if(iscoupler) if(modelid==earthid) if(!istransientmasstransport) EarthMassTransport(femmodel);
149 
150  /*increment counter, or call solution core if count==frequency:*/
151  if (count<frequency){
152  count++; femmodel->parameters->SetParam(count,SealevelriseRunCountEnum);
153  return;
154  }
155 
156  /*call sea-level rise sub cores:*/
157  if(iscoupler){
158  /*transfer cumulated deltathickness forcing from ice caps to earth model: */
159  TransferForcing(femmodel,SealevelriseCumDeltathicknessEnum);
160 
161  /*we have accumulated thicknesses, dump them in deltathcikness: */
162  if(modelid==earthid)InputDuplicatex(femmodel,SealevelriseCumDeltathicknessEnum,SurfaceloadIceThicknessChangeEnum);
163  }
164 
165  /*run cores:*/
166  if(modelid==earthid){
167 
168  /*call masks core: */
169  masks=sealevelrise_core_masks(femmodel);
170 
171  /*call eustatic core (generalized eustatic - Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS) */
172  RSLg_eustatic=sealevelrise_core_eustatic(femmodel,masks,&oceanarea);
173 
174  /*call non-eustatic core (ocean loading tems - 2nd and 5th terms on the RHS of Farrel and Clark) */
175  RSLg=sealevelrise_core_noneustatic(femmodel,masks,RSLg_eustatic,oceanarea);
176 
177  /*compute other elastic geodetic signatures, such as components of 3-D crustal motion: */
178  sealevelrise_core_elastic(&U_esa,&U_north_esa,&U_east_esa,femmodel,RSLg,masks);
179 
180  /*compute viscosus (GIA) geodetic signatures:*/
181  sealevelrise_core_viscous(&U_gia,&N_gia,femmodel,RSLg);
182 
183  /*compute sea-level rise (low-order spherical harmonics coefficients) diagnostics:*/
184  sealevelrise_diagnostics(femmodel,RSLg);
185 
186  /*recover N_esa = U_esa + RSLg:*/
187  N_esa=U_esa->Duplicate(); U_esa->Copy(N_esa); N_esa->AXPY(RSLg,1);
188 
189  /*if we had bottom pressure loading, remove dynamic sea level from geoid: */
190  femmodel->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
191  if(bp_compute_fingerprints){
192  GetVectorFromInputsx(&BPg,femmodel,DslSeaWaterPressureChangeAtSeaFloor,VertexSIdEnum);
193  N_esa->AXPY(BPg,-1);
194  }
195 
196  /*transform these values into rates (as we only run this once each frequency turn:*/
197  N_esa_rate=N_esa->Duplicate(); N_esa->Copy(N_esa_rate); N_esa_rate->Scale(1/(dt*frequency));
198  U_esa_rate=U_esa->Duplicate(); U_esa->Copy(U_esa_rate); U_esa_rate->Scale(1/(dt*frequency));
199  N_gia_rate=N_gia->Duplicate(); N_gia->Copy(N_gia_rate); N_gia_rate->Scale(1/(dt*frequency));
200  U_gia_rate=U_gia->Duplicate(); U_gia->Copy(U_gia_rate); U_gia_rate->Scale(1/(dt*frequency));
201  RSLg_rate=RSLg->Duplicate(); RSLg->Copy(RSLg_rate); RSLg_rate->Scale(1/(dt*frequency));
202 
203  /*get some results into elements:{{{*/
204  InputUpdateFromVectorx(femmodel,U_esa_rate,SealevelUEsaRateEnum,VertexSIdEnum);
205  InputUpdateFromVectorx(femmodel,N_esa_rate,SealevelNEsaRateEnum,VertexSIdEnum);
206  InputUpdateFromVectorx(femmodel,U_gia_rate,UGiaRateEnum,VertexSIdEnum);
207  InputUpdateFromVectorx(femmodel,N_gia_rate,NGiaRateEnum,VertexSIdEnum);
208  InputUpdateFromVectorx(femmodel,RSLg_rate,SealevelRSLRateEnum,VertexSIdEnum);
209  InputUpdateFromVectorx(femmodel,U_esa,SealevelUEsaEnum,VertexSIdEnum);
210  InputUpdateFromVectorx(femmodel,N_esa,SealevelNEsaEnum,VertexSIdEnum);
211  InputUpdateFromVectorx(femmodel,U_gia,UGiaEnum,VertexSIdEnum);
212  InputUpdateFromVectorx(femmodel,N_gia,NGiaEnum,VertexSIdEnum);
213  InputUpdateFromVectorx(femmodel,RSLg,SealevelRSLEnum,VertexSIdEnum);
214  InputUpdateFromVectorx(femmodel,RSLg_eustatic,SealevelRSLEustaticEnum,VertexSIdEnum);
215 
216  if (horiz){
217  InputUpdateFromVectorx(femmodel,U_north_esa,SealevelUNorthEsaEnum,VertexSIdEnum); // north motion
218  InputUpdateFromVectorx(femmodel,U_east_esa,SealevelUEastEsaEnum,VertexSIdEnum); // east motion
219  } /*}}}*/
220  }
221 
222  if(iscoupler){
223  /*transfer sea level back to ice caps:*/
224  TransferSealevel(femmodel,SealevelNEsaRateEnum);
225  TransferSealevel(femmodel,NGiaRateEnum);
226  TransferSealevel(femmodel,SealevelUEsaRateEnum);
227  TransferSealevel(femmodel,UGiaRateEnum);
228 
229  //reset cumdeltathickness to 0:
230  InputUpdateFromConstantx(femmodel->inputs2,femmodel->elements,0.,SealevelriseCumDeltathicknessEnum);
231  }
232 
233  /*reset counter to 1:*/
234  femmodel->parameters->SetParam(1,SealevelriseRunCountEnum); //reset counter.
235 
236  /*free ressources:{{{*/
237  delete RSLg;
238  delete RSLg_rate;
239  delete RSLg_eustatic;
240  delete U_esa;
241  delete U_esa_rate;
242  delete N_esa;
243  delete N_esa_rate;
244  delete BPg;
245 
246  if(horiz){
247  delete U_north_esa;
248  delete U_east_esa;
249  }
250  delete N_gia;
251  delete U_gia;
252  delete N_gia_rate;
253  delete U_gia_rate;
254  //delete masks;
255  /*}}}*/
256 
257 }

◆ dynstr_core()

void dynstr_core ( FemModel *  femmodel)

Definition at line 259 of file sealevelchange_core.cpp.

259  { /*{{{*/
260 
261  /*variables:*/
262  Vector<IssmDouble> *bedrock = NULL;
263  Vector<IssmDouble> *SL = NULL;
264  Vector<IssmDouble> *steric_rate_g = NULL;
265  Vector<IssmDouble> *dynamic_rate_g = NULL;
266  Vector<IssmDouble> *U_esa_rate= NULL;
267  Vector<IssmDouble> *N_esa_rate= NULL;
268  Vector<IssmDouble> *U_gia_rate= NULL;
269  Vector<IssmDouble> *N_gia_rate= NULL;
270 
271  /*parameters: */
272  bool isslr=0;
273  int solution_type;
274  IssmDouble dt;
275  int geodetic=0;
276 
277  /*Retrieve parameters:*/
278  femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
279  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
280  femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
281  femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum);
282 
283  /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/
284  if(solution_type==SealevelriseSolutionEnum)isslr=1;
285 
286  /*Should we be here?:*/
287  if(!isslr)return;
288 
289  /*Verbose: */
290  if(VerboseSolution()) _printf0_(" computing steric sea level rise\n");
291 
292  /*Retrieve geoid viscous and elastic rates, bedrock uplift viscous and elastic rates + steric rate, as vectors:*/
293  GetVectorFromInputsx(&bedrock,femmodel,BedEnum,VertexSIdEnum);
294  GetVectorFromInputsx(&SL,femmodel,SealevelEnum,VertexSIdEnum);
295  GetStericRate(&steric_rate_g,femmodel);
296  GetDynamicRate(&dynamic_rate_g,femmodel);
297  if(geodetic){
298  GetVectorFromInputsx(&U_esa_rate,femmodel,SealevelUEsaRateEnum,VertexSIdEnum);
299  GetVectorFromInputsx(&U_gia_rate,femmodel,UGiaRateEnum,VertexSIdEnum);
300  GetVectorFromInputsx(&N_esa_rate,femmodel,SealevelNEsaRateEnum,VertexSIdEnum);
301  GetVectorFromInputsx(&N_gia_rate,femmodel,NGiaRateEnum,VertexSIdEnum);
302  }
303 
304  /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate dt*/
305  if(geodetic){
306  SL->AXPY(N_gia_rate,dt);
307  SL->AXPY(N_esa_rate,dt);
308  }
309  SL->AXPY(steric_rate_g,dt);
310  SL->AXPY(dynamic_rate_g,dt);
311 
312  /*compute new bedrock position: */
313  if(geodetic){
314  bedrock->AXPY(U_esa_rate,dt);
315  bedrock->AXPY(U_gia_rate,dt);
316  }
317 
318  /*update element inputs:*/
319  InputUpdateFromVectorx(femmodel,bedrock,BedEnum,VertexSIdEnum);
320  InputUpdateFromVectorx(femmodel,SL,SealevelEnum,VertexSIdEnum);
321 
322  /*Free ressources:*/
323  delete bedrock;
324  delete SL;
325  delete steric_rate_g;
326  delete dynamic_rate_g;
327  if(geodetic){
328  delete U_esa_rate;
329  delete U_gia_rate;
330  delete N_esa_rate;
331  delete N_gia_rate;
332  }
333 }

◆ sealevelrise_core_geometry()

void sealevelrise_core_geometry ( FemModel *  femmodel)

Definition at line 351 of file sealevelchange_core.cpp.

351  { /*{{{*/
352 
353  /*Geometry core where we compute indices into tables pre computed in the SealevelRiseAnalysis: */
354 
355  /*parameters: */
356  bool spherical=true;
357  IssmDouble *latitude = NULL;
358  IssmDouble *longitude = NULL;
359  IssmDouble *radius = NULL;
360  IssmDouble *xx = NULL;
361  IssmDouble *yy = NULL;
362  IssmDouble *zz = NULL;
363  int horiz;
364  bool geometrydone = false;
365 
366 
367  /*retrieve parameters:*/
368  femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
369  femmodel->parameters->FindParam(&geometrydone,SealevelriseGeometryDoneEnum);
370 
371  if(geometrydone){
372  if(VerboseSolution()) _printf0_(" geometrical offsets have already been computed, skipping \n");
373  return; //don't need to run this again.
374  }
375 
376  /*Verbose: */
377  if(VerboseSolution()) _printf0_(" computing geometrical offsets into precomputed Green tables \n");
378 
379  /*first, recover lat,long and radius vectors from vertices: */
380  VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
381  if(horiz) VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
382 
383 
384  /*Run sealevelrie geometry routine in elements:*/
385  for(int i=0;i<femmodel->elements->Size();i++){
386  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
387  element->SealevelriseGeometry(latitude,longitude,radius,xx,yy,zz);
388  }
389 
390  /*Free ressources:*/
391  if(horiz){
392  xDelete<IssmDouble>(xx);
393  xDelete<IssmDouble>(yy);
394  xDelete<IssmDouble>(zz);
395  }
396  xDelete<IssmDouble>(longitude);
397  xDelete<IssmDouble>(radius);
398 
399  /*Record the fact that we ran this module already: */
400  femmodel->parameters->SetParam(true,SealevelriseGeometryDoneEnum);
401 
402 
403 }/*}}}*/

◆ sealevelrise_core_masks()

SealevelMasks* sealevelrise_core_masks ( FemModel *  femmodel)

Definition at line 336 of file sealevelchange_core.cpp.

336  { /*{{{*/
337 
338  if(VerboseSolution()) _printf0_(" computing sea level masks\n");
339 
340  /*initialize SealevelMasks structure: */
341  SealevelMasks* masks=new SealevelMasks(femmodel->elements->Size());
342 
343  /*go through elements and fill the masks: */
344  for (int i=0;i<femmodel->elements->Size();i++){
345  Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
346  element->SetSealevelMasks(masks);
347  }
348 
349  return masks;
350 }/*}}}*/

◆ sealevelrise_core_eustatic()

Vector<IssmDouble>* sealevelrise_core_eustatic ( FemModel *  femmodel,
SealevelMasks *  mask,
IssmDouble *  poceanarea 
)

Definition at line 404 of file sealevelchange_core.cpp.

404  { /*{{{*/
405 
406  /*Eustatic core of the SLR solution (terms that are constant with respect to sea-level)*/
407 
408  Vector<IssmDouble> *RSLgi = NULL;
409  IssmDouble RSLgi_oceanaverage = 0;
410 
411  /*parameters: */
412  int gsize;
413  IssmDouble oceanarea;
414 
415  /*outputs:*/
416  IssmDouble eustatic;
417 
418  if(VerboseSolution()) _printf0_(" computing eustatic components on ice\n");
419 
420 
421  /*Figure out size of g-set deflection vector and allocate solution vector: */
422  gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
423 
424  /*Initialize:*/
425  RSLgi = new Vector<IssmDouble>(gsize);
426 
427  /*call the eustatic main module: */
428  femmodel->SealevelriseEustatic(RSLgi,&oceanarea,&eustatic, masks); //this computes
429 
430  /*we need to average RSLgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */
431  RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,masks, oceanarea);
432 
433  /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the
434  * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/
435  RSLgi->Shift(-eustatic-RSLgi_oceanaverage);
436 
437  /*save eustatic value for results: */
438  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic));
439 
440  /*Assign output pointers and return: */
441  *poceanarea=oceanarea;
442  return RSLgi;
443 }/*}}}*/

◆ sealevelrise_core_noneustatic()

Vector<IssmDouble>* sealevelrise_core_noneustatic ( FemModel *  femmodel,
SealevelMasks *  masks,
Vector< IssmDouble > *  RSLg_eustatic,
IssmDouble  oceanarea 
)

Definition at line 444 of file sealevelchange_core.cpp.

444  { /*{{{*/
445 
446  /*sealevelrise_core_noneustatic.cpp //this computes the contributions from Eq.4 of Farrel and Clarke, rhs terms 2 and 5.
447  non eustatic core of the SLR solution */
448 
449  Vector<IssmDouble> *RSLg = NULL;
450  Vector<IssmDouble> *RSLg_old = NULL;
451  Vector<IssmDouble> *BPg = NULL;
452 
453  Vector<IssmDouble> *RSLgo = NULL; //ocean convolution of the perturbation to gravity potential.
454  Vector<IssmDouble> *RSLgo_rot= NULL; // rotational feedback
455  IssmDouble RSLgo_oceanaverage = 0; //average of RSLgo over the ocean.
456 
457  /*parameters: */
458  int count;
459  bool save_results;
460  int gsize;
461  bool converged=true;
462  bool rotation=true;
463  bool verboseconvolution=true;
464  int max_nonlinear_iterations;
465  IssmDouble eps_rel;
466  IssmDouble eps_abs;
467  IssmDouble eustatic;
468  IssmDouble Ixz, Iyz, Izz;
469  int bp_compute_fingerprints= 0;
470 
471  if(VerboseSolution()) _printf0_(" converging on ocean components\n");
472 
473  /*Recover some parameters: */
474  femmodel->parameters->FindParam(&max_nonlinear_iterations,SolidearthSettingsMaxiterEnum);
475  femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum);
476  femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum);
477 
478  /*computational flag: */
479  femmodel->parameters->FindParam(&rotation,SolidearthSettingsRotationEnum);
480 
481  /*Figure out size of g-set deflection vector and allocate solution vector: */
482  gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
483 
484  /*Initialize:*/
485  RSLg = new Vector<IssmDouble>(gsize);
486  RSLg->Assemble();
487  RSLg_eustatic->Copy(RSLg); //first initialize RSLg with the eustatic component computed in sealevelrise_core_eustatic.
488 
489  RSLg_old = new Vector<IssmDouble>(gsize);
490  RSLg_old->Assemble();
491 
492  count=1;
493  converged=false;
494 
495  /*Start loop: */
496  for(;;){
497 
498  //save pointer to old sea level rise
499  delete RSLg_old; RSLg_old=RSLg;
500 
501  /*Initialize solution vector: */
502  RSLg = new Vector<IssmDouble>(gsize); RSLg->Assemble();
503  RSLgo = new Vector<IssmDouble>(gsize); RSLgo->Assemble();
504 
505  /*call the non eustatic module: */
506  femmodel->SealevelriseNonEustatic(RSLgo, RSLg_old, masks, verboseconvolution);
507 
508  /*assemble solution vector: */
509  RSLgo->Assemble();
510 
511  if(rotation){
512 
513  /*call rotational feedback module: */
514  RSLgo_rot = new Vector<IssmDouble>(gsize); RSLgo_rot->Assemble();
515  femmodel->SealevelriseRotationalFeedback(RSLgo_rot,RSLg_old,&Ixz,&Iyz,&Izz, masks);
516  RSLgo_rot->Assemble();
517 
518  /*save changes in inertia tensor as results: */
519  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorXZEnum,Ixz));
520  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorYZEnum,Iyz));
521  femmodel->results->AddResult(new GenericExternalResult<IssmDouble>(femmodel->results->Size()+1,SealevelInertiaTensorZZEnum,Izz));
522 
523  RSLgo->AXPY(RSLgo_rot,1);
524  }
525 
526  /*we need to average RSLgo over the ocean: RHS term 5 in Eq.4 of Farrel and clarke. Only the elements can do that: */
527  RSLgo_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgo,masks, oceanarea);
528 
529  /*RSLg is the sum of the eustatic term, and the ocean terms: */
530  RSLg_eustatic->Copy(RSLg); RSLg->AXPY(RSLgo,1);
531  RSLg->Shift(-RSLgo_oceanaverage);
532 
533  /*convergence criterion:*/
534  slrconvergence(&converged,RSLg,RSLg_old,eps_rel,eps_abs);
535 
536 
537  /*free ressources: */
538  delete RSLgo;
539  delete RSLgo_rot;
540 
541  /*Increase count: */
542  count++;
543  if(converged==true){
544  break;
545  }
546  if(count>=max_nonlinear_iterations){
547  _printf0_(" maximum number of nonlinear iterations (" << max_nonlinear_iterations << ") exceeded\n");
548  converged=true;
549  break;
550  }
551 
552  /*some minor verbosing adjustment:*/
553  if(count>1)verboseconvolution=false;
554 
555  }
556  if(VerboseConvergence()) _printf0_("\n total number of iterations: " << count-1 << "\n");
557 
558 
559  /*if we had bottom pressure loading, add dynamic sea level
560  * to RSL:*/
561  femmodel->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum);
562  if(bp_compute_fingerprints){
563  GetVectorFromInputsx(&BPg,femmodel,DslSeaWaterPressureChangeAtSeaFloor,VertexSIdEnum);
564  RSLg->AXPY(BPg,1);
565  }
566 
567  delete RSLg_old;
568  delete BPg;
569 
570  return RSLg;
571 } /*}}}*/

◆ sealevelrise_core_elastic()

void sealevelrise_core_elastic ( Vector< IssmDouble > **  pU_radial,
Vector< IssmDouble > **  pU_north,
Vector< IssmDouble > **  pU_east,
FemModel *  femmodel,
Vector< IssmDouble > *  RSLg,
SealevelMasks *  masks 
)

Definition at line 572 of file sealevelchange_core.cpp.

572  { /*{{{*/
573 
574  Vector<IssmDouble> *U_esa = NULL;
575  Vector<IssmDouble> *U_north_esa = NULL;
576  Vector<IssmDouble> *U_east_esa = NULL;
577 
578  /*parameters: */
579  int gsize;
580  bool spherical=true;
581 
582  IssmDouble *latitude = NULL;
583  IssmDouble *longitude = NULL;
584  IssmDouble *radius = NULL;
585  IssmDouble *xx = NULL;
586  IssmDouble *yy = NULL;
587  IssmDouble *zz = NULL;
588  int horiz;
589 
590  if(VerboseSolution()) _printf0_(" computing vertical and horizontal geodetic signatures\n");
591 
592  /*retrieve some parameters:*/
593  femmodel->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum);
594 
595  /*find size of vectors:*/
596  gsize = femmodel->nodes->NumberOfDofs(GsetEnum);
597 
598  /*intialize vectors:*/
599  U_esa = new Vector<IssmDouble>(gsize);
600  if (horiz){
601  U_north_esa = new Vector<IssmDouble>(gsize);
602  U_east_esa = new Vector<IssmDouble>(gsize);
603  }
604 
605  /*retrieve geometric information: */
606  VertexCoordinatesx(&latitude,&longitude,&radius,femmodel->vertices,spherical);
607  VertexCoordinatesx(&xx,&yy,&zz,femmodel->vertices);
608 
609  /*call the elastic main modlule:*/
610  femmodel->SealevelriseElastic(U_esa,U_north_esa,U_east_esa,RSLg, masks);
611 
612  /*Assign output pointers:*/
613  *pU_esa=U_esa;
614  if(horiz){
615  *pU_east_esa=U_east_esa;
616  *pU_north_esa=U_north_esa;
617  }
618 
619  /*Free ressources: */
620  xDelete<IssmDouble>(longitude);
621  xDelete<IssmDouble>(latitude);
622  xDelete<IssmDouble>(xx);
623  xDelete<IssmDouble>(yy);
624  xDelete<IssmDouble>(zz);
625  xDelete<IssmDouble>(radius);
626 }

◆ sealevelrise_core_viscous()

void sealevelrise_core_viscous ( Vector< IssmDouble > **  pU_gia,
Vector< IssmDouble > **  pN_gia,
FemModel *  femmodel,
Vector< IssmDouble > *  RSLg 
)

Definition at line 628 of file sealevelchange_core.cpp.

628  { /*{{{*/
629 
630  /*variables:*/
631  Vector<IssmDouble> *U_gia = NULL;
632  Vector<IssmDouble> *N_gia = NULL;
633 
634  /*parameters:*/
635  int frequency;
636  IssmDouble dt;
637 
638  if(VerboseSolution()) _printf0_(" computing viscous components\n");
639 
640  /*retrieve some parameters:*/
641  femmodel->parameters->FindParam(&frequency,SolidearthSettingsRunFrequencyEnum);
642  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
643 
644  /*recover GIA rates:*/
645  GetVectorFromInputsx(&U_gia,femmodel,UGiaRateEnum,VertexSIdEnum);
646  GetVectorFromInputsx(&N_gia,femmodel,NGiaRateEnum,VertexSIdEnum);
647 
648  /*we just loaded rates, that's not what's being asked, scale by time:*/
649  U_gia->Scale(frequency*dt);
650  N_gia->Scale(frequency*dt);
651 
652  /*Assign output pointers:*/
653  *pU_gia=U_gia;
654  *pN_gia=N_gia;
655 
656 }

◆ sealevelrise_diagnostics()

void sealevelrise_diagnostics ( FemModel *  femmodel,
Vector< IssmDouble > *  RSLg 
)

Definition at line 658 of file sealevelchange_core.cpp.

658  { /*{{{*/
659 
660  /*compute spherical harmonics deg 1 and deg 2 coefficeints:*/
661 
662 }

◆ objectivefunction()

IssmDouble objectivefunction ( IssmDouble  search_scalar,
FemModel *  femmodel 
)

◆ GetStericRate()

void GetStericRate ( Vector< IssmDouble > **  psteric_rate_g,
FemModel *  femmodel 
)

Definition at line 707 of file sealevelchange_core.cpp.

707  { /*{{{*/
708 
709  int dslmodel=-1;
710  IssmDouble time;
711 
712  /*variables:*/
713  Vector<IssmDouble> *steric_rate_g = NULL;
714 
715  /*Update steric rates before retrieving them on Vertex SID set:*/
716  femmodel->parameters->FindParam(&dslmodel,DslModelEnum);
717  femmodel->parameters->FindParam(&time,TimeEnum);
718  if(dslmodel==1){
719  TransientInput2* transient_input = femmodel->inputs2->GetTransientInput(DslGlobalAverageThermostericSeaLevelChangeEnum);
720  TriaInput2* tria_input=transient_input->GetTriaInput(time);
721  Input2* tria_input_copy=tria_input->copy();
722  tria_input_copy->ChangeEnum(DslStericRateEnum);
723  femmodel->inputs2->AddInput(tria_input_copy);
724  }
725  else if (dslmodel==2){
726  int modelid;
727 
728  /*Recover modelid:*/
729  femmodel->parameters->FindParam(&modelid,DslModelidEnum);
730 
731  modelid--; //from matlab.
732 
733  /*find the DslGlobalAverageThermostericSeaLevelChangeEnum dataset of transient inputs:*/
734  DatasetInput2* dataset_input=femmodel->inputs2->GetDatasetInput2(DslGlobalAverageThermostericSeaLevelChangeEnum);
735 
736  /*Go find the modelid'th transient input:*/
737  TriaInput2* tria_input=dataset_input->GetTriaInputByOffset(modelid);
738 
739  /*Plug into DslStericRate input: */
740  Input2* tria_input_copy=tria_input->copy();
741  tria_input_copy->ChangeEnum(DslStericRateEnum);
742  femmodel->inputs2->AddInput(tria_input_copy);
743  }
744  else _error_("not implemented yet");
745 
746  GetVectorFromInputsx(&steric_rate_g,femmodel,DslStericRateEnum,VertexSIdEnum);
747  *psteric_rate_g=steric_rate_g;
748 }

◆ GetDynamicRate()

void GetDynamicRate ( Vector< IssmDouble > **  pdynamic_rate_g,
FemModel *  femmodel 
)

Definition at line 664 of file sealevelchange_core.cpp.

664  { /*{{{*/
665 
666  int dslmodel=-1;
667  IssmDouble time;
668 
669  /*variables:*/
670  Vector<IssmDouble> *dynamic_rate_g = NULL;
671 
672  /*Update steric rates before retrieving them on Vertex SID set:*/
673  femmodel->parameters->FindParam(&dslmodel,DslModelEnum);
674  femmodel->parameters->FindParam(&time,TimeEnum);
675  if(dslmodel==1){
676  TransientInput2* transient_input = femmodel->inputs2->GetTransientInput(DslSeaSurfaceHeightChangeAboveGeoidEnum);
677  TriaInput2* tria_input=transient_input->GetTriaInput(time);
678  Input2* tria_input_copy=tria_input->copy();
679  tria_input_copy->ChangeEnum(DslDynamicRateEnum);
680  femmodel->inputs2->AddInput(tria_input_copy);
681  }
682  else if(dslmodel==2){
683 
684  int modelid;
685 
686  /*Recover modelid:*/
687  femmodel->parameters->FindParam(&modelid,DslModelidEnum);
688  modelid--; //from matlab.
689 
690  /*find the DslSeaSurfaceHeightChangeAboveGeoidEnum dataset of transient inputs:*/
691  DatasetInput2* dataset_input=femmodel->inputs2->GetDatasetInput2(DslSeaSurfaceHeightChangeAboveGeoidEnum);
692 
693  /*Go find the modelid'th transient input:*/
694  TriaInput2* tria_input=dataset_input->GetTriaInputByOffset(modelid);
695 
696  /*Plug into DslDynamicRate input: */
697  Input2* tria_input_copy=tria_input->copy();
698  tria_input_copy->ChangeEnum(DslDynamicRateEnum);
699  femmodel->inputs2->AddInput(tria_input_copy);
700  }
701  else _error_("not implemented yet");
702 
703  GetVectorFromInputsx(&dynamic_rate_g,femmodel,DslDynamicRateEnum,VertexSIdEnum);
704  *pdynamic_rate_g=dynamic_rate_g;
705 }

◆ GradJSearch()

int GradJSearch ( IssmDouble *  search_vector,
FemModel *  femmodel,
int  step 
)

◆ ProcessArguments()

void ProcessArguments ( int *  solution,
char **  pbinname,
char **  poutbinname,
char **  ptoolkitsname,
char **  plockname,
char **  prestartname,
char **  prootpath,
int  argc,
char **  argv 
)

Definition at line 10 of file ProcessArguments.cpp.

10  {
11 
12  /*Check input arguments*/
13  if(argc<2) _error_("Usage error: no solution requested");
14  if(argc<3) _error_("Usage error: missing execution directory");
15  if(argc<4) _error_("Usage error: missing model name");
16 
17  /*Get some arguments*/
18  *solution_type = StringToEnumx(argv[1]);
19  char* rootpatharg = argv[2];
20  char* modelname = argv[3];
21 
22  /*Recover myrank and length of string "my_rank" */
23  int my_rank = IssmComm::GetRank();
24  int rank_length = (my_rank == 0 ? 1 : (int)(log10(static_cast<double>(my_rank))+1));
25 
26  /*Create rootpath from argument*/
27  char* rootpath = xNew<char>(strlen(rootpatharg)+2); sprintf(rootpath,"%s/",rootpatharg);
28 
29  /*Create all file paths*/
30  char* binfilename = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".bin") +1); sprintf(binfilename, "%s%s%s",rootpath,modelname,".bin");
31  char* outbinfilename = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".outbin") +1); sprintf(outbinfilename,"%s%s%s",rootpath,modelname,".outbin");
32  char* toolkitsfilename = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".toolkits") +1); sprintf(toolkitsfilename,"%s%s%s",rootpath,modelname,".toolkits");
33  char* lockfilename = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".lock") +1); sprintf(lockfilename, "%s%s%s",rootpath,modelname,".lock");
34  char* restartfilename = xNew<char>(strlen(rootpath)+strlen(modelname)+strlen(".rst.")+rank_length+1);
35  sprintf(restartfilename,"%s%s%s%i",rootpath,modelname,".rst.",my_rank);
36 
37  /*Clean up and assign output pointer*/
38  *pbinfilename=binfilename;
39  *poutbinfilename=outbinfilename;
40  *ptoolkitsfilename=toolkitsfilename;
41  *plockfilename=lockfilename;
42  *prestartfilename=restartfilename;
43  *prootpath=rootpath;
44 
45 }

◆ WriteLockFile()

void WriteLockFile ( char *  filename)

Definition at line 10 of file WriteLockFile.cpp.

10  {
11 
12  /*recover my_rank:*/
13  int my_rank=IssmComm::GetRank();
14 
15  /* output: */
16  FILE* fid=NULL;
17 
18  /* Open lock file and write 1 into it: */
19  if(my_rank==0){
20  fid=fopen(filename,"w");
21  if(fid==NULL) _error_("error message: could not open lock file " << filename);
22 
23  /*Close file: */
24  if(fclose(fid)!=0) _error_("could not close lock file " << filename);
25  }
26 
27 }

◆ ResetBoundaryConditions()

void ResetBoundaryConditions ( FemModel *  femmodel,
int  analysis_type 
)

Definition at line 9 of file ResetBoundaryConditions.cpp.

9  {
10 
11  /*variables: */
12  Vector<IssmDouble>* yg = NULL;
13 
14  if(VerboseSolution()) _printf0_(" updating boundary conditions...\n");
15  _assert_(femmodel->analysis_type_list[femmodel->analysis_counter]==analysis_type);
16 
17  /*set current analysis: */
18  femmodel->SetCurrentConfiguration(analysis_type);
19  int index = femmodel->AnalysisIndex(analysis_type);
20 
21  /*retrieve boundary conditions from element inputs :*/
22  GetSolutionFromInputsx(&yg,femmodel);
23 
24  /*update spcs using this new vector of constraints: */
25  UpdateDynamicConstraintsx(femmodel->constraints,femmodel->nodes,femmodel->parameters,yg);
26 
27  /*Free ressources:*/
28  delete yg;
29 }

◆ PrintBanner()

void PrintBanner ( void  )

◆ TransferForcing()

void TransferForcing ( FemModel *  femmodel,
int  forcingenum 
)

Definition at line 752 of file sealevelchange_core.cpp.

752  { /*{{{*/
753 
754  /*forcing being transferred from models to earth: */
755  IssmDouble** forcings=NULL;
756  IssmDouble* forcing=NULL;
757  Vector<IssmDouble>* forcingglobal=NULL;
758  int* nvs=NULL;
759 
760  /*transition vectors:*/
761  IssmDouble** transitions=NULL;
762  int ntransitions;
763  int* transitions_m=NULL;
764  int* transitions_n=NULL;
765  int nv;
766 
767  /*communicators:*/
768  ISSM_MPI_Comm tocomm;
769  ISSM_MPI_Comm* fromcomms=NULL;
770  ISSM_MPI_Status status;
771  int my_rank;
772  int modelid,earthid;
773  int nummodels;
774 
775  /*Recover some parameters: */
776  femmodel->parameters->FindParam(&modelid,ModelIdEnum);
777  femmodel->parameters->FindParam(&earthid,EarthIdEnum);
778  femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
779  my_rank=IssmComm::GetRank();
780 
781  /*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */
782  if(modelid==earthid){
783  GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
784  if(!parcoms)_error_("TransferForcing error message: could not find IcecapToEarthComm communicator");
785  fromcomms=parcoms->GetParameterValue();
786  }
787  else {
788  GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
789  if(!parcom)_error_("TransferForcing error message: could not find IcecapToEarthComm communicator");
790  tocomm=parcom->GetParameterValue();
791  }
792 
793  /*For each icecap, retrieve the forcing vector that will be sent to the earth model: */
794  if(modelid!=earthid){
795  nv=femmodel->vertices->NumberOfVertices();
796  GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum);
797  }
798 
799  /*Send the forcing to the earth model:{{{*/
800  if(my_rank==0){
801  if(modelid==earthid){
802  forcings=xNew<IssmDouble*>(nummodels-1);
803  nvs=xNew<int>(nummodels-1);
804  for(int i=0;i<earthid;i++){
805  ISSM_MPI_Recv(nvs+i, 1, ISSM_MPI_INT, 0,i, fromcomms[i], &status);
806  forcings[i]=xNew<IssmDouble>(nvs[i]);
807  ISSM_MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status);
808  }
809 
810  }
811  else{
812  ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, modelid, tocomm);
813  ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, tocomm);
814  }
815  }
816  /*}}}*/
817 
818  /*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/
819  if(modelid==earthid){
820 
821  /*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions.
822  *First, build the global delta thickness vector in the earth model: */
823  nv=femmodel->vertices->NumberOfVertices();
824  GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
825 
826  /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/
827  femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
828 
829  if(ntransitions!=earthid)_error_("TransferForcing error message: number of transition vectors is not equal to the number of icecaps!");
830 
831  /*Go through all the delta thicknesses coming from each ice cap: */
832  if(my_rank==0){
833  for(int i=0;i<earthid;i++){
834 
835  IssmDouble* forcingfromcap= forcings[i]; //careful, this only exists on rank 0 of the earth model!
836  IssmDouble* transition=transitions[i];
837  int M=transitions_m[i];
838 
839  /*build index to plug values: */
840  int* index=xNew<int>(M); for(int i=0;i<M;i++)index[i]=reCast<int>(transition[i])-1; //matlab indexing!
841 
842  /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular
843  * ice cap: */
844  forcingglobal->SetValues(M,index,forcingfromcap,ADD_VAL);
845  xDelete<int>(index);
846  }
847  }
848 
849  /*Assemble vector:*/
850  forcingglobal->Assemble();
851 
852  /*Plug into elements:*/
853  InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum);
854  }
855  /*}}}*/
856 
857  /*Free ressources:{{{*/
858  if(forcings){
859  for(int i=0;i<nummodels-1;i++){
860  IssmDouble* temp=forcings[i]; xDelete<IssmDouble>(temp);
861  }
862  xDelete<IssmDouble*>(forcings);
863  }
864  if(forcing)xDelete<IssmDouble>(forcing);
865  if(forcingglobal)delete forcingglobal;
866  if(transitions){
867  for(int i=0;i<earthid;i++){
868  IssmDouble* temp=transitions[i];
869  xDelete<IssmDouble>(temp);
870  }
871  xDelete<IssmDouble*>(transitions);
872  xDelete<int>(transitions_m);
873  xDelete<int>(transitions_n);
874  }
875  if(nvs)xDelete<int>(nvs);
876  /*}}}*/
877 
878 } /*}}}*/

◆ TransferSealevel()

void TransferSealevel ( FemModel *  femmodel,
int  forcingenum 
)

Definition at line 879 of file sealevelchange_core.cpp.

879  { /*{{{*/
880 
881  /*forcing being transferred from earth to ice caps: */
882  IssmDouble* forcing=NULL;
883  IssmDouble* forcingglobal=NULL;
884 
885  /*transition vectors:*/
886  IssmDouble** transitions=NULL;
887  int ntransitions;
888  int* transitions_m=NULL;
889  int* transitions_n=NULL;
890  int nv;
891 
892  /*communicators:*/
893  ISSM_MPI_Comm fromcomm;
894  ISSM_MPI_Comm* tocomms=NULL;
895  ISSM_MPI_Status status;
896  int my_rank;
897  int modelid,earthid;
898  int nummodels;
899  int numcoms;
900 
901  /*Recover some parameters: */
902  femmodel->parameters->FindParam(&modelid,ModelIdEnum);
903  femmodel->parameters->FindParam(&earthid,EarthIdEnum);
904  femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
905  my_rank=IssmComm::GetRank();
906 
907  /*retrieve the inter communicators that will be used to send data from earth to ice caps:*/
908  if(modelid==earthid){
909  GenericParam<ISSM_MPI_Comm*>* parcoms = dynamic_cast<GenericParam<ISSM_MPI_Comm*>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
910  if(!parcoms)_error_("TransferSealevel error message: could not find IcecapToEarthComm communicator");
911  tocomms=parcoms->GetParameterValue();
912  //femmodel->parameters->FindParam((int**)(&tocomms),&numcoms,IcecapToEarthCommEnum);
913  }
914  else{
915  GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(IcecapToEarthCommEnum));
916  if(!parcom)_error_("TransferSealevel error message: could not find IcecapToEarthComm communicator");
917  fromcomm=parcom->GetParameterValue();
918  //femmodel->parameters->FindParam((int*)(&fromcomm), IcecapToEarthCommEnum);
919  }
920 
921  /*Retrieve sea-level on earth model: */
922  if(modelid==earthid){
923  nv=femmodel->vertices->NumberOfVertices();
924  GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
925  }
926 
927  /*Send the forcing to the ice caps:{{{*/
928  if(my_rank==0){
929 
930  if(modelid==earthid){
931 
932  /*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */
933  femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
934 
935  if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!");
936 
937  for(int i=0;i<earthid;i++){
938  nv=transitions_m[i];
939  forcing=xNew<IssmDouble>(nv);
940  IssmDouble* transition=transitions[i];
941  for(int j=0;j<nv;j++){
942  forcing[j]=forcingglobal[reCast<int>(transition[j])-1];
943  }
944  ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, i, tocomms[i]);
945  ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, i, tocomms[i]);
946  }
947  }
948  else{
949  ISSM_MPI_Recv(&nv, 1, ISSM_MPI_INT, 0, modelid, fromcomm, &status);
950  forcing=xNew<IssmDouble>(nv);
951  ISSM_MPI_Recv(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, fromcomm, &status);
952  }
953  }
954  /*}}}*/
955 
956  /*On each ice cap, spread the forcing across cpus, and update the elements dataset accordingly: {{{*/
957  if(modelid!=earthid){
958 
959  ISSM_MPI_Bcast(&nv,1,ISSM_MPI_INT,0,IssmComm::GetComm());
960  if(my_rank!=0)forcing=xNew<IssmDouble>(nv);
961  ISSM_MPI_Bcast(forcing,nv,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
962 
963  /*Plug into elements:*/
964  InputUpdateFromVectorx(femmodel,forcing,forcingenum,VertexSIdEnum);
965  }
966  /*}}}*/
967 
968  /*Free ressources:{{{*/
969  if(forcingglobal)xDelete<IssmDouble>(forcingglobal);
970  if(forcing)xDelete<IssmDouble>(forcing);
971  if(transitions){
972  for(int i=0;i<ntransitions;i++){
973  IssmDouble* temp=transitions[i];
974  xDelete<IssmDouble>(temp);
975  }
976  xDelete<IssmDouble*>(transitions);
977  xDelete<int>(transitions_m);
978  xDelete<int>(transitions_n);
979  }
980  /*}}}*/
981 
982 } /*}}}*/

◆ EarthMassTransport()

void EarthMassTransport ( FemModel *  femmodel)

Definition at line 983 of file sealevelchange_core.cpp.

983  { /*{{{*/
984 
985  IssmDouble time,dt;
986  Vector<IssmDouble> *oldthickness = NULL;
987  Vector<IssmDouble> *newthickness = NULL;
988  Vector<IssmDouble> *deltathickness = NULL;
989  Vector<IssmDouble> *cumdeltathickness = NULL;
990  int nv;
991 
992  if(VerboseSolution()) _printf0_(" computing earth mass transport\n");
993 
994  /*This module has to be recoded! We do not grab spc thicknesses from an earth model
995  * anymore! Thicknesses come from a mass transport module applied to each basin!
996  * Commeting out for now:*/
997  _error_("EarthMassTransport error message: not supported anymore!");
998 
999  /*This mass transport module for the Earth is because we might have thickness variations as spcs
1000  * specified in the md.slr class, outside of what we will get from the icecaps. That's why we get t
1001  * the thickness variations from SealevelriseSpcthicknessEnum.*/
1002 
1003  /*No mass transport module was called, so we are just going to retrieve the geometry thickness
1004  * at this time step, at prior time step, and plug the difference as deltathickness: */
1005  /*femmodel->parameters->FindParam(&time,TimeEnum);
1006  femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum);
1007  nv=femmodel->vertices->NumberOfVertices();
1008 
1009  GetVectorFromInputsx(&newthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum);
1010  GetVectorFromInputsx(&oldthickness,femmodel,SealevelriseSpcthicknessEnum,VertexSIdEnum,time-dt);*/
1011 
1012  /*compute deltathickness: */
1013  /*deltathickness = new Vector<IssmDouble>(nv);
1014  newthickness->Copy(deltathickness); deltathickness->AXPY(oldthickness,-1); */
1015 
1016  /*plug into elements:*/
1017 // InputUpdateFromVectorx(femmodel,deltathickness,SurfaceloadIceThicknessChangeEnum,VertexSIdEnum);
1018 
1019  /*add to cumulated delta thickness: */
1020  /*GetVectorFromInputsx(&cumdeltathickness,femmodel,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);
1021  cumdeltathickness->AXPY(deltathickness,1);
1022  InputUpdateFromVectorx(femmodel,cumdeltathickness,SealevelriseCumDeltathicknessEnum,VertexSIdEnum);*/
1023 
1024  /*free ressources:*/
1025  /*delete oldthickness;
1026  delete newthickness;
1027  delete deltathickness;
1028  delete cumdeltathickness;*/
1029 
1030 } /*}}}*/

◆ slrconvergence()

void slrconvergence ( bool *  pconverged,
Vector< IssmDouble > *  RSLg,
Vector< IssmDouble > *  RSLg_old,
IssmDouble  eps_rel,
IssmDouble  eps_abs 
)

Definition at line 1031 of file sealevelchange_core.cpp.

1031  { /*{{{*/
1032 
1033  bool converged=true;
1034  IssmDouble ndS,nS;
1035  Vector<IssmDouble> *dRSLg = NULL;
1036 
1037  //compute norm(du) and norm(u) if requested
1038  dRSLg=RSLg_old->Duplicate(); RSLg_old->Copy(dRSLg); dRSLg->AYPX(RSLg,-1.0);
1039  ndS=dRSLg->Norm(NORM_TWO);
1040 
1041  if (xIsNan<IssmDouble>(ndS)) _error_("convergence criterion is NaN!");
1042 
1043  if(!xIsNan<IssmDouble>(eps_rel)){
1044  nS=RSLg_old->Norm(NORM_TWO);
1045  if (xIsNan<IssmDouble>(nS)) _error_("convergence criterion is NaN!");
1046  }
1047 
1048  //clean up
1049  delete dRSLg;
1050 
1051  //print
1052  if(!xIsNan<IssmDouble>(eps_rel)){
1053  if((ndS/nS)<eps_rel){
1054  if(VerboseConvergence()) _printf0_(setw(50) << left << " convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " < " << eps_rel*100 << " %\n");
1055  }
1056  else{
1057  if(VerboseConvergence()) _printf0_(setw(50) << left << " convergence criterion: norm(dS)/norm(S)" << ndS/nS*100 << " > " << eps_rel*100 << " %\n");
1058  converged=false;
1059  }
1060  }
1061  if(!xIsNan<IssmDouble>(eps_abs)){
1062  if(ndS<eps_abs){
1063  if(VerboseConvergence()) _printf0_(setw(50) << left << " convergence criterion: norm(dS)" << ndS << " < " << eps_abs << " \n");
1064  }
1065  else{
1066  if(VerboseConvergence()) _printf0_(setw(50) << left << " convergence criterion: norm(dS)" << ndS << " > " << eps_abs << " \n");
1067  converged=false;
1068  }
1069  }
1070 
1071  /*assign output*/
1072  *pconverged=converged;
1073 
1074 } /*}}}*/

◆ CorePointerFromSolutionEnum()

void CorePointerFromSolutionEnum ( void(**)(FemModel *)  psolutioncore,
Parameters *  parameters,
int  solutiontype 
)

Definition at line 18 of file CorePointerFromSolutionEnum.cpp.

18  {
19 
20  /*output: */
21  void (*solutioncore)(FemModel*)=NULL;
22 
23  switch(solutiontype){
24 
25  case StressbalanceSolutionEnum:
26  solutioncore=&stressbalance_core;
27  break;
28  case SteadystateSolutionEnum:
29  solutioncore=&steadystate_core;
30  break;
31  case ThermalSolutionEnum:
32  solutioncore=&thermal_core;
33  break;
34  case BalancethicknessSolutionEnum:
35  solutioncore=&balancethickness_core;
36  break;
37  case Balancethickness2SolutionEnum:
38  solutioncore=&balancethickness2_core;
39  break;
40  case BalancethicknessSoftSolutionEnum:
41  solutioncore=&dummy_core;
42  break;
43  case BalancevelocitySolutionEnum:
44  solutioncore=&balancevelocity_core;
45  break;
46  case HydrologySolutionEnum:
47  solutioncore=&hydrology_core;
48  break;
49  case SurfaceSlopeSolutionEnum:
50  solutioncore=&surfaceslope_core;
51  break;
52  case BedSlopeSolutionEnum:
53  solutioncore=&bedslope_core;
54  break;
55  case TransientSolutionEnum:
56  solutioncore=&transient_core;
57  break;
58  case MasstransportSolutionEnum:
59  solutioncore=&masstransport_core;
60  break;
61  case SealevelriseSolutionEnum:
62  solutioncore=&sealevelchange_core;
63  break;
64  case EsaSolutionEnum:
65  solutioncore=&esa_core;
66  break;
67  case GiaSolutionEnum:
68  #if _HAVE_GIA_
69  solutioncore=&gia_core;
70  #else
71  _error_("ISSM not compiled with Gia capability");
72  #endif
73  break;
74  case DamageEvolutionSolutionEnum:
75  solutioncore=&damage_core;
76  break;
77  case LoveSolutionEnum:
78  #if _HAVE_LOVE_
79  solutioncore=&love_core;
80  #else
81  _error_("ISSM not compiled with Love capability");
82  #endif
83  break;
84 
85  default:
86  _error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!");
87  break;
88  }
89 
90  /*Assign output pointer:*/
91  _assert_(psolutioncore);
92  *psolutioncore=solutioncore;
93 }

◆ WrapperCorePointerFromSolutionEnum()

void WrapperCorePointerFromSolutionEnum ( void(**)(FemModel *)  psolutioncore,
Parameters *  parameters,
int  solutiontype,
bool  nodakotacore = false 
)

Definition at line 17 of file WrapperCorePointerFromSolutionEnum.cpp.

17  {
18 
19  /*output: */
20  void (*solutioncore)(FemModel*)=NULL;
21 
22  /*parameters: */
23  bool control_analysis;
24  bool dakota_analysis;
25  int inversiontype;
26 
27  /* retrieve some parameters that tell us whether wrappers are allowed, or whether we return
28  * a pure core. Wrappers can be dakota_core (which samples many solution_cores) or control_core (which
29  * carries out adjoint based inversion on a certain core: */
30  parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum);
31  parameters->FindParam(&control_analysis,InversionIscontrolEnum);
32  parameters->FindParam(&inversiontype,InversionTypeEnum);
33 
34  if(nodakotacore)dakota_analysis=false;
35 
36  if(dakota_analysis){
37  #ifdef _HAVE_DAKOTA_
38  solutioncore=dakota_core;
39  #else
40  _error_("ISSM was not compiled with dakota support, cannot carry out dakota analysis!");
41  #endif
42  }
43  else if(control_analysis){
44  switch(inversiontype){
45  case 0: solutioncore=control_core; break;
46  case 1: solutioncore=controltao_core; break;
47  case 2: solutioncore=controlm1qn3_core; break;
48  case 3: solutioncore=controlvalidation_core; break;
49  case 4: solutioncore=controladm1qn3_core; break;
50  default: _error_("control type not supported");
51  }
52  }
53  else CorePointerFromSolutionEnum(&solutioncore,parameters,solutiontype); /*This means we retrieve a core solution that is not a wrapper*/
54 
55  /*Assign output pointer:*/
56  _assert_(psolutioncore);
57  *psolutioncore=solutioncore;
58 
59 }

◆ AdjointCorePointerFromSolutionEnum()

void AdjointCorePointerFromSolutionEnum ( void(**)(FemModel *)  padjointcore,
int  solutiontype 
)

Definition at line 18 of file AdjointCorePointerFromSolutionEnum.cpp.

18  {
19 
20  /*output: */
21  void (*adjointcore)(FemModel*)=NULL;
22 
23  switch(solutiontype){
24 
25  case StressbalanceSolutionEnum:
26  adjointcore=&adjointstressbalance_core;
27  break;
28  case SteadystateSolutionEnum:
29  adjointcore=&adjointstressbalance_core;
30  break;
31  case BalancethicknessSolutionEnum:
32  adjointcore=&adjointbalancethickness_core;
33  break;
34  case Balancethickness2SolutionEnum:
35  adjointcore=&adjointbalancethickness2_core;
36  break;
37  case BalancethicknessSoftSolutionEnum:
38  adjointcore=&dummy_core;
39  break;
40  default:
41  _error_("No adjoint has been implemented for solution " << EnumToStringx(solutiontype) << " yet");
42  break;
43  }
44 
45  /*Assign output pointer:*/
46  _assert_(padjointcore);
47  *padjointcore=adjointcore;
48 
49 }
DataSet::Size
int Size()
Definition: DataSet.cpp:399
FreeSurfaceTopAnalysisEnum
@ FreeSurfaceTopAnalysisEnum
Definition: EnumDefinitions.h:1072
GetVectorFromControlInputsx
void GetVectorFromControlInputsx(Vector< IssmDouble > **pvector, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters, const char *data)
Definition: GetVectorFromControlInputsx.cpp:9
BalancethicknessAnalysisEnum
@ BalancethicknessAnalysisEnum
Definition: EnumDefinitions.h:981
controlvalidation_core
void controlvalidation_core(FemModel *femmodel)
Definition: controlvalidation_core.cpp:155
SedimentHeadOldEnum
@ SedimentHeadOldEnum
Definition: EnumDefinitions.h:699
SmbRequestedOutputsEnum
@ SmbRequestedOutputsEnum
Definition: EnumDefinitions.h:382
HydrologySolutionEnum
@ HydrologySolutionEnum
Definition: EnumDefinitions.h:1105
DslSeaWaterPressureChangeAtSeaFloor
@ DslSeaWaterPressureChangeAtSeaFloor
Definition: EnumDefinitions.h:542
DslSeaSurfaceHeightChangeAboveGeoidEnum
@ DslSeaSurfaceHeightChangeAboveGeoidEnum
Definition: EnumDefinitions.h:541
Matlitho::viscosity
IssmDouble * viscosity
Definition: Matlitho.h:20
FemModel::AnalysisIndex
int AnalysisIndex(int)
Definition: FemModel.cpp:220
DslGlobalAverageThermostericSeaLevelChangeEnum
@ DslGlobalAverageThermostericSeaLevelChangeEnum
Definition: EnumDefinitions.h:540
EffectivePressureSubstepEnum
@ EffectivePressureSubstepEnum
Definition: EnumDefinitions.h:549
depthaverage_core
void depthaverage_core(FemModel *femmodel)
Definition: depthaverage_core.cpp:12
TimesteppingFinalTimeEnum
@ TimesteppingFinalTimeEnum
Definition: EnumDefinitions.h:430
SmbRunoffSubstepEnum
@ SmbRunoffSubstepEnum
Definition: EnumDefinitions.h:773
StressbalanceAnalysis
Definition: StressbalanceAnalysis.h:11
SaveResultsEnum
@ SaveResultsEnum
Definition: EnumDefinitions.h:302
SealevelUEsaRateEnum
@ SealevelUEsaRateEnum
Definition: EnumDefinitions.h:686
movingfront_core
void movingfront_core(FemModel *femmodel)
Definition: movingfront_core.cpp:12
BaseEnum
@ BaseEnum
Definition: EnumDefinitions.h:495
EnthalpyAnalysisEnum
@ EnthalpyAnalysisEnum
Definition: EnumDefinitions.h:1052
HydrologyShaktiAnalysis::UpdateGapHeight
void UpdateGapHeight(FemModel *femmodel)
Definition: HydrologyShaktiAnalysis.cpp:457
ControlInputSetGradientx
void ControlInputSetGradientx(Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters, IssmDouble *gradient)
Definition: ControlInputSetGradientx.cpp:9
OptPars::nsize
int nsize
Definition: OptPars.h:17
SolidearthSettingsAbstolEnum
@ SolidearthSettingsAbstolEnum
Definition: EnumDefinitions.h:306
SMBd18opddEnum
@ SMBd18opddEnum
Definition: EnumDefinitions.h:1243
SmbMassBalanceEnum
@ SmbMassBalanceEnum
Definition: EnumDefinitions.h:748
SealevelUEastEsaEnum
@ SealevelUEastEsaEnum
Definition: EnumDefinitions.h:684
TransientNumRequestedOutputsEnum
@ TransientNumRequestedOutputsEnum
Definition: EnumDefinitions.h:456
SmbEnum
@ SmbEnum
Definition: EnumDefinitions.h:358
SmbAnalysisEnum
@ SmbAnalysisEnum
Definition: EnumDefinitions.h:1274
_assert_
#define _assert_(ignore)
Definition: exceptions.h:37
FemModel::vertices
Vertices * vertices
Definition: FemModel.h:49
LoveMu0Enum
@ LoveMu0Enum
Definition: EnumDefinitions.h:236
adjointstressbalance_core
void adjointstressbalance_core(FemModel *femmodel)
Definition: adjointstressbalance_core.cpp:12
Analysis::Core
virtual void Core(FemModel *femmodel)=0
Balancethickness2SolutionEnum
@ Balancethickness2SolutionEnum
Definition: EnumDefinitions.h:980
TransientIsmasstransportEnum
@ TransientIsmasstransportEnum
Definition: EnumDefinitions.h:449
GetPassiveVectorFromControlInputsx
void GetPassiveVectorFromControlInputsx(IssmPDouble **pvector, int *pN, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters, const char *data)
Definition: GetVectorFromControlInputsx.cpp:127
IssmDouble
double IssmDouble
Definition: types.h:37
OutputResultsx
void OutputResultsx(FemModel *femmodel)
Definition: OutputResultsx.cpp:17
extrudefrombase_core
void extrudefrombase_core(FemModel *femmodel)
Definition: extrudefrombase_core.cpp:12
StressbalanceSIAAnalysis
Definition: StressbalanceSIAAnalysis.h:11
InputToDepthaverageOutEnum
@ InputToDepthaverageOutEnum
Definition: EnumDefinitions.h:204
AutodiffDependentObjectsEnum
@ AutodiffDependentObjectsEnum
Definition: EnumDefinitions.h:43
DatasetInput2
Definition: DatasetInput2.h:14
sealevelrise_core_masks
SealevelMasks * sealevelrise_core_masks(FemModel *femmodel)
Definition: sealevelchange_core.cpp:336
InversionNumControlParametersEnum
@ InversionNumControlParametersEnum
Definition: EnumDefinitions.h:223
ThermalRequestedOutputsEnum
@ ThermalRequestedOutputsEnum
Definition: EnumDefinitions.h:424
dummy_core
void dummy_core(FemModel *femmodel)
Definition: dummy_core.cpp:7
GroundinglineMigrationEnum
@ GroundinglineMigrationEnum
Definition: EnumDefinitions.h:161
LoveHrEnum
@ LoveHrEnum
Definition: EnumDefinitions.h:1148
ExtrudeFromBaseAnalysisEnum
@ ExtrudeFromBaseAnalysisEnum
Definition: EnumDefinitions.h:1058
MasstransportAnalysisEnum
@ MasstransportAnalysisEnum
Definition: EnumDefinitions.h:1163
ContactEnum
@ ContactEnum
Definition: EnumDefinitions.h:1014
BalancevelocityAnalysisEnum
@ BalancevelocityAnalysisEnum
Definition: EnumDefinitions.h:987
EsaNmotionEnum
@ EsaNmotionEnum
Definition: EnumDefinitions.h:561
_printf0_
#define _printf0_(StreamArgs)
Definition: Print.h:29
AdjointEnum
@ AdjointEnum
Definition: EnumDefinitions.h:464
DslModelidEnum
@ DslModelidEnum
Definition: EnumDefinitions.h:126
TransientRequestedOutputsEnum
@ TransientRequestedOutputsEnum
Definition: EnumDefinitions.h:457
HydrologySheetThicknessOldEnum
@ HydrologySheetThicknessOldEnum
Definition: EnumDefinitions.h:622
SmbAveragingEnum
@ SmbAveragingEnum
Definition: EnumDefinitions.h:349
PetscOptionsDetermineSolverType
void PetscOptionsDetermineSolverType(int *psolver_type)
Definition: PetscOptionsDetermineSolverType.cpp:20
transient_core
void transient_core(FemModel *femmodel)
Definition: transient_core.cpp:19
SurfaceSlopeXEnum
@ SurfaceSlopeXEnum
Definition: EnumDefinitions.h:829
StressbalanceAnalysisEnum
@ StressbalanceAnalysisEnum
Definition: EnumDefinitions.h:1285
BasalforcingsEnum
@ BasalforcingsEnum
Definition: EnumDefinitions.h:64
AppCtx
Definition: control_core.cpp:16
LoveKernelsImagEnum
@ LoveKernelsImagEnum
Definition: EnumDefinitions.h:1149
FemModel::InputMakeDiscontinuous
void InputMakeDiscontinuous(int enum_in)
Definition: FemModel.cpp:1611
StepEnum
@ StepEnum
Definition: EnumDefinitions.h:403
SurfaceOldEnum
@ SurfaceOldEnum
Definition: EnumDefinitions.h:824
VerboseConvergence
bool VerboseConvergence(void)
Definition: Verbosity.cpp:26
HydrologydcEnum
@ HydrologydcEnum
Definition: EnumDefinitions.h:1106
InversionNstepsEnum
@ InversionNstepsEnum
Definition: EnumDefinitions.h:222
LoveNfreqEnum
@ LoveNfreqEnum
Definition: EnumDefinitions.h:237
FormFunction
IssmDouble FormFunction(IssmDouble *X, void *usr)
Definition: control_core.cpp:116
thermal_core
void thermal_core(FemModel *femmodel)
Definition: thermal_core.cpp:13
DataSet::AddObject
int AddObject(Object *object)
Definition: DataSet.cpp:252
HydrologyRequestedOutputsEnum
@ HydrologyRequestedOutputsEnum
Definition: EnumDefinitions.h:173
FlowequationIsSSAEnum
@ FlowequationIsSSAEnum
Definition: EnumDefinitions.h:142
DrivingStressYEnum
@ DrivingStressYEnum
Definition: EnumDefinitions.h:539
LevelsetReinitFrequencyEnum
@ LevelsetReinitFrequencyEnum
Definition: EnumDefinitions.h:228
ExtrapolationAnalysis
Definition: ExtrapolationAnalysis.h:11
_printf_
#define _printf_(StreamArgs)
Definition: Print.h:22
SealevelUEsaEnum
@ SealevelUEsaEnum
Definition: EnumDefinitions.h:685
InversionControlParametersEnum
@ InversionControlParametersEnum
Definition: EnumDefinitions.h:209
AdjointCorePointerFromSolutionEnum
void AdjointCorePointerFromSolutionEnum(void(**padjointcore)(FemModel *), int solutiontype)
Definition: AdjointCorePointerFromSolutionEnum.cpp:18
FemModel::InitTransientInputx
void InitTransientInputx(int *transientinput_enum, int numoutputs)
Definition: FemModel.cpp:5202
AutodiffXpEnum
@ AutodiffXpEnum
Definition: EnumDefinitions.h:57
SealevelUNorthEsaEnum
@ SealevelUNorthEsaEnum
Definition: EnumDefinitions.h:687
HydraulicPotentialEnum
@ HydraulicPotentialEnum
Definition: EnumDefinitions.h:597
MaskOceanLevelsetEnum
@ MaskOceanLevelsetEnum
Definition: EnumDefinitions.h:640
Balancethickness2AnalysisEnum
@ Balancethickness2AnalysisEnum
Definition: EnumDefinitions.h:979
ExtrudeFromTopAnalysisEnum
@ ExtrudeFromTopAnalysisEnum
Definition: EnumDefinitions.h:1059
love_core
void love_core(FemModel *femmodel)
Definition: love_core.cpp:11
TransientSolutionEnum
@ TransientSolutionEnum
Definition: EnumDefinitions.h:1317
LevelsetfunctionSlopeXEnum
@ LevelsetfunctionSlopeXEnum
Definition: EnumDefinitions.h:635
adjointbalancethickness2_core
void adjointbalancethickness2_core(FemModel *femmodel)
Definition: adjointbalancethickness2_core.cpp:12
BalancethicknessSoftSolutionEnum
@ BalancethicknessSoftSolutionEnum
Definition: EnumDefinitions.h:984
FemModel::parameters
Parameters * parameters
Definition: FemModel.h:46
TimeEnum
@ TimeEnum
Definition: EnumDefinitions.h:427
ADD_VAL
@ ADD_VAL
Definition: toolkitsenums.h:14
sealevelrise_diagnostics
void sealevelrise_diagnostics(FemModel *femmodel, Vector< IssmDouble > *RSLg)
Definition: sealevelchange_core.cpp:658
IssmComm::GetComm
static ISSM_MPI_Comm GetComm(void)
Definition: IssmComm.cpp:30
TimesteppingTimeStepEnum
@ TimesteppingTimeStepEnum
Definition: EnumDefinitions.h:433
CorePointerFromSolutionEnum
void CorePointerFromSolutionEnum(void(**psolutioncore)(FemModel *), Parameters *parameters, int solutiontype)
Definition: CorePointerFromSolutionEnum.cpp:18
AutodiffFosForwardIndexEnum
@ AutodiffFosForwardIndexEnum
Definition: EnumDefinitions.h:45
DepthAverageAnalysisEnum
@ DepthAverageAnalysisEnum
Definition: EnumDefinitions.h:1036
AdaptiveTimesteppingEnum
@ AdaptiveTimesteppingEnum
Definition: EnumDefinitions.h:969
BalancevelocitySolutionEnum
@ BalancevelocitySolutionEnum
Definition: EnumDefinitions.h:988
LevelsetfunctionSlopeYEnum
@ LevelsetfunctionSlopeYEnum
Definition: EnumDefinitions.h:636
ThermalNumRequestedOutputsEnum
@ ThermalNumRequestedOutputsEnum
Definition: EnumDefinitions.h:419
DamageEvolutionAnalysisEnum
@ DamageEvolutionAnalysisEnum
Definition: EnumDefinitions.h:1026
BedEnum
@ BedEnum
Definition: EnumDefinitions.h:499
solutionsequence_schurcg
void solutionsequence_schurcg(FemModel *femmodel)
Definition: solutionsequence_schurcg.cpp:834
LoveKernelsRealEnum
@ LoveKernelsRealEnum
Definition: EnumDefinitions.h:1150
HydrologyShaktiAnalysis
Definition: HydrologyShaktiAnalysis.h:11
steadystateconvergence
bool steadystateconvergence(Vector< IssmDouble > *tg, Vector< IssmDouble > *tg_old, Vector< IssmDouble > *ug, Vector< IssmDouble > *ug_old, IssmDouble reltol)
Definition: steadystate_core.cpp:91
TransientInput2
Definition: TransientInput2.h:13
InversionIscontrolEnum
@ InversionIscontrolEnum
Definition: EnumDefinitions.h:218
FemModel::CostFunctionx
void CostFunctionx(IssmDouble *pJ, IssmDouble **pJlist, int *pn)
Definition: FemModel.cpp:1062
GROUNDINGLINECORE
#define GROUNDINGLINECORE
Definition: Profiler.h:25
SmbMassBalanceTransientEnum
@ SmbMassBalanceTransientEnum
Definition: EnumDefinitions.h:750
LoveLiEnum
@ LoveLiEnum
Definition: EnumDefinitions.h:1153
HydrologydcIsefficientlayerEnum
@ HydrologydcIsefficientlayerEnum
Definition: EnumDefinitions.h:185
GetVectorFromInputsx
void GetVectorFromInputsx(IssmDouble **pvector, int *pvector_size, FemModel *femmodel, int name)
Definition: GetVectorFromInputsx.cpp:81
JEnum
@ JEnum
Definition: EnumDefinitions.h:1134
FemModel::results
Results * results
Definition: FemModel.h:48
AutodiffDriverEnum
@ AutodiffDriverEnum
Definition: EnumDefinitions.h:44
DrivingStressXEnum
@ DrivingStressXEnum
Definition: EnumDefinitions.h:538
OptPars::maxiter
int * maxiter
Definition: OptPars.h:15
Vertices::NumberOfVertices
int NumberOfVertices(void)
Definition: Vertices.cpp:255
TransientIsslrEnum
@ TransientIsslrEnum
Definition: EnumDefinitions.h:452
HydrologydcEplThicknessSubstepEnum
@ HydrologydcEplThicknessSubstepEnum
Definition: EnumDefinitions.h:605
GeothermalFluxx
void GeothermalFluxx(FemModel *femmodel)
Definition: GeothermalFluxx.cpp:9
BalancethicknessSolutionEnum
@ BalancethicknessSolutionEnum
Definition: EnumDefinitions.h:985
UGiaRateEnum
@ UGiaRateEnum
Definition: EnumDefinitions.h:594
EsaYmotionEnum
@ EsaYmotionEnum
Definition: EnumDefinitions.h:568
OceanExchangeDatax
void OceanExchangeDatax(FemModel *femmodel, bool init_stage)
Definition: OceanExchangeDatax.cpp:11
balancevelocity_core
void balancevelocity_core(FemModel *femmodel)
Definition: balancevelocity_core.cpp:12
HydrologyHeadOldEnum
@ HydrologyHeadOldEnum
Definition: EnumDefinitions.h:616
Matlitho::issolid
IssmDouble * issolid
Definition: Matlitho.h:27
Matlitho::numlayers
int numlayers
Definition: Matlitho.h:18
LoveKernelsEnum
@ LoveKernelsEnum
Definition: EnumDefinitions.h:235
GiaModelEnum
@ GiaModelEnum
Definition: EnumDefinitions.h:155
SealevelRSLRateEnum
@ SealevelRSLRateEnum
Definition: EnumDefinitions.h:683
Vector::Duplicate
Vector< doubletype > * Duplicate(void)
Definition: Vector.h:230
FreeSurfaceBaseAnalysisEnum
@ FreeSurfaceBaseAnalysisEnum
Definition: EnumDefinitions.h:1071
FSSolverEnum
@ FSSolverEnum
Definition: EnumDefinitions.h:1061
DamageEvolutionSolutionEnum
@ DamageEvolutionSolutionEnum
Definition: EnumDefinitions.h:1027
SealevelNEsaEnum
@ SealevelNEsaEnum
Definition: EnumDefinitions.h:678
Material
Definition: Material.h:21
EnthalpyAnalysis
Definition: EnthalpyAnalysis.h:12
HydrologyShaktiAnalysisEnum
@ HydrologyShaktiAnalysisEnum
Definition: EnumDefinitions.h:1103
IceMaskNodeActivationEnum
@ IceMaskNodeActivationEnum
Definition: EnumDefinitions.h:627
BasalforcingsPicoEnum
@ BasalforcingsPicoEnum
Definition: EnumDefinitions.h:990
DslStericRateEnum
@ DslStericRateEnum
Definition: EnumDefinitions.h:543
EsaAnalysisEnum
@ EsaAnalysisEnum
Definition: EnumDefinitions.h:1053
MasstransportStabilizationEnum
@ MasstransportStabilizationEnum
Definition: EnumDefinitions.h:249
SettingsSbCouplingFrequencyEnum
@ SettingsSbCouplingFrequencyEnum
Definition: EnumDefinitions.h:339
EnthalpyEnum
@ EnthalpyEnum
Definition: EnumDefinitions.h:551
BedSlopeXEnum
@ BedSlopeXEnum
Definition: EnumDefinitions.h:500
esa_core
void esa_core(FemModel *femmodel)
Definition: esa_core.cpp:12
AmrRestartEnum
@ AmrRestartEnum
Definition: EnumDefinitions.h:29
MasstransportRequestedOutputsEnum
@ MasstransportRequestedOutputsEnum
Definition: EnumDefinitions.h:248
HydrologydcEplThicknessTransientEnum
@ HydrologydcEplThicknessTransientEnum
Definition: EnumDefinitions.h:606
ConstantsYtsEnum
@ ConstantsYtsEnum
Definition: EnumDefinitions.h:104
HydrologyPismAnalysisEnum
@ HydrologyPismAnalysisEnum
Definition: EnumDefinitions.h:1102
SteadystateNumRequestedOutputsEnum
@ SteadystateNumRequestedOutputsEnum
Definition: EnumDefinitions.h:400
FemModel::OutputControlsx
void OutputControlsx(Results **presults)
Definition: FemModel.cpp:2171
IcecapToEarthCommEnum
@ IcecapToEarthCommEnum
Definition: EnumDefinitions.h:200
FemModel::TimeAdaptx
void TimeAdaptx(IssmDouble *pdt)
Definition: FemModel.cpp:2895
SealevelriseSolutionEnum
@ SealevelriseSolutionEnum
Definition: EnumDefinitions.h:1267
AutodiffFovForwardIndicesEnum
@ AutodiffFovForwardIndicesEnum
Definition: EnumDefinitions.h:47
SolidearthSettingsHorizEnum
@ SolidearthSettingsHorizEnum
Definition: EnumDefinitions.h:323
VerboseModule
bool VerboseModule(void)
Definition: Verbosity.cpp:23
TransientIsgroundinglineEnum
@ TransientIsgroundinglineEnum
Definition: EnumDefinitions.h:447
DamageStabilizationEnum
@ DamageStabilizationEnum
Definition: EnumDefinitions.h:119
UGiaEnum
@ UGiaEnum
Definition: EnumDefinitions.h:593
EsaUmotionEnum
@ EsaUmotionEnum
Definition: EnumDefinitions.h:566
SealevelriseRunCountEnum
@ SealevelriseRunCountEnum
Definition: EnumDefinitions.h:331
SteadystateRequestedOutputsEnum
@ SteadystateRequestedOutputsEnum
Definition: EnumDefinitions.h:402
DamageEvolutionNumRequestedOutputsEnum
@ DamageEvolutionNumRequestedOutputsEnum
Definition: EnumDefinitions.h:113
EarthMassTransport
void EarthMassTransport(FemModel *femmodel)
Definition: sealevelchange_core.cpp:983
grd_core
void grd_core(FemModel *femmodel)
Definition: sealevelchange_core.cpp:84
VyEnum
@ VyEnum
Definition: EnumDefinitions.h:850
solutionsequence_linear
void solutionsequence_linear(FemModel *femmodel)
Definition: solutionsequence_linear.cpp:10
HydrologydcEplThicknessOldEnum
@ HydrologydcEplThicknessOldEnum
Definition: EnumDefinitions.h:604
steadystate_core
void steadystate_core(FemModel *femmodel)
Definition: steadystate_core.cpp:20
FlowequationFeFSEnum
@ FlowequationFeFSEnum
Definition: EnumDefinitions.h:137
InversionMaxiterPerStepEnum
@ InversionMaxiterPerStepEnum
Definition: EnumDefinitions.h:220
EplHeadEnum
@ EplHeadEnum
Definition: EnumDefinitions.h:553
SmoothAnalysisEnum
@ SmoothAnalysisEnum
Definition: EnumDefinitions.h:1276
solutionsequence_nonlinear
void solutionsequence_nonlinear(FemModel *femmodel, bool conserve_loads)
Definition: solutionsequence_nonlinear.cpp:11
Input2::ChangeEnum
void ChangeEnum(int newenumtype)
Definition: Input2.h:26
NORM_TWO
@ NORM_TWO
Definition: toolkitsenums.h:15
FemModel::RequestedDependentsx
void RequestedDependentsx(void)
Definition: FemModel.cpp:2220
solutionsequence_shakti_nonlinear
void solutionsequence_shakti_nonlinear(FemModel *femmodel)
Definition: solutionsequence_shakti_nonlinear.cpp:11
ResetBoundaryConditions
void ResetBoundaryConditions(FemModel *femmodel, int analysis_type)
Definition: ResetBoundaryConditions.cpp:9
hydrology_core
void hydrology_core(FemModel *femmodel)
Definition: hydrology_core.cpp:12
HydrologyDCInefficientAnalysisEnum
@ HydrologyDCInefficientAnalysisEnum
Definition: EnumDefinitions.h:1099
AdjointpEnum
@ AdjointpEnum
Definition: EnumDefinitions.h:465
AppCtx::nsize
int nsize
Definition: control_core.cpp:18
SmbStepsPerStepEnum
@ SmbStepsPerStepEnum
Definition: EnumDefinitions.h:389
LoveKiEnum
@ LoveKiEnum
Definition: EnumDefinitions.h:1151
BrentSearch
void BrentSearch(IssmDouble **pJ, OptPars optpars, IssmDouble *X0, IssmDouble(*f)(IssmDouble *, void *), IssmDouble(*g)(IssmDouble **, IssmDouble *, void *), void *usr)
Definition: BrentSearch.cpp:23
Vector::Norm
doubletype Norm(NormMode norm_type)
Definition: Vector.h:342
DependentObject::AddValue
void AddValue(IssmDouble in_value)
Definition: DependentObject.cpp:106
Vector::GetSize
void GetSize(int *pM)
Definition: Vector.h:185
control_core
void control_core(FemModel *femmodel)
Definition: control_core.cpp:22
SmbSolutionEnum
@ SmbSolutionEnum
Definition: EnumDefinitions.h:1275
SealevelriseTransitionsEnum
@ SealevelriseTransitionsEnum
Definition: EnumDefinitions.h:332
sealevelrise_core_geometry
void sealevelrise_core_geometry(FemModel *femmodel)
Definition: sealevelchange_core.cpp:351
HydrologyGlaDSAnalysisEnum
@ HydrologyGlaDSAnalysisEnum
Definition: EnumDefinitions.h:1100
Element
Definition: Element.h:41
AutodiffJacobianEnum
@ AutodiffJacobianEnum
Definition: EnumDefinitions.h:978
Vector::Scale
void Scale(doubletype scale_factor)
Definition: Vector.h:355
OptPars::xmax
IssmDouble xmax
Definition: OptPars.h:13
solutionsequence_la
void solutionsequence_la(FemModel *femmodel)
Definition: solutionsequence_la.cpp:10
TriaInput2::copy
Input2 * copy()
Definition: TriaInput2.cpp:72
Domain2DhorizontalEnum
@ Domain2DhorizontalEnum
Definition: EnumDefinitions.h:534
DslModelEnum
@ DslModelEnum
Definition: EnumDefinitions.h:125
SolidearthSettingsReltolEnum
@ SolidearthSettingsReltolEnum
Definition: EnumDefinitions.h:327
SealevelriseCumDeltathicknessEnum
@ SealevelriseCumDeltathicknessEnum
Definition: EnumDefinitions.h:688
ISSM_MPI_DOUBLE
#define ISSM_MPI_DOUBLE
Definition: issmmpi.h:125
ThermalSolutionEnum
@ ThermalSolutionEnum
Definition: EnumDefinitions.h:1303
DslComputeFingerprintsEnum
@ DslComputeFingerprintsEnum
Definition: EnumDefinitions.h:128
MeltingAnalysisEnum
@ MeltingAnalysisEnum
Definition: EnumDefinitions.h:1182
InputExtrudex
void InputExtrudex(FemModel *femmodel, int input_enum, int start)
Definition: InputExtrudex.cpp:10
HydrologyGlaDSAnalysis::UpdateEffectivePressure
void UpdateEffectivePressure(FemModel *femmodel)
Definition: HydrologyGlaDSAnalysis.cpp:547
LoveShNmaxEnum
@ LoveShNmaxEnum
Definition: EnumDefinitions.h:239
sealevelrise_core_geometry
void sealevelrise_core_geometry(FemModel *femmodel)
Definition: sealevelchange_core.cpp:351
LoveAllowLayerDeletionEnum
@ LoveAllowLayerDeletionEnum
Definition: EnumDefinitions.h:231
StressbalanceVerticalAnalysis
Definition: StressbalanceVerticalAnalysis.h:11
Parameters::SetParam
void SetParam(bool boolean, int enum_type)
Definition: Parameters.cpp:441
LATaylorHoodEnum
@ LATaylorHoodEnum
Definition: EnumDefinitions.h:1139
SedimentHeadTransientEnum
@ SedimentHeadTransientEnum
Definition: EnumDefinitions.h:701
ISSM_MPI_INT
#define ISSM_MPI_INT
Definition: issmmpi.h:127
GenericParam::GetParameterValue
P & GetParameterValue()
Definition: GenericParam.h:62
SurfaceAreax
void SurfaceAreax(IssmDouble *pS, FemModel *femmodel)
Definition: SurfaceAreax.cpp:11
ModelIdEnum
@ ModelIdEnum
Definition: EnumDefinitions.h:274
SolidearthSettingsMaxiterEnum
@ SolidearthSettingsMaxiterEnum
Definition: EnumDefinitions.h:324
Vector::AXPY
void AXPY(Vector *X, doubletype a)
Definition: Vector.h:256
SMBCORE
#define SMBCORE
Definition: Profiler.h:24
VerboseAutodiff
bool VerboseAutodiff(void)
Definition: Verbosity.cpp:29
DependentObject
Definition: DependentObject.h:15
AdjointBalancethicknessAnalysisEnum
@ AdjointBalancethicknessAnalysisEnum
Definition: EnumDefinitions.h:971
GiaSolutionEnum
@ GiaSolutionEnum
Definition: EnumDefinitions.h:1085
sealevelchange_core
void sealevelchange_core(FemModel *femmodel)
Definition: sealevelchange_core.cpp:17
FemModel::HydrologyIDSupdateDomainx
void HydrologyIDSupdateDomainx(IssmDouble *pIDScount)
Definition: FemModel.cpp:5072
LoveForcingTypeEnum
@ LoveForcingTypeEnum
Definition: EnumDefinitions.h:232
L2ProjectionBaseAnalysisEnum
@ L2ProjectionBaseAnalysisEnum
Definition: EnumDefinitions.h:1136
VertexSIdEnum
@ VertexSIdEnum
Definition: EnumDefinitions.h:1325
GiaAnalysisEnum
@ GiaAnalysisEnum
Definition: EnumDefinitions.h:1084
BasalforcingsGroundediceMeltingRateEnum
@ BasalforcingsGroundediceMeltingRateEnum
Definition: EnumDefinitions.h:478
ExtrapolationVariableEnum
@ ExtrapolationVariableEnum
Definition: EnumDefinitions.h:135
InputToExtrudeEnum
@ InputToExtrudeEnum
Definition: EnumDefinitions.h:205
EplHeadTransientEnum
@ EplHeadTransientEnum
Definition: EnumDefinitions.h:558
TimesteppingTypeEnum
@ TimesteppingTypeEnum
Definition: EnumDefinitions.h:436
solutionsequence_glads_nonlinear
void solutionsequence_glads_nonlinear(FemModel *femmodel)
Definition: solutionsequence_glads_nonlinear.cpp:11
Vector::SetValues
void SetValues(int ssize, int *list, doubletype *values, InsMode mode)
Definition: Vector.h:153
DomainTypeEnum
@ DomainTypeEnum
Definition: EnumDefinitions.h:124
extrudefromtop_core
void extrudefromtop_core(FemModel *femmodel)
Definition: extrudefromtop_core.cpp:12
VxAverageEnum
@ VxAverageEnum
Definition: EnumDefinitions.h:845
DAMAGECORE
#define DAMAGECORE
Definition: Profiler.h:21
FemModel::nodes
Nodes * nodes
Definition: FemModel.h:56
SealevelInertiaTensorXZEnum
@ SealevelInertiaTensorXZEnum
Definition: EnumDefinitions.h:1261
InputUpdateFromVectorx
void InputUpdateFromVectorx(FemModel *femmodel, Vector< IssmDouble > *vector, int name, int type)
Definition: InputUpdateFromVectorx.cpp:9
AdjointHorizAnalysisEnum
@ AdjointHorizAnalysisEnum
Definition: EnumDefinitions.h:972
xDelete
void xDelete(T **&aT_pp)
Definition: MemOps.h:97
EarthIdEnum
@ EarthIdEnum
Definition: EnumDefinitions.h:129
AdjointyEnum
@ AdjointyEnum
Definition: EnumDefinitions.h:467
KillIcebergsx
void KillIcebergsx(FemModel *femmodel)
Definition: KillIcebergsx.cpp:11
SmbMassBalanceSubstepEnum
@ SmbMassBalanceSubstepEnum
Definition: EnumDefinitions.h:749
FemModel::analysis_counter
int analysis_counter
Definition: FemModel.h:37
HydrologySheetThicknessEnum
@ HydrologySheetThicknessEnum
Definition: EnumDefinitions.h:621
ResetFSBasalBoundaryConditionx
void ResetFSBasalBoundaryConditionx(Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters)
Definition: ResetFSBasalBoundaryConditionx.cpp:9
VzEnum
@ VzEnum
Definition: EnumDefinitions.h:853
FemModel::materials
Materials * materials
Definition: FemModel.h:45
EnumToStringx
const char * EnumToStringx(int enum_in)
Definition: EnumToStringx.cpp:15
TriaInput2
Definition: TriaInput2.h:8
TransientIssmbEnum
@ TransientIssmbEnum
Definition: EnumDefinitions.h:453
HydrologyPismAnalysis
Definition: HydrologyPismAnalysis.h:11
SealevelRSLEustaticEnum
@ SealevelRSLEustaticEnum
Definition: EnumDefinitions.h:681
smb_core
void smb_core(FemModel *femmodel)
Definition: smb_core.cpp:12
ThermalIsenthalpyEnum
@ ThermalIsenthalpyEnum
Definition: EnumDefinitions.h:417
InputToSmoothEnum
@ InputToSmoothEnum
Definition: EnumDefinitions.h:207
NGiaRateEnum
@ NGiaRateEnum
Definition: EnumDefinitions.h:592
IssmComm::GetSize
static int GetSize(void)
Definition: IssmComm.cpp:46
slrconvergence
void slrconvergence(bool *pconverged, Vector< IssmDouble > *RSLg, Vector< IssmDouble > *RSLg_old, IssmDouble eps_rel, IssmDouble eps_abs)
Definition: sealevelchange_core.cpp:1031
FloatingiceMeltingRatex
void FloatingiceMeltingRatex(FemModel *femmodel)
Definition: FloatingiceMeltingRatex.cpp:10
SMBgradientscomponentsEnum
@ SMBgradientscomponentsEnum
Definition: EnumDefinitions.h:1248
LoveHiEnum
@ LoveHiEnum
Definition: EnumDefinitions.h:1147
SmoothThicknessMultiplierEnum
@ SmoothThicknessMultiplierEnum
Definition: EnumDefinitions.h:397
GsetEnum
@ GsetEnum
Definition: EnumDefinitions.h:1093
FemModel::analysis_type_list
int * analysis_type_list
Definition: FemModel.h:38
controladm1qn3_core
void controladm1qn3_core(FemModel *femmodel)
Definition: controladm1qn3_core.cpp:657
BasalforcingsFloatingiceMeltingRateEnum
@ BasalforcingsFloatingiceMeltingRateEnum
Definition: EnumDefinitions.h:476
Domain3DsurfaceEnum
@ Domain3DsurfaceEnum
Definition: EnumDefinitions.h:1039
TransientIsgiaEnum
@ TransientIsgiaEnum
Definition: EnumDefinitions.h:446
alpha
IssmDouble alpha(IssmDouble x, IssmDouble y, IssmDouble z, int testid)
Definition: fsanalyticals.cpp:221
StressbalanceSolutionEnum
@ StressbalanceSolutionEnum
Definition: EnumDefinitions.h:1288
solutionsequence_thermal_nonlinear
void solutionsequence_thermal_nonlinear(FemModel *femmodel)
Definition: solutionsequence_thermal_nonlinear.cpp:11
MASSTRANSPORTCORE
#define MASSTRANSPORTCORE
Definition: Profiler.h:23
TransientIsoceancouplingEnum
@ TransientIsoceancouplingEnum
Definition: EnumDefinitions.h:451
sealevelrise_core_eustatic
Vector< IssmDouble > * sealevelrise_core_eustatic(FemModel *femmodel, SealevelMasks *masks, IssmDouble *poceanarea)
Definition: sealevelchange_core.cpp:404
gia_core
void gia_core(FemModel *femmodel)
Definition: gia_core.cpp:13
MasstransportSolutionEnum
@ MasstransportSolutionEnum
Definition: EnumDefinitions.h:1164
AdjointxEnum
@ AdjointxEnum
Definition: EnumDefinitions.h:466
AutodiffFosReverseIndexEnum
@ AutodiffFosReverseIndexEnum
Definition: EnumDefinitions.h:46
MOVINGFRONTCORE
#define MOVINGFRONTCORE
Definition: Profiler.h:22
FemModel::AverageTransientInputx
void AverageTransientInputx(int *transientinput_enum, int *averagedinput_enum, IssmDouble init_time, IssmDouble end_time, int numoutputs, int averaging_method)
Definition: FemModel.cpp:5246
MeshElementtypeEnum
@ MeshElementtypeEnum
Definition: EnumDefinitions.h:271
ISSM_MPI_Send
int ISSM_MPI_Send(void *buf, int count, ISSM_MPI_Datatype datatype, int dest, int tag, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:484
FemModel::inputs2
Inputs2 * inputs2
Definition: FemModel.h:47
GIACORE
#define GIACORE
Definition: Profiler.h:26
Matlitho::lame_lambda
IssmDouble * lame_lambda
Definition: Matlitho.h:21
Gradjx
void Gradjx(Vector< IssmDouble > **pgradient, IssmDouble **pnorm_list, Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters)
Definition: Gradjx.cpp:9
Damagex
void Damagex(FemModel *femmodel)
Definition: Damagex.cpp:10
QmuIsdakotaEnum
@ QmuIsdakotaEnum
Definition: EnumDefinitions.h:288
sealevelrise_core_viscous
void sealevelrise_core_viscous(Vector< IssmDouble > **pU_gia, Vector< IssmDouble > **pN_gia, FemModel *femmodel, Vector< IssmDouble > *RSLg)
Definition: sealevelchange_core.cpp:628
FemModel::constraints
Constraints * constraints
Definition: FemModel.h:52
AutodiffIsautodiffEnum
@ AutodiffIsautodiffEnum
Definition: EnumDefinitions.h:50
InputDuplicatex
void InputDuplicatex(FemModel *femmodel, int original_enum, int new_enum)
Definition: InputDuplicatex.cpp:10
GetMaskOfIceVerticesLSMx
void GetMaskOfIceVerticesLSMx(FemModel *femmodel)
Definition: SetActiveNodesLSMx.cpp:82
FixedTimesteppingEnum
@ FixedTimesteppingEnum
Definition: EnumDefinitions.h:1066
TransferForcing
void TransferForcing(FemModel *femmodel, int forcingenum)
Definition: sealevelchange_core.cpp:752
SolutionTypeEnum
@ SolutionTypeEnum
Definition: EnumDefinitions.h:398
FemModel::UpdateConstraintsExtrudeFromBasex
void UpdateConstraintsExtrudeFromBasex()
Definition: FemModel.cpp:3009
ISSM_MPI_Status
int ISSM_MPI_Status
Definition: issmmpi.h:121
WaterColumnOldEnum
@ WaterColumnOldEnum
Definition: EnumDefinitions.h:858
SurfaceEnum
@ SurfaceEnum
Definition: EnumDefinitions.h:823
TransientIshydrologyEnum
@ TransientIshydrologyEnum
Definition: EnumDefinitions.h:448
HydrologyAveragingEnum
@ HydrologyAveragingEnum
Definition: EnumDefinitions.h:162
THERMALCORE
#define THERMALCORE
Definition: Profiler.h:18
TransientAmrFrequencyEnum
@ TransientAmrFrequencyEnum
Definition: EnumDefinitions.h:442
SealevelriseAnalysisEnum
@ SealevelriseAnalysisEnum
Definition: EnumDefinitions.h:1266
Results::AddResult
int AddResult(ExternalResult *result)
Definition: Results.cpp:33
LoveKrEnum
@ LoveKrEnum
Definition: EnumDefinitions.h:1152
dynstr_core
void dynstr_core(FemModel *femmodel)
Definition: sealevelchange_core.cpp:259
InputUpdateFromConstantx
void InputUpdateFromConstantx(FemModel *femmodel, bool constant, int name)
Definition: InputUpdateFromConstantx.cpp:10
TransientIscouplerEnum
@ TransientIscouplerEnum
Definition: EnumDefinitions.h:443
FemModel::loads
Loads * loads
Definition: FemModel.h:54
Profiler::Stop
void Stop(int tagenum, bool dontmpisync=true)
Definition: Profiler.cpp:179
SealevelriseGeometryDoneEnum
@ SealevelriseGeometryDoneEnum
Definition: EnumDefinitions.h:309
SettingsRecordingFrequencyEnum
@ SettingsRecordingFrequencyEnum
Definition: EnumDefinitions.h:337
BaseSlopeXEnum
@ BaseSlopeXEnum
Definition: EnumDefinitions.h:497
HydrologyNumRequestedOutputsEnum
@ HydrologyNumRequestedOutputsEnum
Definition: EnumDefinitions.h:170
TransientIsthermalEnum
@ TransientIsthermalEnum
Definition: EnumDefinitions.h:455
Domain3DEnum
@ Domain3DEnum
Definition: EnumDefinitions.h:536
masstransport_core
void masstransport_core(FemModel *femmodel)
Definition: masstransport_core.cpp:12
controlm1qn3_core
void controlm1qn3_core(FemModel *femmodel)
Definition: controlm1qn3_core.cpp:304
FemModel::StackTransientInputx
void StackTransientInputx(int *input_enum, int *transientinput_enum, IssmDouble hydrotime, int numoutputs)
Definition: FemModel.cpp:5213
OptPars
Definition: OptPars.h:10
DslDynamicRateEnum
@ DslDynamicRateEnum
Definition: EnumDefinitions.h:544
FemModel::elements
Elements * elements
Definition: FemModel.h:44
LevelsetAnalysis
Definition: LevelsetAnalysis.h:11
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
DependentObject::Responsex
void Responsex(IssmDouble *poutput_value, FemModel *femmodel)
Definition: DependentObject.cpp:93
HydrologyShreveAnalysisEnum
@ HydrologyShreveAnalysisEnum
Definition: EnumDefinitions.h:1104
GiaMmeUgiaEnum
@ GiaMmeUgiaEnum
Definition: EnumDefinitions.h:546
LoveR0Enum
@ LoveR0Enum
Definition: EnumDefinitions.h:238
TransientIsstressbalanceEnum
@ TransientIsstressbalanceEnum
Definition: EnumDefinitions.h:454
GiaMmeNgiaEnum
@ GiaMmeNgiaEnum
Definition: EnumDefinitions.h:545
SLRCORE
#define SLRCORE
Definition: Profiler.h:28
Object::ObjectEnum
virtual int ObjectEnum()=0
TimesteppingStartTimeEnum
@ TimesteppingStartTimeEnum
Definition: EnumDefinitions.h:432
surfaceslope_core
void surfaceslope_core(FemModel *femmodel)
Definition: surfaceslope_core.cpp:12
dakota_core
void dakota_core(FemModel *femmodel)
Definition: dakota_core.cpp:313
Vector::Assemble
void Assemble(void)
Definition: Vector.h:142
BaseSlopeYEnum
@ BaseSlopeYEnum
Definition: EnumDefinitions.h:498
StringToEnumx
int StringToEnumx(const char *string_in, bool notfounderror=true)
Definition: StringToEnumx.cpp:14
Input2
Definition: Input2.h:18
Matlitho::density
IssmDouble * density
Definition: Matlitho.h:25
SealevelRSLEnum
@ SealevelRSLEnum
Definition: EnumDefinitions.h:680
Matlitho::isburgers
IssmDouble * isburgers
Definition: Matlitho.h:26
FemModel
Definition: FemModel.h:31
controltao_core
void controltao_core(FemModel *femmodel)
Definition: controltao_core.cpp:213
EsaEmotionEnum
@ EsaEmotionEnum
Definition: EnumDefinitions.h:560
bmb_core
void bmb_core(FemModel *femmodel)
Definition: bmb_core.cpp:12
LevelsetKillIcebergsEnum
@ LevelsetKillIcebergsEnum
Definition: EnumDefinitions.h:227
ThermalAnalysisEnum
@ ThermalAnalysisEnum
Definition: EnumDefinitions.h:1302
LoveSolutionEnum
@ LoveSolutionEnum
Definition: EnumDefinitions.h:1155
TemperatureEnum
@ TemperatureEnum
Definition: EnumDefinitions.h:831
InversionStepThresholdEnum
@ InversionStepThresholdEnum
Definition: EnumDefinitions.h:225
MatlithoEnum
@ MatlithoEnum
Definition: EnumDefinitions.h:1170
SealevelEnum
@ SealevelEnum
Definition: EnumDefinitions.h:675
SedimentHeadSubstepEnum
@ SedimentHeadSubstepEnum
Definition: EnumDefinitions.h:700
HydrologydcEplThicknessEnum
@ HydrologydcEplThicknessEnum
Definition: EnumDefinitions.h:603
Inputs2::GetTransientInput
TransientInput2 * GetTransientInput(int enum_type)
Definition: Inputs2.cpp:406
StressbalanceSIAAnalysisEnum
@ StressbalanceSIAAnalysisEnum
Definition: EnumDefinitions.h:1287
SmbRunoffEnum
@ SmbRunoffEnum
Definition: EnumDefinitions.h:772
SmbAnalysis
Definition: SmbAnalysis.h:11
EnthalpyAnalysis::Core
void Core(FemModel *femmodel)
Definition: EnthalpyAnalysis.cpp:531
HydrologyshaktiEnum
@ HydrologyshaktiEnum
Definition: EnumDefinitions.h:1108
GenericExternalResult
Definition: GenericExternalResult.h:21
OptPars::nsteps
int nsteps
Definition: OptPars.h:16
IssmComm::GetRank
static int GetRank(void)
Definition: IssmComm.cpp:34
adjointbalancethickness_core
void adjointbalancethickness_core(FemModel *femmodel)
Definition: adjointbalancethickness_core.cpp:12
LACrouzeixRaviartEnum
@ LACrouzeixRaviartEnum
Definition: EnumDefinitions.h:1138
LoveFrequenciesEnum
@ LoveFrequenciesEnum
Definition: EnumDefinitions.h:233
FemModel::ResetLevelset
void ResetLevelset()
Definition: FemModel.cpp:2545
StressbalanceRequestedOutputsEnum
@ StressbalanceRequestedOutputsEnum
Definition: EnumDefinitions.h:411
FemModel::CheckPoint
void CheckPoint(void)
Definition: FemModel.cpp:238
TransientIsesaEnum
@ TransientIsesaEnum
Definition: EnumDefinitions.h:445
MasstransportIsfreesurfaceEnum
@ MasstransportIsfreesurfaceEnum
Definition: EnumDefinitions.h:244
sealevelrise_core_elastic
void sealevelrise_core_elastic(Vector< IssmDouble > **pU_esa, Vector< IssmDouble > **pU_north_esa, Vector< IssmDouble > **pU_east_esa, FemModel *femmodel, Vector< IssmDouble > *RSLg, SealevelMasks *masks)
Definition: sealevelchange_core.cpp:572
SealevelInertiaTensorZZEnum
@ SealevelInertiaTensorZZEnum
Definition: EnumDefinitions.h:1263
HydrologyModelEnum
@ HydrologyModelEnum
Definition: EnumDefinitions.h:169
InversionControlScalingFactorsEnum
@ InversionControlScalingFactorsEnum
Definition: EnumDefinitions.h:210
InputToDepthaverageInEnum
@ InputToDepthaverageInEnum
Definition: EnumDefinitions.h:203
Vector::AYPX
void AYPX(Vector *X, doubletype a)
Definition: Vector.h:267
GLheightadvectionAnalysisEnum
@ GLheightadvectionAnalysisEnum
Definition: EnumDefinitions.h:1077
OptPars::xmin
IssmDouble xmin
Definition: OptPars.h:12
ThicknessOldEnum
@ ThicknessOldEnum
Definition: EnumDefinitions.h:841
StressbalanceNumRequestedOutputsEnum
@ StressbalanceNumRequestedOutputsEnum
Definition: EnumDefinitions.h:408
HydrologyPismAnalysis::UpdateWaterColumn
void UpdateWaterColumn(FemModel *femmodel)
Definition: HydrologyPismAnalysis.cpp:92
SealevelriseRequestedOutputsEnum
@ SealevelriseRequestedOutputsEnum
Definition: EnumDefinitions.h:328
_error_
#define _error_(StreamArgs)
Definition: exceptions.h:49
EsaXmotionEnum
@ EsaXmotionEnum
Definition: EnumDefinitions.h:567
FlowequationIsFSEnum
@ FlowequationIsFSEnum
Definition: EnumDefinitions.h:138
Profiler::Start
void Start(int tagenum, bool dontmpisync=true)
Definition: Profiler.cpp:139
SurfaceSlopeSolutionEnum
@ SurfaceSlopeSolutionEnum
Definition: EnumDefinitions.h:1298
GetSolutionFromInputsx
void GetSolutionFromInputsx(Vector< IssmDouble > **psolution, FemModel *femmodel)
Definition: GetSolutionFromInputsx.cpp:9
SolidearthSettingsRotationEnum
@ SolidearthSettingsRotationEnum
Definition: EnumDefinitions.h:330
VerboseSolution
bool VerboseSolution(void)
Definition: Verbosity.cpp:24
DataSet::GetObjectByOffset
Object * GetObjectByOffset(int offset)
Definition: DataSet.cpp:334
VxEnum
@ VxEnum
Definition: EnumDefinitions.h:846
HydrologypismEnum
@ HydrologypismEnum
Definition: EnumDefinitions.h:1107
SurfaceloadIceThicknessChangeEnum
@ SurfaceloadIceThicknessChangeEnum
Definition: EnumDefinitions.h:691
WatercolumnEnum
@ WatercolumnEnum
Definition: EnumDefinitions.h:859
TransientIsmovingfrontEnum
@ TransientIsmovingfrontEnum
Definition: EnumDefinitions.h:450
FemModel::Restart
void Restart(void)
Definition: FemModel.cpp:549
EsaSolutionEnum
@ EsaSolutionEnum
Definition: EnumDefinitions.h:1054
InversionNumCostFunctionsEnum
@ InversionNumCostFunctionsEnum
Definition: EnumDefinitions.h:224
SettingsOutputFrequencyEnum
@ SettingsOutputFrequencyEnum
Definition: EnumDefinitions.h:336
PentaEnum
@ PentaEnum
Definition: EnumDefinitions.h:1231
ISSM_MPI_Comm
int ISSM_MPI_Comm
Definition: issmmpi.h:118
AdolcParamEnum
@ AdolcParamEnum
Definition: EnumDefinitions.h:12
SolidearthSettingsRunFrequencyEnum
@ SolidearthSettingsRunFrequencyEnum
Definition: EnumDefinitions.h:321
SmbRunoffTransientEnum
@ SmbRunoffTransientEnum
Definition: EnumDefinitions.h:774
AppCtx::femmodel
FemModel * femmodel
Definition: control_core.cpp:17
LoveLrEnum
@ LoveLrEnum
Definition: EnumDefinitions.h:1154
DamageEvolutionRequestedOutputsEnum
@ DamageEvolutionRequestedOutputsEnum
Definition: EnumDefinitions.h:114
GroundinglineMigrationx
void GroundinglineMigrationx(Elements *elements, Nodes *nodes, Vertices *vertices, Loads *loads, Materials *materials, Parameters *parameters)
Definition: GroundinglineMigrationx.cpp:10
AutodiffNumDependentsEnum
@ AutodiffNumDependentsEnum
Definition: EnumDefinitions.h:52
FemModel::RequestedOutputsx
void RequestedOutputsx(Results **presults, char **requested_outputs, int numoutputs, bool save_results=true)
Definition: FemModel.cpp:2267
FemModel::SetCurrentConfiguration
void SetCurrentConfiguration(int configuration_type)
Definition: FemModel.cpp:634
SmbNumRequestedOutputsEnum
@ SmbNumRequestedOutputsEnum
Definition: EnumDefinitions.h:379
VyAverageEnum
@ VyAverageEnum
Definition: EnumDefinitions.h:849
Matlitho
Definition: Matlitho.h:14
STRESSBALANCECORE
#define STRESSBALANCECORE
Definition: Profiler.h:20
damage_core
void damage_core(FemModel *femmodel)
Definition: damage_core.cpp:12
HydrologyGlaDSAnalysis
Definition: HydrologyGlaDSAnalysis.h:11
Matlitho::burgers_viscosity
IssmDouble * burgers_viscosity
Definition: Matlitho.h:23
HydrologyStepsPerStepEnum
@ HydrologyStepsPerStepEnum
Definition: EnumDefinitions.h:175
Parameters::FindParam
void FindParam(bool *pinteger, int enum_type)
Definition: Parameters.cpp:262
levelsetfunctionslope_core
void levelsetfunctionslope_core(FemModel *femmodel)
Definition: levelsetfunctionslope_core.cpp:12
ThicknessEnum
@ ThicknessEnum
Definition: EnumDefinitions.h:840
balancethickness2_core
void balancethickness2_core(FemModel *femmodel)
Definition: balancethickness2_core.cpp:12
InversionTypeEnum
@ InversionTypeEnum
Definition: EnumDefinitions.h:226
OptPars::cm_jump
IssmDouble * cm_jump
Definition: OptPars.h:14
SetControlInputsFromVectorx
void SetControlInputsFromVectorx(FemModel *femmodel, IssmDouble *vector)
Definition: SetControlInputsFromVectorx.cpp:9
InputToL2ProjectEnum
@ InputToL2ProjectEnum
Definition: EnumDefinitions.h:206
SteadystateMaxiterEnum
@ SteadystateMaxiterEnum
Definition: EnumDefinitions.h:399
Matlitho::radius
IssmDouble * radius
Definition: Matlitho.h:19
bedslope_core
void bedslope_core(FemModel *femmodel)
Definition: bedslope_core.cpp:12
HydrologyGlaDSEnum
@ HydrologyGlaDSEnum
Definition: EnumDefinitions.h:1101
VelEnum
@ VelEnum
Definition: EnumDefinitions.h:844
Matlitho::lame_mu
IssmDouble * lame_mu
Definition: Matlitho.h:22
solutionsequence_hydro_nonlinear
void solutionsequence_hydro_nonlinear(FemModel *femmodel)
Definition: solutionsequence_hydro_nonlinear.cpp:12
Inputs2::GetDatasetInput2
DatasetInput2 * GetDatasetInput2(int enum_type)
Definition: Inputs2.cpp:438
EplHeadOldEnum
@ EplHeadOldEnum
Definition: EnumDefinitions.h:554
FormFunctionGradient
IssmDouble FormFunctionGradient(IssmDouble **pG, IssmDouble *X, void *usr)
Definition: control_core.cpp:197
SteadystateSolutionEnum
@ SteadystateSolutionEnum
Definition: EnumDefinitions.h:1283
Calvingx
void Calvingx(FemModel *femmodel)
Definition: Calvingx.cpp:9
EsaRequestedOutputsEnum
@ EsaRequestedOutputsEnum
Definition: EnumDefinitions.h:133
VertexCoordinatesx
void VertexCoordinatesx(IssmDouble **px, IssmDouble **py, IssmDouble **pz, Vertices *vertices, bool spherical)
Definition: VertexCoordinatesx.cpp:11
SMBpddSicopolisEnum
@ SMBpddSicopolisEnum
Definition: EnumDefinitions.h:1253
InputDepthAverageAtBasex
void InputDepthAverageAtBasex(FemModel *femmodel, int original_enum, int new_enum)
Definition: InputDepthAverageAtBasex.cpp:10
BedSlopeSolutionEnum
@ BedSlopeSolutionEnum
Definition: EnumDefinitions.h:992
ISSM_MPI_Recv
int ISSM_MPI_Recv(void *buf, int count, ISSM_MPI_Datatype datatype, int source, int tag, ISSM_MPI_Comm comm, ISSM_MPI_Status *status)
Definition: issmmpi.cpp:342
EffectivePressureEnum
@ EffectivePressureEnum
Definition: EnumDefinitions.h:548
SealevelriseCumDeltathicknessOldEnum
@ SealevelriseCumDeltathicknessOldEnum
Definition: EnumDefinitions.h:689
HydrologyHeadEnum
@ HydrologyHeadEnum
Definition: EnumDefinitions.h:615
Vector::ToMPISerial
doubletype * ToMPISerial(void)
Definition: Vector.h:277
DependentObject::GetValue
IssmDouble GetValue(void)
Definition: DependentObject.cpp:102
SedimentHeadEnum
@ SedimentHeadEnum
Definition: EnumDefinitions.h:698
solutionsequence_fct
void solutionsequence_fct(FemModel *femmodel)
Definition: solutionsequence_fct.cpp:388
Vector::Copy
void Copy(Vector *to)
Definition: Vector.h:319
LoveG0Enum
@ LoveG0Enum
Definition: EnumDefinitions.h:234
EplHeadSubstepEnum
@ EplHeadSubstepEnum
Definition: EnumDefinitions.h:557
solutionsequence_adjoint_linear
void solutionsequence_adjoint_linear(FemModel *femmodel)
Definition: solutionsequence_adjoint_linear.cpp:10
MasstransportNumRequestedOutputsEnum
@ MasstransportNumRequestedOutputsEnum
Definition: EnumDefinitions.h:246
UpdateDynamicConstraintsx
void UpdateDynamicConstraintsx(Constraints *constraints, Nodes *nodes, Parameters *parameters, Vector< IssmDouble > *yg)
Definition: UpdateDynamicConstraintsx.cpp:10
GiaModelidEnum
@ GiaModelidEnum
Definition: EnumDefinitions.h:156
HydrologyshreveEnum
@ HydrologyshreveEnum
Definition: EnumDefinitions.h:1109
IssmPDouble
IssmDouble IssmPDouble
Definition: types.h:38
TransientInput2::GetTriaInput
TriaInput2 * GetTriaInput()
Definition: TransientInput2.cpp:270
Parameters::FindParamObject
Param * FindParamObject(int enum_type)
Definition: Parameters.cpp:588
DataSet
Declaration of DataSet class.
Definition: DataSet.h:14
Vector::Shift
void Shift(doubletype shift)
Definition: Vector.h:309
SteadystateReltolEnum
@ SteadystateReltolEnum
Definition: EnumDefinitions.h:401
stressbalance_core
void stressbalance_core(FemModel *femmodel)
Definition: stressbalance_core.cpp:13
TransferSealevel
void TransferSealevel(FemModel *femmodel, int forcingenum)
Definition: sealevelchange_core.cpp:879
ESACORE
#define ESACORE
Definition: Profiler.h:27
Domain2DverticalEnum
@ Domain2DverticalEnum
Definition: EnumDefinitions.h:535
BaseOldEnum
@ BaseOldEnum
Definition: EnumDefinitions.h:496
AutodiffNumIndependentsEnum
@ AutodiffNumIndependentsEnum
Definition: EnumDefinitions.h:53
FemModel::WriteMeshInResults
void WriteMeshInResults(void)
Definition: FemModel.cpp:3610
Nodes::NumberOfDofs
int NumberOfDofs(int setenum)
Definition: Nodes.cpp:314
SealevelInertiaTensorYZEnum
@ SealevelInertiaTensorYZEnum
Definition: EnumDefinitions.h:1262
FemModel::profiler
Profiler * profiler
Definition: FemModel.h:42
VerboseControl
bool VerboseControl(void)
Definition: Verbosity.cpp:27
Inputs2::AddInput
void AddInput(Input2 *in_input)
Definition: Inputs2.cpp:147
HydrologyGlaDSAnalysis::SetChannelCrossSectionOld
void SetChannelCrossSectionOld(FemModel *femmodel)
Definition: HydrologyGlaDSAnalysis.cpp:444
Vector< IssmDouble >
NGiaEnum
@ NGiaEnum
Definition: EnumDefinitions.h:591
GetStericRate
void GetStericRate(Vector< IssmDouble > **psteric_rate_g, FemModel *femmodel)
Definition: sealevelchange_core.cpp:707
HYDROLOGYCORE
#define HYDROLOGYCORE
Definition: Profiler.h:19
FemModel::ReMesh
void ReMesh(void)
Definition: FemModel.cpp:3101
FemModel::UpdateConstraintsExtrudeFromTopx
void UpdateConstraintsExtrudeFromTopx()
Definition: FemModel.cpp:3018
FrontalForcingsx
void FrontalForcingsx(FemModel *femmodel)
Definition: FrontalForcingsx.cpp:9
SealevelNEsaRateEnum
@ SealevelNEsaRateEnum
Definition: EnumDefinitions.h:679
FlowequationIsSIAEnum
@ FlowequationIsSIAEnum
Definition: EnumDefinitions.h:141
balancethickness_core
void balancethickness_core(FemModel *femmodel)
Definition: balancethickness_core.cpp:12
Matlitho::burgers_mu
IssmDouble * burgers_mu
Definition: Matlitho.h:24
EffectivePressureTransientEnum
@ EffectivePressureTransientEnum
Definition: EnumDefinitions.h:550
AggressiveMigrationEnum
@ AggressiveMigrationEnum
Definition: EnumDefinitions.h:973
SolidearthSettingsComputesealevelchangeEnum
@ SolidearthSettingsComputesealevelchangeEnum
Definition: EnumDefinitions.h:320
HydraulicPotentialOldEnum
@ HydraulicPotentialOldEnum
Definition: EnumDefinitions.h:598
GetDynamicRate
void GetDynamicRate(Vector< IssmDouble > **pdynamic_rate_g, FemModel *femmodel)
Definition: sealevelchange_core.cpp:664
Analysis
Definition: Analysis.h:30
FemModel::UpdateVertexPositionsx
int UpdateVertexPositionsx(void)
Definition: FemModel.cpp:3055
ISSM_MPI_Gather
int ISSM_MPI_Gather(void *sendbuf, int sendcnt, ISSM_MPI_Datatype sendtype, void *recvbuf, int recvcnt, ISSM_MPI_Datatype recvtype, int root, ISSM_MPI_Comm comm)
Definition: issmmpi.cpp:242
SealevelMasks
Definition: SealevelMasks.h:10
BedSlopeYEnum
@ BedSlopeYEnum
Definition: EnumDefinitions.h:501
AdjointBalancethickness2AnalysisEnum
@ AdjointBalancethickness2AnalysisEnum
Definition: EnumDefinitions.h:970
FlowequationIsHOEnum
@ FlowequationIsHOEnum
Definition: EnumDefinitions.h:139
LoveShNminEnum
@ LoveShNminEnum
Definition: EnumDefinitions.h:240
FourierLoveCorex
void FourierLoveCorex(IssmDouble *LoveKr, IssmDouble *LoveKi, IssmDouble *LoveHr, IssmDouble *LoveHi, IssmDouble *LoveLr, IssmDouble *LoveLi, IssmDouble *LoveKernelsReal, IssmDouble *LoveKernelsImag, int nfreq, IssmDouble *frequencies, int sh_nmax, int sh_nmin, IssmDouble g0, IssmDouble r0, IssmDouble mu0, bool allow_layer_deletion, int forcing_type, bool verbosemod, int numlayers, IssmDouble *radius, IssmDouble *viscosity, IssmDouble *lame_lambda, IssmDouble *lame_mu, IssmDouble *burgers_viscosity, IssmDouble *burgers_mu, IssmDouble *density, IssmDouble *isburgers, IssmDouble *issolid)
Definition: FourierLoveCorex.cpp:21
SMBpddEnum
@ SMBpddEnum
Definition: EnumDefinitions.h:1252
DatasetInput2::GetTriaInputByOffset
TriaInput2 * GetTriaInputByOffset(int i)
Definition: DatasetInput2.cpp:222
NumModelsEnum
@ NumModelsEnum
Definition: EnumDefinitions.h:276
GenericParam
Definition: GenericParam.h:26
sealevelrise_core_noneustatic
Vector< IssmDouble > * sealevelrise_core_noneustatic(FemModel *femmodel, SealevelMasks *masks, Vector< IssmDouble > *RSLg_eustatic, IssmDouble oceanarea)
Definition: sealevelchange_core.cpp:444
BasalforcingsPicoIsplumeEnum
@ BasalforcingsPicoIsplumeEnum
Definition: EnumDefinitions.h:83
SurfaceSlopeYEnum
@ SurfaceSlopeYEnum
Definition: EnumDefinitions.h:830
FlowequationIsL1L2Enum
@ FlowequationIsL1L2Enum
Definition: EnumDefinitions.h:140
TransientIsdamageevolutionEnum
@ TransientIsdamageevolutionEnum
Definition: EnumDefinitions.h:444
femmodel
FemModel * femmodel
Definition: esmfbinders.cpp:16
Generated on Thu Jul 2 2020 08:09:23 for Ice Sheet System Model by   doxygen 1.8.19