Index: /issm/trunk-jpl/src/c/CMakeLists.txt =================================================================== --- /issm/trunk-jpl/src/c/CMakeLists.txt (revision 16517) +++ /issm/trunk-jpl/src/c/CMakeLists.txt (revision 16518) @@ -72,5 +72,5 @@ $ENV{ISSM_DIR}/src/c/modules/VecMergex $ENV{ISSM_DIR}/src/c/modules/Mergesolutionfromftogx - $ENV{ISSM_DIR}/src/c/analyses + $ENV{ISSM_DIR}/src/c/cores $ENV{ISSM_DIR}/src/c/solutionsequences $ENV{ISSM_DIR}/src/c/solutionsequences) @@ -214,12 +214,12 @@ ./modules/VecMergex/VecMergex.cpp ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp - ./analyses/ProcessArguments.cpp - ./analyses/ResetBoundaryConditions.cpp - ./analyses/AnalysisConfiguration.cpp - ./analyses/WrapperCorePointerFromSolutionEnum.cpp - ./analyses/CorePointerFromSolutionEnum.cpp - ./analyses/EnvironmentInit.cpp - ./analyses/EnvironmentFinalize.cpp - ./analyses/ad_core.cpp + ./cores/ProcessArguments.cpp + ./cores/ResetBoundaryConditions.cpp + ./cores/AnalysisConfiguration.cpp + ./cores/WrapperCorePointerFromSolutionEnum.cpp + ./cores/CorePointerFromSolutionEnum.cpp + ./cores/EnvironmentInit.cpp + ./cores/EnvironmentFinalize.cpp + ./cores/ad_core.cpp ./solutionsequences/solutionsequence_linear.cpp ./solutionsequences/solutionsequence_nonlinear.cpp @@ -238,13 +238,13 @@ ./modules/AverageOntoPartitionx/AverageOntoPartitionx.cpp ./modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp - ./analyses/dakota_core.cpp - ./analyses/DakotaSpawnCore.cpp) + ./cores/dakota_core.cpp + ./cores/DakotaSpawnCore.cpp) #}}} #Transient sources {{{ set(transient_sources ./modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp - ./analyses/transient_core.cpp) + ./cores/transient_core.cpp) #}}} #Steadystate sources {{{ -set(steadystate_sources ./analyses/steadystate_core.cpp) +set(steadystate_sources ./cores/steadystate_core.cpp) #}}} #Prognostic sources {{{ @@ -253,5 +253,5 @@ ./modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp ./modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp - ./analyses/prognostic_core.cpp) + ./cores/prognostic_core.cpp) #}}} #Thermal sources {{{ @@ -271,5 +271,5 @@ ./modules/ConstraintsStatex/ThermalIsPresent.cpp ./modules/ResetConstraintsx/ThermalConstraintsReset.cpp - ./analyses/thermal_core.cpp + ./cores/thermal_core.cpp ./solutionsequences/solutionsequence_thermal_nonlinear.cpp) #}}} @@ -298,11 +298,11 @@ ./shared/Numerics/BrentSearch.cpp ./shared/Numerics/OptimalSearch.cpp - ./analyses/control_core.cpp - ./analyses/controltao_core.cpp - ./analyses/objectivefunction.cpp - ./analyses/gradient_core.cpp - ./analyses/adjointdiagnostic_core.cpp - ./analyses/adjointbalancethickness_core.cpp - ./analyses/AdjointCorePointerFromSolutionEnum.cpp + ./cores/control_core.cpp + ./cores/controltao_core.cpp + ./cores/objectivefunction.cpp + ./cores/gradient_core.cpp + ./cores/adjointdiagnostic_core.cpp + ./cores/adjointbalancethickness_core.cpp + ./cores/AdjointCorePointerFromSolutionEnum.cpp ./solutionsequences/solutionsequence_adjoint_linear.cpp) #}}} @@ -323,5 +323,5 @@ ./modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp ./modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp - ./analyses/hydrology_core.cpp + ./cores/hydrology_core.cpp ./solutionsequences/solutionsequence_hydro_nonlinear.cpp) #}}} @@ -339,5 +339,5 @@ ./modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp ./modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp - ./analyses/diagnostic_core.cpp + ./cores/diagnostic_core.cpp ./solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp) #}}} @@ -347,6 +347,6 @@ ./modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp ./modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp - ./analyses/balancethickness_core.cpp - ./analyses/dummy_core.cpp) + ./cores/balancethickness_core.cpp + ./cores/dummy_core.cpp) #}}} #Slope sources {{{ @@ -359,9 +359,9 @@ ./modules/ModelProcessorx/SurfaceSlope/CreateConstraintsSurfaceSlope.cpp ./modules/ModelProcessorx/SurfaceSlope/CreateLoadsSurfaceSlope.cpp - ./analyses/surfaceslope_core.cpp - ./analyses/bedslope_core.cpp) + ./cores/surfaceslope_core.cpp + ./cores/bedslope_core.cpp) #}}} #Gia sources {{{ -set(gia_sources ./analyses/gia_core.cpp +set(gia_sources ./cores/gia_core.cpp ./modules/ModelProcessorx/Gia/UpdateElementsGia.cpp ./modules/ModelProcessorx/Gia/CreateNodesGia.cpp Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 16517) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 16518) @@ -330,12 +330,12 @@ ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.cpp\ ./modules/Mergesolutionfromftogx/Mergesolutionfromftogx.h\ - ./analyses/ProcessArguments.cpp\ - ./analyses/ResetBoundaryConditions.cpp\ - ./analyses/AnalysisConfiguration.cpp\ - ./analyses/WrapperCorePointerFromSolutionEnum.cpp\ - ./analyses/CorePointerFromSolutionEnum.cpp\ - ./analyses/EnvironmentInit.cpp\ - ./analyses/EnvironmentFinalize.cpp\ - ./analyses/ad_core.cpp\ + ./cores/ProcessArguments.cpp\ + ./cores/ResetBoundaryConditions.cpp\ + ./cores/AnalysisConfiguration.cpp\ + ./cores/WrapperCorePointerFromSolutionEnum.cpp\ + ./cores/CorePointerFromSolutionEnum.cpp\ + ./cores/ad_core.cpp\ + ./main/EnvironmentInit.cpp\ + ./main/EnvironmentFinalize.cpp\ ./solutionsequences/solutionsequence_linear.cpp\ ./solutionsequences/solutionsequence_nonlinear.cpp\ @@ -362,14 +362,14 @@ ./modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp\ ./modules/ModelProcessorx/Dakota/UpdateElementsAndMaterialsDakota.cpp\ - ./analyses/dakota_core.cpp\ - ./analyses/DakotaSpawnCore.h\ - ./analyses/DakotaSpawnCore.cpp + ./cores/dakota_core.cpp\ + ./cores/DakotaSpawnCore.h\ + ./cores/DakotaSpawnCore.cpp #}}} #Transient sources {{{ transient_sources = ./modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp \ - ./analyses/transient_core.cpp + ./cores/transient_core.cpp #}}} #Steadystate sources {{{ -steadystate_sources = ./analyses/steadystate_core.cpp +steadystate_sources = ./cores/steadystate_core.cpp #}}} #Masstransport sources {{{ @@ -394,7 +394,7 @@ ./modules/ModelProcessorx/ExtrudeFromTop/CreateConstraintsExtrudeFromTop.cpp\ ./modules/ModelProcessorx/ExtrudeFromTop/CreateLoadsExtrudeFromTop.cpp\ - ./analyses/masstransport_core.cpp\ - ./analyses/extrudefrombase_core.cpp\ - ./analyses/extrudefromtop_core.cpp + ./cores/masstransport_core.cpp\ + ./cores/extrudefrombase_core.cpp\ + ./cores/extrudefromtop_core.cpp #}}} @@ -414,5 +414,5 @@ ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.h\ ./modules/PostprocessingEnthalpyx/PostprocessingEnthalpyx.cpp\ - ./analyses/thermal_core.cpp\ + ./cores/thermal_core.cpp\ ./solutionsequences/solutionsequence_thermal_nonlinear.cpp #}}} @@ -460,11 +460,11 @@ ./shared/Numerics/BrentSearch.cpp\ ./shared/Numerics/OptimalSearch.cpp \ - ./analyses/control_core.cpp\ - ./analyses/controltao_core.cpp\ - ./analyses/objectivefunction.cpp\ - ./analyses/gradient_core.cpp\ - ./analyses/adjointstressbalance_core.cpp\ - ./analyses/adjointbalancethickness_core.cpp\ - ./analyses/AdjointCorePointerFromSolutionEnum.cpp\ + ./cores/control_core.cpp\ + ./cores/controltao_core.cpp\ + ./cores/objectivefunction.cpp\ + ./cores/gradient_core.cpp\ + ./cores/adjointstressbalance_core.cpp\ + ./cores/adjointbalancethickness_core.cpp\ + ./cores/AdjointCorePointerFromSolutionEnum.cpp\ ./solutionsequences/solutionsequence_adjoint_linear.cpp @@ -486,5 +486,5 @@ ./modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp \ ./modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp \ - ./analyses/hydrology_core.cpp\ + ./cores/hydrology_core.cpp\ ./solutionsequences/solutionsequence_hydro_nonlinear.cpp #}}} @@ -503,5 +503,5 @@ ./modules/ModelProcessorx/StressbalanceSIA/CreateConstraintsStressbalanceSIA.cpp \ ./modules/ModelProcessorx/StressbalanceSIA/CreateLoadsStressbalanceSIA.cpp \ - ./analyses/stressbalance_core.cpp\ + ./cores/stressbalance_core.cpp\ ./solutionsequences/solutionsequence_stokescoupling_nonlinear.cpp #}}} @@ -515,7 +515,7 @@ ./modules/ModelProcessorx/Balancevelocity/CreateConstraintsBalancevelocity.cpp\ ./modules/ModelProcessorx/Balancevelocity/CreateLoadsBalancevelocity.cpp\ - ./analyses/balancethickness_core.cpp \ - ./analyses/balancevelocity_core.cpp \ - ./analyses/dummy_core.cpp + ./cores/balancethickness_core.cpp \ + ./cores/balancevelocity_core.cpp \ + ./cores/dummy_core.cpp #}}} #Slope sources {{{ @@ -528,6 +528,6 @@ ./modules/ModelProcessorx/SurfaceSlope/CreateConstraintsSurfaceSlope.cpp\ ./modules/ModelProcessorx/SurfaceSlope/CreateLoadsSurfaceSlope.cpp\ - ./analyses/surfaceslope_core.cpp\ - ./analyses/bedslope_core.cpp + ./cores/surfaceslope_core.cpp\ + ./cores/bedslope_core.cpp #}}} #MeshDeformation sources {{{ @@ -536,8 +536,8 @@ ./modules/ModelProcessorx/MeshDeformation/CreateConstraintsMeshDeformation.cpp\ ./modules/ModelProcessorx/MeshDeformation/CreateLoadsMeshDeformation.cpp\ - ./analyses/meshdeformation_core.cpp + ./cores/meshdeformation_core.cpp #}}} #Gia sources {{{ -gia_sources = ./analyses/gia_core.cpp\ +gia_sources = ./cores/gia_core.cpp\ ./modules/ModelProcessorx/Gia/UpdateElementsGia.cpp\ ./modules/ModelProcessorx/Gia/CreateNodesGia.cpp \ @@ -557,5 +557,5 @@ #}}} #Damage sources {{{ -damage_sources = ./analyses/damage_core.cpp\ +damage_sources = ./cores/damage_core.cpp\ ./modules/ModelProcessorx/Damage/UpdateElementsDamage.cpp\ ./modules/ModelProcessorx/Damage/CreateNodesDamage.cpp \ Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16517) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 16518) @@ -10,5 +10,5 @@ #include -#include "../analyses/analyses.h" +#include "../cores/cores.h" #include "../shared/io/io.h" #include "./classes.h" Index: /issm/trunk-jpl/src/c/cores/AdjointCorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/AdjointCorePointerFromSolutionEnum.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/AdjointCorePointerFromSolutionEnum.cpp (revision 16518) @@ -0,0 +1,46 @@ +/*!\file: AdjointCorePointerFromSolutionEnum.cpp + * \brief: return type of analyses, number of analyses and core solution function. + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void AdjointCorePointerFromSolutionEnum(void (**padjointcore)(FemModel*),int solutiontype){ + + /*output: */ + void (*adjointcore)(FemModel*)=NULL; + + switch(solutiontype){ + + case StressbalanceSolutionEnum: + adjointcore=&adjointstressbalance_core; + break; + case SteadystateSolutionEnum: + adjointcore=&adjointstressbalance_core; + break; + case BalancethicknessSolutionEnum: + adjointcore=&adjointbalancethickness_core; + break; + case BalancethicknessSoftSolutionEnum: + adjointcore=&dummy_core; + break; + default: + _error_("No adjoint has been implemented for solution " << EnumToStringx(solutiontype) << " yet"); + break; + } + + /*Assign output pointer:*/ + _assert_(padjointcore); + *padjointcore=adjointcore; + +} Index: /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/AnalysisConfiguration.cpp (revision 16518) @@ -0,0 +1,144 @@ +/*!\file: AnalysisConfiguration.cpp + * \brief: return type of analyses, number of analyses + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void AnalysisConfiguration(int** panalyses,int* pnumanalyses, int solutiontype){ + + /*output: */ + int numanalyses; + int* analyses=NULL; + + /*Analyses lists*/ + switch(solutiontype){ + + case StressbalanceSolutionEnum: + numanalyses=5; + analyses=xNew(numanalyses); + analyses[0]=StressbalanceAnalysisEnum; + analyses[1]=StressbalanceVerticalAnalysisEnum; + analyses[2]=StressbalanceSIAAnalysisEnum; + analyses[3]=SurfaceSlopeAnalysisEnum; + analyses[4]=BedSlopeAnalysisEnum; + break; + + case SteadystateSolutionEnum: + numanalyses=8; + analyses=xNew(numanalyses); + analyses[0]=StressbalanceAnalysisEnum; + analyses[1]=StressbalanceVerticalAnalysisEnum; + analyses[2]=StressbalanceSIAAnalysisEnum; + analyses[3]=SurfaceSlopeAnalysisEnum; + analyses[4]=BedSlopeAnalysisEnum; + analyses[5]=EnthalpyAnalysisEnum; + analyses[6]=ThermalAnalysisEnum; + analyses[7]=MeltingAnalysisEnum; + break; + + case ThermalSolutionEnum: + numanalyses=3; + analyses=xNew(numanalyses); + analyses[0]=ThermalAnalysisEnum; + analyses[1]=MeltingAnalysisEnum; + analyses[2]=EnthalpyAnalysisEnum; + break; + + case HydrologySolutionEnum: + numanalyses=5; + analyses=xNew(numanalyses); + analyses[0]=HydrologyShreveAnalysisEnum; + analyses[1]=HydrologyDCInefficientAnalysisEnum; + analyses[2]=HydrologyDCEfficientAnalysisEnum; + analyses[3]=SurfaceSlopeAnalysisEnum; + analyses[4]=BedSlopeAnalysisEnum; + break; + + case MasstransportSolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=MasstransportAnalysisEnum; + break; + + case BalancethicknessSolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=BalancethicknessAnalysisEnum; + break; + + case BalancethicknessSoftSolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=BalancethicknessAnalysisEnum; + break; + + case BalancevelocitySolutionEnum: + numanalyses=2; + analyses=xNew(numanalyses); + analyses[0]=BalancevelocityAnalysisEnum; + analyses[1]=SurfaceSlopeAnalysisEnum; + break; + + case SurfaceSlopeSolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=SurfaceSlopeAnalysisEnum; + break; + + case BedSlopeSolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=BedSlopeAnalysisEnum; + break; + + case GiaSolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=GiaAnalysisEnum; + break; + + case DamageEvolutionSolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=DamageEvolutionAnalysisEnum; + break; + + case TransientSolutionEnum: + numanalyses=13; + analyses=xNew(numanalyses); + analyses[ 0]=StressbalanceAnalysisEnum; + analyses[ 1]=StressbalanceVerticalAnalysisEnum; + analyses[ 2]=StressbalanceSIAAnalysisEnum; + analyses[ 3]=SurfaceSlopeAnalysisEnum; + analyses[ 4]=BedSlopeAnalysisEnum; + analyses[ 5]=ThermalAnalysisEnum; + analyses[ 6]=MeltingAnalysisEnum; + analyses[ 7]=EnthalpyAnalysisEnum; + analyses[ 8]=MasstransportAnalysisEnum; + analyses[ 9]=FreeSurfaceBaseAnalysisEnum; + analyses[10]=FreeSurfaceTopAnalysisEnum; + analyses[11]=ExtrudeFromBaseAnalysisEnum; + analyses[12]=ExtrudeFromTopAnalysisEnum; + break; + + default: + _error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!"); + break; + } + + /*Assign output pointers:*/ + if(pnumanalyses) *pnumanalyses=numanalyses; + if(panalyses) *panalyses=analyses; + else xDelete(analyses); +} Index: /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/CorePointerFromSolutionEnum.cpp (revision 16518) @@ -0,0 +1,127 @@ +/*!\file: CorePointerFromSolutionEnum.cpp + * \brief: return type of analyses, number of analyses and core solution function. + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype){ + + /*output: */ + void (*solutioncore)(FemModel*)=NULL; + + switch(solutiontype){ + + case StressbalanceSolutionEnum: + #ifdef _HAVE_STRESSBALANCE_ + solutioncore=&stressbalance_core; + #else + _error_("ISSM was not compiled with stressbalance capabilities. Exiting"); + #endif + break; + case SteadystateSolutionEnum: + #ifdef _HAVE_STEADYSTATE_ + solutioncore=&steadystate_core; + #else + _error_("ISSM was not compiled with steady state capabilities. Exiting"); + #endif + break; + case ThermalSolutionEnum: + #ifdef _HAVE_THERMAL_ + solutioncore=&thermal_core; + #else + _error_("ISSM was not compiled with thermal capabilities. Exiting"); + #endif + break; + case BalancethicknessSolutionEnum: + #ifdef _HAVE_BALANCED_ + solutioncore=&balancethickness_core; + #else + _error_("ISSM was not compiled with balanced capabilities. Exiting"); + #endif + break; + case BalancethicknessSoftSolutionEnum: + #ifdef _HAVE_BALANCED_ + solutioncore=&dummy_core; + #else + _error_("ISSM was not compiled with balanced capabilities. Exiting"); + #endif + break; + case BalancevelocitySolutionEnum: + #ifdef _HAVE_BALANCED_ + solutioncore=&balancevelocity_core; + #else + _error_("ISSM was not compiled with balanced capabilities. Exiting"); + #endif + break; + case HydrologySolutionEnum: + #ifdef _HAVE_HYDROLOGY_ + solutioncore=&hydrology_core; + #else + _error_("ISSM was not compiled with hydrology capabilities. Exiting"); + #endif + break; + case SurfaceSlopeSolutionEnum: + #ifdef _HAVE_SLOPE_ + solutioncore=&surfaceslope_core; + #else + _error_("ISSM was not compiled with slope capabilities. Exiting"); + #endif + break; + case BedSlopeSolutionEnum: + #ifdef _HAVE_SLOPE_ + solutioncore=&bedslope_core; + #else + _error_("ISSM was not compiled with slope capabilities. Exiting"); + #endif + break; + case TransientSolutionEnum: + #ifdef _HAVE_TRANSIENT_ + solutioncore=&transient_core; + #else + _error_("ISSM was not compiled with transient capabilities. Exiting"); + #endif + break; + case MasstransportSolutionEnum: + #ifdef _HAVE_MASSTRANSPORT_ + solutioncore=&masstransport_core; + #else + _error_("ISSM was not compiled with masstransport capabilities. Exiting"); + #endif + break; + + case GiaSolutionEnum: + #ifdef _HAVE_GIA_ + solutioncore=&gia_core; + #else + _error_("ISSM was not compiled with gia capabilities. Exiting"); + #endif + break; + case DamageEvolutionSolutionEnum: + #ifdef _HAVE_DAMAGE_ + solutioncore=&damage_core; + #else + _error_("ISSM was not compiled with damage evolution capabilities. Exiting"); + #endif + break; + + default: + _error_("solution type: " << EnumToStringx(solutiontype) << " not supported yet!"); + break; + } + + /*Assign output pointer:*/ + _assert_(psolutioncore); + *psolutioncore=solutioncore; + +} Index: /issm/trunk-jpl/src/c/cores/DakotaSpawnCore.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/DakotaSpawnCore.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/DakotaSpawnCore.cpp (revision 16518) @@ -0,0 +1,187 @@ +/*!\file: DakotaSpawnCore.cpp + + * \brief: run core ISSM solution using Dakota inputs coming from CPU 0. + * \sa qmu.cpp DakotaPlugin.cpp + * + * This routine needs to be understood simultaneously with qmu.cpp and DakotaPlugin. + * DakotaSpawnCoreParallel is called by all CPUS, with CPU 0 holding Dakota variable values, along + * with variable descriptors. + * + * DakotaSpawnCoreParallel takes care of broadcasting the variables and their descriptors across the MPI + * ring. Once this is done, we use the variables to modify the inputs for the solution core. + * For ex, if "rho_ice" is provided, for ex 920, we include "rho_ice" in the inputs, then + * call the core with the modified inputs. This is the way we get Dakota to explore the parameter + * spce of the core. + * + * Once the core is called, we process the results of the core, and using the processed results, + * we compute response functions. The responses are computed on all CPUS, but they are targeted + * for CPU 0, which will get these values back to the Dakota engine. + * + */ + +/*Includes and prototypes: {{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./cores.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +void DakotaMPI_Bcast(double** pvariables, char*** pvariables_descriptors,int* pnumvariables, int* pnumresponses); +void DakotaFree(double** pvariables,char*** pvariables_descriptors,char*** presponses_descriptors,int numvariables,int numresponses); +/*}}}*/ + +/*Notice the d_, which prefixes anything that is being provided to us by the Dakota pluggin. Careful. some things are ours, some are dakotas!: */ +int DakotaSpawnCore(double* d_responses, int d_numresponses, double* d_variables, char** d_variables_descriptors,int d_numvariables, void* void_femmodel,int counter){ + + char **responses_descriptors = NULL; //these are our! there are only numresponsedescriptors of them, not d_numresponses!!! + int numresponsedescriptors; + int solution_type; + bool control_analysis = false; + void (*solutioncore)(FemModel*) = NULL; + FemModel *femmodel = NULL; + bool nodakotacore = true; + + /*If counter==-1 on cpu0, it means that the dakota runs are done. In which case, bail out and return 0: */ + ISSM_MPI_Bcast(&counter,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + if(counter==-1)return 0; + + /*cast void_femmodel to FemModel: */ + femmodel=reinterpret_cast(void_femmodel); + + /*retrieve parameters: */ + femmodel->parameters->FindParam(&responses_descriptors,&numresponsedescriptors,QmuResponsedescriptorsEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&control_analysis,InversionIscontrolEnum); + + if(VerboseQmu()) _printf0_("qmu iteration: " << counter << "\n"); + + /* only cpu 0, running dakota is providing us with variables and variables_descriptors and numresponses: broadcast onto other cpus: */ + DakotaMPI_Bcast(&d_variables,&d_variables_descriptors,&d_numvariables,&d_numresponses); + + /*Modify core inputs in objects contained in femmodel, to reflect the dakota variables inputs: */ + InputUpdateFromDakotax(femmodel,d_variables,d_variables_descriptors,d_numvariables); + + /*Determine solution sequence: */ + if(VerboseQmu()) _printf0_("Starting " << EnumToStringx(solution_type) << " core:\n"); + WrapperCorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type,nodakotacore); + + /*Run the core solution sequence: */ + solutioncore(femmodel); + + /*compute responses: */ + if(VerboseQmu()) _printf0_("compute dakota responses:\n"); + femmodel->DakotaResponsesx(d_responses,responses_descriptors,numresponsedescriptors,d_numresponses); + + /*Free ressources:*/ + DakotaFree(&d_variables,&d_variables_descriptors,&responses_descriptors, d_numvariables, numresponsedescriptors); + + return 1; //this is critical! do not return 0, otherwise, dakota_core will stop running! +} + +void DakotaMPI_Bcast(double** pvariables, char*** pvariables_descriptors,int* pnumvariables, int* pnumresponses){ /*{{{*/ + + /* * \brief: broadcast variables_descriptors, variables, numvariables and numresponses + * from cpu 0 to all other cpus. + */ + + int i; + int my_rank; + + /*inputs and outputs: */ + double* variables=NULL; + char** variables_descriptors=NULL; + int numvariables; + int numresponses; + + /*intermediary: */ + char* string=NULL; + int string_length; + + /*recover my_rank:*/ + my_rank=IssmComm::GetRank(); + + /*recover inputs from pointers: */ + variables=*pvariables; + variables_descriptors=*pvariables_descriptors; + numvariables=*pnumvariables; + numresponses=*pnumresponses; + + /*numvariables: */ + ISSM_MPI_Bcast(&numvariables,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*variables:*/ + if(my_rank!=0)variables=xNew(numvariables); + ISSM_MPI_Bcast(variables,numvariables,MPI_DOUBLE,0,IssmComm::GetComm()); + + /*variables_descriptors: */ + if(my_rank!=0){ + variables_descriptors=xNew(numvariables); + } + for(i=0;i(string_length); + ISSM_MPI_Bcast(string,string_length,ISSM_MPI_CHAR,0,IssmComm::GetComm()); + if(my_rank!=0)variables_descriptors[i]=string; + } + + /*numresponses: */ + ISSM_MPI_Bcast(&numresponses,1,ISSM_MPI_INT,0,IssmComm::GetComm()); + + /*Assign output pointers:*/ + *pnumvariables=numvariables; + *pvariables=variables; + *pvariables_descriptors=variables_descriptors; + *pnumresponses=numresponses; +} /*}}}*/ +void DakotaFree(double** pvariables,char*** pvariables_descriptors,char*** presponses_descriptors,int numvariables,int numresponses){ /*{{{*/ + + /*\brief DakotaFree: free allocations on other cpus, not done by Dakota.*/ + + int i; + int my_rank; + + double *variables = NULL; + char **variables_descriptors = NULL; + char **responses_descriptors = NULL; + char *string = NULL; + + /*recover pointers: */ + variables=*pvariables; + variables_descriptors=*pvariables_descriptors; + responses_descriptors=*presponses_descriptors; + + /*recover my_rank:*/ + my_rank=IssmComm::GetRank(); + + /*Free variables and variables_descriptors only on cpu !=0*/ + if(my_rank!=0){ + xDelete(variables); + for(i=0;i(string); + } + xDelete(variables_descriptors); + } + + //responses descriptors on every cpu + for(i=0;i(string); + } + //rest of dynamic allocations. + xDelete(responses_descriptors); + + /*Assign output pointers:*/ + *pvariables=variables; + *pvariables_descriptors=variables_descriptors; + *presponses_descriptors=responses_descriptors; +} /*}}}*/ Index: /issm/trunk-jpl/src/c/cores/DakotaSpawnCore.h =================================================================== --- /issm/trunk-jpl/src/c/cores/DakotaSpawnCore.h (revision 16518) +++ /issm/trunk-jpl/src/c/cores/DakotaSpawnCore.h (revision 16518) @@ -0,0 +1,4 @@ +#ifndef _DAKOTA_SPAWN_CORE_ +#define _DAKOTA_SPAWN_CORE_ +int DakotaSpawnCore(double* responses, int numresponses, double* variables, char** variables_descriptors,int numvariables, void* femmodel,int counter); +#endif Index: /issm/trunk-jpl/src/c/cores/ProcessArguments.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/ProcessArguments.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/ProcessArguments.cpp (revision 16518) @@ -0,0 +1,57 @@ +/*!\file: ProcessArguments.cpp + * \brief: process arguments + */ + +#include +#include + +#include "../shared/shared.h" + +void ProcessArguments(int* solution_type,char** pbinfilename,char** poutbinfilename,char** ptoolkitsfilename,char** plockfilename,char** prootpath, int argc,char **argv){ + + char *modelname = NULL; + char *binfilename = NULL; + char *outbinfilename = NULL; + char *toolkitsfilename = NULL; + char *lockfilename = NULL; + char *rootpath = NULL; + char *rootpatharg = NULL; + + /*Check input arguments*/ + if(argc<2)_error_("Usage error: no solution requested"); + if(argc<3)_error_("Usage error: missing execution directory"); + if(argc<4)_error_("Usage error: missing model name"); + + /*Get requested solution*/ + *solution_type=StringToEnumx(argv[1]); + + rootpatharg=argv[2]; + if(strcmp(strstr(rootpatharg,"/"),"/")!=0){ + rootpath = xNew(strlen(rootpatharg)+2); sprintf(rootpath,"%s/",rootpatharg); + } + else{ + rootpath = xNew(strlen(rootpatharg)+1); sprintf(rootpath,"%s",rootpatharg); + } + + modelname=argv[3]; + if(strstr(modelname,rootpath)==NULL){ + binfilename = xNew(strlen(rootpath)+strlen(modelname)+strlen(".bin") +1); sprintf(binfilename, "%s%s%s",rootpath,modelname,".bin"); + outbinfilename = xNew(strlen(rootpath)+strlen(modelname)+strlen(".outbin")+1); sprintf(outbinfilename,"%s%s%s",rootpath,modelname,".outbin"); + toolkitsfilename = xNew(strlen(rootpath)+strlen(modelname)+strlen(".toolkits") +1); sprintf(toolkitsfilename, "%s%s%s",rootpath,modelname,".toolkits"); + lockfilename = xNew(strlen(rootpath)+strlen(modelname)+strlen(".lock") +1); sprintf(lockfilename, "%s%s%s",rootpath,modelname,".lock"); + } + else{ + binfilename = xNew(strlen(modelname)+strlen(".bin") +1); sprintf(binfilename, "%s%s",modelname,".bin"); + outbinfilename = xNew(strlen(modelname)+strlen(".outbin")+1); sprintf(outbinfilename,"%s%s",modelname,".outbin"); + toolkitsfilename = xNew(strlen(modelname)+strlen(".toolkits") +1); sprintf(toolkitsfilename, "%s%s",modelname,".toolkits"); + lockfilename = xNew(strlen(modelname)+strlen(".lock") +1); sprintf(lockfilename, "%s%s",modelname,".lock"); + } + + /*Clean up and assign output pointer*/ + *pbinfilename=binfilename; + *poutbinfilename=outbinfilename; + *ptoolkitsfilename=toolkitsfilename; + *plockfilename=lockfilename; + *prootpath=rootpath; + +} Index: /issm/trunk-jpl/src/c/cores/ResetBoundaryConditions.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/ResetBoundaryConditions.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/ResetBoundaryConditions.cpp (revision 16518) @@ -0,0 +1,27 @@ +/*!\file: ResetBoundaryConditions.cpp + * \brief: change boundary conditions of a model, using a solution vector from another analysis + */ + +#include "../classes/classes.h" +#include "../modules/modules.h" +#include "../shared/io/io.h" + +void ResetBoundaryConditions(FemModel* femmodel, int analysis_type){ + + /*variables: */ + Vector* yg = NULL; + + if(VerboseSolution()) _printf0_(" updating boundary conditions...\n"); + + /*set current analysis: */ + femmodel->SetCurrentConfiguration(analysis_type); + + /*retrieve boundary conditions from element inputs :*/ + GetSolutionFromInputsx(&yg,femmodel); + + /*update spcs using this new vector of constraints: */ + UpdateDynamicConstraintsx(femmodel->constraints,femmodel->nodes,femmodel->parameters,yg); + + /*Free ressources:*/ + delete yg; +} Index: /issm/trunk-jpl/src/c/cores/WrapperCorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/WrapperCorePointerFromSolutionEnum.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/WrapperCorePointerFromSolutionEnum.cpp (revision 16518) @@ -0,0 +1,57 @@ +/*!\file: WrapperCorePointerFromSolutionEnum.cpp + * \brief: return type of analyses, number of analyses and core solution function. + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void WrapperCorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype,bool nodakotacore){ + + /*output: */ + void (*solutioncore)(FemModel*)=NULL; + + /*parameters: */ + bool control_analysis=false; + bool tao_analysis=false; + bool dakota_analysis=false; + + /* retrieve some parameters that tell us whether wrappers are allowed, or whether we return + * a pure core. Wrappers can be dakota_core (which samples many solution_cores) or control_core (which + * carries out adjoint based inversion on a certain core: */ + parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + parameters->FindParam(&control_analysis,InversionIscontrolEnum); + parameters->FindParam(&tao_analysis,InversionTaoEnum); + + if(nodakotacore)dakota_analysis=false; + + if(dakota_analysis){ + #ifdef _HAVE_DAKOTA_ + solutioncore=dakota_core; + #else + _error_("ISSM was not compiled with dakota support, cannot carry out dakota analysis!"); + #endif + } + else if(control_analysis){ + #ifdef _HAVE_CONTROL_ + if(tao_analysis) solutioncore=controltao_core; + else solutioncore=control_core; + #else + _error_("ISSM was not compiled with control support, cannot carry out control analysis!"); + #endif + } + else CorePointerFromSolutionEnum(&solutioncore,parameters,solutiontype); /*This means we retrieve a core solution that is not a wrapper*/ + + /*Assign output pointer:*/ + _assert_(psolutioncore); + *psolutioncore=solutioncore; + +} Index: /issm/trunk-jpl/src/c/cores/ad_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/ad_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/ad_core.cpp (revision 16518) @@ -0,0 +1,356 @@ +/*!\file ad_core + * \brief: compute outputs from the AD mode, using our dependents and independents, and drivers available in Adolc. + */ + +/*Includes: {{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" +/*}}}*/ + +void ad_core(FemModel* femmodel){ + + /*diverse: */ + int i; + int dummy; + int num_dependents=0; + int num_independents=0; + bool isautodiff = false; + char *driver = NULL; + size_t tape_stats[11]; + + /*state variables: */ + IssmDouble *axp = NULL; + double *xp = NULL; + int my_rank=IssmComm::GetRank(); + + /*AD mode on?: */ + femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); + + if(isautodiff){ + + #ifdef _HAVE_ADOLC_ + if(VerboseAutodiff())_printf0_(" start ad core\n"); + + /*First, stop tracing: */ + trace_off(); + + /*retrieve parameters: */ + femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); + femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum); + + /*if no dependents, no point in running a driver: */ + if(!(num_dependents*num_independents)) return; + + /*for adolc to run in parallel, we 0 out on rank~=0:*/ + if (my_rank!=0){ + num_dependents=0; num_independents=0; + } + + /*Print tape statistics so that user can kill this run if something is off already:*/ + tapestats(my_rank,tape_stats); //reading of tape statistics + if(VerboseAutodiff()){ + int commSize=IssmComm::GetSize(); + int *sstats=new int[7]; + sstats[0]=tape_stats[NUM_OPERATIONS]; + sstats[1]=tape_stats[OP_FILE_ACCESS]; + sstats[2]=tape_stats[NUM_LOCATIONS]; + sstats[3]=tape_stats[LOC_FILE_ACCESS]; + sstats[4]=tape_stats[NUM_VALUES]; + sstats[5]=tape_stats[VAL_FILE_ACCESS]; + sstats[6]=tape_stats[TAY_STACK_SIZE]; + int *rstats=NULL; + if (my_rank==0) rstats=new int[commSize*7]; + ISSM_MPI_Gather(sstats,7,ISSM_MPI_INT,rstats,7,ISSM_MPI_INT,0,IssmComm::GetComm()); + if (my_rank==0) { + int offset=50; + int rOffset=(commSize/10)+1; + _printf_(" ADOLC statistics: \n"); + _printf_(" "<file":"") << "\n"); + _printf_(" locations: entry size " << sizeof(locint) << " Bytes\n"); + _printf_(" "<file":"") << "\n"); + _printf_(" constant values: entry size " << sizeof(double) << " Bytes\n"); + _printf_(" "<file":"") << "\n"); + _printf_(" Taylor stack: entry size " << sizeof(revreal) << " Bytes\n"); + _printf_(" "<tape_stats[TAY_BUFFER_SIZE]?" ->file":"") << "\n"); + delete []rstats; + } + delete [] sstats; + } + + /*retrieve state variable: */ + femmodel->parameters->FindParam(&axp,&dummy,AutodiffXpEnum); + + /* driver argument */ + xp=xNew(num_independents); + for(i=0;i(axp[i]); + } + + /*get the EDF pointer:*/ + ext_diff_fct *anEDF_for_solverx_p=dynamic_cast * >(femmodel->parameters->FindParamObject(AdolcParamEnum))->GetParameterValue().myEDF_for_solverx_p; + + /*Branch according to AD driver: */ + femmodel->parameters->FindParam(&driver,AutodiffDriverEnum); + + /* these are always needed regardless of the interpreter */ + anEDF_for_solverx_p->dp_x=xNew(anEDF_for_solverx_p->max_n); + anEDF_for_solverx_p->dp_y=xNew(anEDF_for_solverx_p->max_m); + + if (strcmp(driver,"fos_forward")==0){ /*{{{*/ + + int anIndepIndex; + double *tangentDir = NULL; + double *jacTimesTangentDir = NULL; + double *theOutput = NULL; + + /*retrieve direction index: */ + femmodel->parameters->FindParam(&anIndepIndex,AutodiffFosForwardIndexEnum); + + if (anIndepIndex<0 || anIndepIndex>=num_independents) _error_("index value for AutodiffFosForwardIndexEnum should be in [0,num_independents-1]"); + + tangentDir=xNewZeroInit(num_independents); + tangentDir[anIndepIndex]=1.0; + + jacTimesTangentDir=xNew(num_dependents); + theOutput=xNew(num_dependents); + + /*set the forward method function pointer: */ +#ifdef _HAVE_GSL_ + anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx; +#endif + + /*allocate the space for the parameters to invoke the EDF fos_forward:*/ + anEDF_for_solverx_p->dp_X=xNew(anEDF_for_solverx_p->max_n); + anEDF_for_solverx_p->dp_Y=xNew(anEDF_for_solverx_p->max_m); + + /*call driver: */ + fos_forward(my_rank,num_dependents,num_independents, 0, xp, tangentDir, theOutput, jacTimesTangentDir ); + + /*add to results*/ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,jacTimesTangentDir,num_dependents,1,1,0.0)); + + /*free resources :*/ + xDelete(theOutput); + xDelete(jacTimesTangentDir); + xDelete(tangentDir); + } /*}}}*/ + else if ((strcmp(driver,"fov_forward")==0) || (strcmp(driver,"fov_forward_all")==0)){ /*{{{*/ + + int tangentDirNum; + int dummy; + int *indepIndices = NULL; + double **jacTimesSeed = NULL; + double **seed = NULL; + double *theOutput = NULL; + std::set anIndexSet; + + /*retrieve directions:*/ + if (strcmp(driver,"fov_forward_all")==0){ + tangentDirNum=num_independents; + indepIndices=xNewZeroInit(tangentDirNum); + for(i=0;iparameters->FindParam(&indepIndices,&tangentDirNum,&dummy,AutodiffFovForwardIndicesEnum); + } + + /*Some checks: */ + if (tangentDirNum<1 || tangentDirNum>num_independents) _error_("tangentDirNum should be in [1,num_independents]"); + + /* full Jacobian or Jacobian projection:*/ + jacTimesSeed=xNew(num_dependents,tangentDirNum); + + /*set the forward method function pointers: */ +#ifdef _HAVE_GSL_ + anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx; +#endif + // anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; + + /*allocate the space for the parameters to invoke EDF fov_forward:*/ + anEDF_for_solverx_p->dpp_X=xNew(anEDF_for_solverx_p->max_n, tangentDirNum); + anEDF_for_solverx_p->dpp_Y=xNew(anEDF_for_solverx_p->max_m, tangentDirNum); + + /*seed matrix: */ + seed=xNewZeroInit(num_independents,tangentDirNum); + + /*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/ + for (int i=0; inum_independents) { + _error_("indepIndices values must be in [0,num_independents-1]"); + } + if (anIndexSet.find(indepIndices[i])!=anIndexSet.end()) { + _error_("duplicate indepIndices values are not allowed until we implement Jacobian decompression"); + } + anIndexSet.insert(indepIndices[i]); + /* now populate the seed matrix from the set of independent indices; + * simple setup with a single 1.0 per column and at most a single 1.0 per row*/ + seed[indepIndices[i]][i]=1.0; + } + + /*allocate output: */ + theOutput=xNew(num_dependents); + + /*call driver: */ + fov_forward(my_rank,num_dependents,num_independents, tangentDirNum, xp, seed, theOutput, jacTimesSeed ); + /*Free resources: */ + xDelete(theOutput); + xDelete(indepIndices); + xDelete(seed); + + /*add to results: */ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,*jacTimesSeed,num_dependents*tangentDirNum,1,1,0.0)); + + /*Free resources: */ + xDelete(jacTimesSeed); + xDelete(indepIndices); + } /*}}}*/ + else if (strcmp(driver,"fos_reverse")==0) { /*{{{*/ + + int aDepIndex=0; + double *aWeightVector=NULL; + double *weightVectorTimesJac=NULL; + + /*retrieve direction index: */ + femmodel->parameters->FindParam(&aDepIndex,AutodiffFosReverseIndexEnum); + aWeightVector=xNewZeroInit(num_dependents); + if (my_rank==0) { + if (aDepIndex<0 || aDepIndex>=num_dependents) _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]"); + aWeightVector[aDepIndex]=1.0; + } + weightVectorTimesJac=xNew(num_independents); + + /*set the forward method function pointer: */ +#ifdef _HAVE_GSL_ + anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx; +#endif +#ifdef _HAVE_MUMPS_ + anEDF_for_solverx_p->fos_reverse_iArr=fos_reverse_mumpsSolveEDF; +#endif + + /*allocate the space for the parameters to invoke the EDF fos_reverse :*/ + anEDF_for_solverx_p->dp_U=xNew(anEDF_for_solverx_p->max_m); + anEDF_for_solverx_p->dp_Z=xNew(anEDF_for_solverx_p->max_n); + + /*call driver: */ + fos_reverse(my_rank,num_dependents,num_independents, aWeightVector, weightVectorTimesJac ); + if(VerboseAutodiff())_printf0_(" done with fos_reverse\n"); + + /*add to results*/ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,weightVectorTimesJac,num_independents,1,1,0.0)); + + /*free resources :*/ + xDelete(weightVectorTimesJac); + xDelete(aWeightVector); + } /*}}}*/ + else if ((strcmp(driver,"fov_reverse")==0) || (strcmp(driver,"fov_reverse_all")==0)){ /*{{{*/ + + int* depIndices=NULL; + int weightNum; + int dummy; + double **weightsTimesJac=NULL; + double **weights=NULL; + std::set anIndexSet; + + /*retrieve directions:*/ + if (strcmp(driver,"fov_reverse_all")==0){ + weightNum=num_dependents; + depIndices=xNewZeroInit(weightNum); + for(i=0;iparameters->FindParam(&depIndices,&weightNum,&dummy,AutodiffFovForwardIndicesEnum); + } + + /*Some checks: */ + if (weightNum<1 || weightNum>num_dependents) _error_("tangentDirNum should be in [1,num_dependents]"); + + /* full Jacobian or Jacobian projection:*/ + weightsTimesJac=xNew(weightNum,num_independents); + + /*set the forward method function pointers: */ + #ifdef _HAVE_GSL_ + anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; + #endif + + /*allocate the space for the parameters to invoke the EDF fos_reverse :*/ + anEDF_for_solverx_p->dpp_U=xNew(weightNum,anEDF_for_solverx_p->max_m); + anEDF_for_solverx_p->dpp_Z=xNew(weightNum,anEDF_for_solverx_p->max_n); + + /*seed matrix: */ + weights=xNewZeroInit(weightNum,num_dependents); + + /*collect indices in a set to prevent accidental duplicates as long as we don't do compression:*/ + for (int i=0; inum_dependents) { + _error_("depIndices values must be in [0,num_dependents-1]"); + } + if (anIndexSet.find(depIndices[i])!=anIndexSet.end()) { + _error_("duplicate depIndices values are not allowed until we implement Jacobian decompression"); + } + anIndexSet.insert(depIndices[i]); + /* now populate the seed matrix from the set of independent indices; + * simple setup with a single 1.0 per column and at most a single 1.0 per row*/ + weights[depIndices[i]][i]=1.0; + } + + /*call driver: */ + fov_reverse(my_rank,num_dependents,num_independents, weightNum, weights, weightsTimesJac ); + + /*add to results: */ + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,AutodiffJacobianEnum,*weightsTimesJac,weightNum*num_independents,1,1,0.0)); + + /*Free resources: */ + xDelete(weights); + xDelete(weightsTimesJac); + xDelete(depIndices); + } /*}}}*/ + else _error_("driver: " << driver << " not yet supported!"); + + /* delete the allocated space for the parameters:*/ + xDelete(anEDF_for_solverx_p->dp_x); + xDelete(anEDF_for_solverx_p->dp_X); + xDelete(anEDF_for_solverx_p->dpp_X); + xDelete(anEDF_for_solverx_p->dp_y); + xDelete(anEDF_for_solverx_p->dp_Y); + xDelete(anEDF_for_solverx_p->dpp_Y); + xDelete(anEDF_for_solverx_p->dp_U); + xDelete(anEDF_for_solverx_p->dpp_U); + xDelete(anEDF_for_solverx_p->dp_Z); + xDelete(anEDF_for_solverx_p->dpp_Z); + + if(VerboseAutodiff())_printf0_(" end AD core\n"); + + /*Free resources: */ + xDelete(xp); + xDelete(axp); + xDelete(driver); + #else + _error_("Should not be requesting AD drivers when an AD library is not available!"); + #endif + } +} Index: /issm/trunk-jpl/src/c/cores/adjointbalancethickness_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/adjointbalancethickness_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/adjointbalancethickness_core.cpp (revision 16518) @@ -0,0 +1,39 @@ +/*!\file: adjointbalancethickness_core.cpp + * \brief compute inverse method adjoint state + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void adjointbalancethickness_core(FemModel* femmodel){ + + /*parameters: */ + bool save_results; + + /*retrieve parameters:*/ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + /*compute thickness */ + if(VerboseSolution()) _printf0_(" computing thickness\n"); + femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); + solutionsequence_linear(femmodel); + + /*Call SurfaceAreax, because some it might be needed by PVector*/ + SurfaceAreax(NULL,femmodel); + + /*compute adjoint*/ + if(VerboseSolution()) _printf0_(" computing adjoint\n"); + femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum,AdjointBalancethicknessAnalysisEnum); + solutionsequence_adjoint_linear(femmodel); + + /*Save results*/ + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + int outputs[1] = {AdjointEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1); + } +} Index: /issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/adjointstressbalance_core.cpp (revision 16518) @@ -0,0 +1,48 @@ +/*!\file: adjointstressbalance_core.cpp + * \brief compute inverse method adjoint state + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void adjointstressbalance_core(FemModel* femmodel){ + + /*parameters: */ + bool isFS; + bool save_results; + bool conserve_loads = true; + + /*retrieve parameters:*/ + femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + /*Compute velocities*/ + if(VerboseSolution()) _printf0_(" computing velocities\n"); + femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); + solutionsequence_nonlinear(femmodel,conserve_loads); + + /*Call SurfaceAreax, because some it might be needed by PVector*/ + SurfaceAreax(NULL,femmodel); + + /*Compute adjoint*/ + if(VerboseSolution()) _printf0_(" computing adjoint\n"); + femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum,AdjointHorizAnalysisEnum); + solutionsequence_adjoint_linear(femmodel); + + /*Save results*/ + if(save_results || true){ + if(VerboseSolution()) _printf0_(" saving results\n"); + if (isFS){ + int outputs[4] = {AdjointxEnum,AdjointyEnum,AdjointzEnum,AdjointpEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4); + } + else{ + int outputs[2] = {AdjointxEnum,AdjointyEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2); + } + } +} Index: /issm/trunk-jpl/src/c/cores/balancethickness_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/balancethickness_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/balancethickness_core.cpp (revision 16518) @@ -0,0 +1,32 @@ +/*!\file: balancethickness_core.cpp + * \brief: core of the balancethickness solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void balancethickness_core(FemModel* femmodel){ + + /*parameters: */ + bool save_results; + + /*activate formulation: */ + femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); + + /*recover parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + if(VerboseSolution()) _printf0_("call computational core:\n"); + solutionsequence_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + int outputs = ThicknessEnum; + femmodel->RequestedOutputsx(&femmodel->results,&outputs,1); + } + +} Index: /issm/trunk-jpl/src/c/cores/balancevelocity_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/balancevelocity_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/balancevelocity_core.cpp (revision 16518) @@ -0,0 +1,37 @@ +/*!\file: balancevelocity_core.cpp + * \brief: core of the balancevelocity solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void balancevelocity_core(FemModel* femmodel){ + + /*parameters: */ + bool save_results; + + /*recover parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + if(VerboseSolution()) _printf0_("computing smoothed slopes:\n"); + femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeXAnalysisEnum); + solutionsequence_linear(femmodel); + femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SmoothedSurfaceSlopeYAnalysisEnum); + solutionsequence_linear(femmodel); + //surfaceslope_core(femmodel); + + if(VerboseSolution()) _printf0_("call computational core:\n"); + femmodel->SetCurrentConfiguration(BalancevelocityAnalysisEnum); + solutionsequence_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + int outputs[3] = {SurfaceSlopeXEnum,SurfaceSlopeYEnum,VelEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3); + } + +} Index: /issm/trunk-jpl/src/c/cores/bedslope_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/bedslope_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/bedslope_core.cpp (revision 16518) @@ -0,0 +1,44 @@ +/*!\file: bedslope_core.cpp + * \brief: core of the slope solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void bedslope_core(FemModel* femmodel){ + + /*parameters: */ + bool save_results; + int meshtype; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&meshtype,MeshTypeEnum); + + if(VerboseSolution()) _printf0_(" computing slope\n"); + + /*Call on core computations: */ + femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeXAnalysisEnum); + solutionsequence_linear(femmodel); + if(meshtype!=Mesh2DverticalEnum){ + femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeYAnalysisEnum); + solutionsequence_linear(femmodel); + } + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + if(meshtype!=Mesh2DverticalEnum){ + int outputs[2] = {BedSlopeXEnum,BedSlopeYEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2); + } + else{ + int outputs[1] = {BedSlopeXEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1); + } + } + +} Index: /issm/trunk-jpl/src/c/cores/control_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/control_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/control_core.cpp (revision 16518) @@ -0,0 +1,129 @@ +/*!\file: control_core.cpp + * \brief: core of the control solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +/*Local prototypes*/ +bool controlconvergence(IssmDouble J, IssmDouble tol_cm); + +void control_core(FemModel* femmodel){ + + int i; + + /*parameters: */ + int num_controls; + int nsteps; + IssmDouble tol_cm; + int solution_type; + bool isFS; + bool dakota_analysis = false; + + int *control_type = NULL; + IssmDouble *maxiter = NULL; + IssmDouble *cm_jump = NULL; + + /*intermediary: */ + IssmDouble search_scalar = 1.; + OptArgs optargs; + OptPars optpars; + + /*Solution and Adjoint core pointer*/ + void (*solutioncore)(FemModel*) = NULL; + void (*adjointcore)(FemModel*) = NULL; + + /*output: */ + IssmDouble* J=NULL; + + /*Recover parameters used throughout the solution*/ + femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum); + femmodel->parameters->FindParam(&maxiter,NULL,InversionMaxiterPerStepEnum); + femmodel->parameters->FindParam(&cm_jump,NULL,InversionStepThresholdEnum); + femmodel->parameters->FindParam(&tol_cm,InversionCostFunctionThresholdEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); + femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + femmodel->parameters->SetParam(false,SaveResultsEnum); + + /*out of solution_type, figure out solution core and adjoint function pointer*/ + CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); + AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type); + + /*Launch once a complete solution to set up all inputs*/ + if(VerboseControl()) _printf0_(" preparing initial solution\n"); + if(isFS) solutioncore(femmodel); + + /*Initialize cost function: */ + J=xNew(nsteps); + + /*Initialize some of the BrentSearch arguments: */ + optargs.femmodel=femmodel; + optpars.xmin=0; optpars.xmax=1; + + /*Start looping: */ + for(int n=0;n(maxiter[n]); optpars.cm_jump=cm_jump[n]; + BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,&optargs); + + if(VerboseControl()) _printf0_(" updating parameter using optimized search scalar\n"); //true means update save controls + InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,true); + + if(controlconvergence(J[n],tol_cm)) break; + } + + if(VerboseControl()) _printf0_(" preparing final solution\n"); + femmodel->parameters->SetParam(true,SaveResultsEnum); + solutioncore(femmodel); + + /*some results not computed by steadystate_core or stressbalance_core: */ + if(!dakota_analysis){ //do not save this if we are running the control core from a qmu run! + femmodel->OutputControlsx(&femmodel->results); + + #ifdef _HAVE_ADOLC_ + IssmPDouble* J_passive=xNew(nsteps); + for(i=0;i(J[i]); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,J_passive,nsteps,1,1,0)); + xDelete(J_passive); + #else + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,J,nsteps,1,1,0)); + #endif + } + + /*Free ressources: */ + xDelete(control_type); + xDelete(maxiter); + xDelete(cm_jump); + xDelete(J); +} +bool controlconvergence(IssmDouble J, IssmDouble tol_cm){ + + bool converged=false; + + /*Has convergence been reached?*/ + if (!xIsNan(tol_cm) && J +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +#if defined (_HAVE_TAO_) && defined (_HAVE_PETSC_) && (_PETSC_MAJOR_ == 3 && _PETSC_MINOR_ > 1) +#include + +/*Local prototype*/ +int FormFunctionGradient(TaoSolver,Vec,IssmDouble*,Vec,void*); +int IssmMonitor(TaoSolver,void*); +typedef struct { + FemModel* femmodel; + double* J; +} AppCtx; + +void controltao_core(FemModel* femmodel){ + + /*TAO*/ + int ierr; + int num_controls,solution_type; + int nsteps,maxiter; + AppCtx user; + TaoSolver tao = 0; + int *control_list = NULL; + Vector *X = NULL; + Vector *G = NULL; + Vector *XL = NULL; + Vector *XU = NULL; + + /*Initialize TAO*/ + int argc; char **args=NULL; + PetscGetArgs(&argc,&args); + ierr = TaoInitialize(&argc,&args,(char*)0,""); + if(ierr) _error_("Could not initialize Tao"); + + /*Recover some parameters*/ + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&num_controls,InversionNumControlParametersEnum); + femmodel->parameters->FindParam(&control_list,NULL,InversionControlParametersEnum); + femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum); + femmodel->parameters->SetParam(false,SaveResultsEnum); + maxiter=nsteps*10; + + /*Prepare Toolkit*/ + ToolkitsOptionsFromAnalysis(femmodel->parameters,DefaultAnalysisEnum); + + /*Initialize TAO*/ + TaoCreate(IssmComm::GetComm(),&tao); + if(VerboseControl()) _printf0_(" Initializing the Toolkit for Advanced Optimization (TAO)\n"); + TaoSetFromOptions(tao); + TaoSetType(tao,"tao_blmvm"); + //TaoSetType(tao,"tao_cg"); + //TaoSetType(tao,"tao_lmvm"); + + /*Prepare all TAO parameters*/ + TaoSetMonitor(tao,IssmMonitor,&user,NULL); + TaoSetMaximumFunctionEvaluations(tao,maxiter); + TaoSetMaximumIterations(tao,nsteps); + TaoSetTolerances(tao,0.,0.,0.,0.,0.); + + GetVectorFromControlInputsx(&X, femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"value"); + GetVectorFromControlInputsx(&XL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"lowerbound"); + GetVectorFromControlInputsx(&XU,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,"upperbound"); + TaoSetInitialVector(tao,X->pvector->vector); + TaoSetVariableBounds(tao,XL->pvector->vector,XU->pvector->vector); + delete XL; + delete XU; + + user.J=xNewZeroInit(maxiter+5); + user.femmodel=femmodel; + TaoSetObjectiveAndGradientRoutine(tao,FormFunctionGradient,(void*)&user); + + /*Solver optimization problem*/ + if(VerboseControl()) _printf0_(" Starting optimization\n"); + TaoSolve(tao); + TaoView(tao,PETSC_VIEWER_STDOUT_WORLD); + + /*Save results*/ + TaoGetSolutionVector(tao,&X->pvector->vector); + G=new Vector(0); VecFree(&G->pvector->vector); + TaoGetGradientVector(tao,&G->pvector->vector); + SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X); + ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,G); + femmodel->OutputControlsx(&femmodel->results); + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,user.J,maxiter+3,1,1,0)); + + /*Finalize*/ + if(VerboseControl()) _printf0_(" preparing final solution\n"); + femmodel->parameters->SetParam(true,SaveResultsEnum); + void (*solutioncore)(FemModel*)=NULL; + CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); + solutioncore(femmodel); + + /*Clean up and return*/ + xDelete(control_list); + xDelete(user.J); + delete X; + TaoDestroy(&tao); + TaoFinalize(); +} +int FormFunctionGradient(TaoSolver tao, Vec Xpetsc, IssmDouble *fcn,Vec G,void *userCtx){ + + /*Retreive arguments*/ + int solution_type; + AppCtx *user = (AppCtx *)userCtx; + FemModel *femmodel = user->femmodel; + Vector *gradient = NULL; + Vector *X = NULL; + + /*Convert input to Vec*/ + X=new Vector(Xpetsc); + + /*Set new variable*/ + //VecView(X,PETSC_VIEWER_STDOUT_WORLD); + SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X); + delete X; + + /*Recover some parameters*/ + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + + /*Compute solution and adjoint*/ + void (*solutioncore)(FemModel*)=NULL; + void (*adjointcore)(FemModel*)=NULL; + CorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type); + AdjointCorePointerFromSolutionEnum(&adjointcore,solution_type); + solutioncore(femmodel); + adjointcore(femmodel); + + /*Compute objective function*/ + femmodel->CostFunctionx(fcn); + + /*Compute gradient*/ + Gradjx(&gradient,NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + VecCopy(gradient->pvector->vector,G); delete gradient; + VecScale(G,-1.); + + /*Clean-up and return*/ + return 0; +} +int IssmMonitor(TaoSolver tao, void *userCtx){ + + int its,num_responses; + IssmDouble f,gnorm,cnorm,xdiff; + AppCtx *user = (AppCtx *)userCtx; + FemModel *femmodel = user->femmodel; + int *responses = NULL; + + femmodel->parameters->FindParam(&responses,&num_responses,InversionCostFunctionsEnum); + + TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, NULL); + if(its==0) _printf0_("Iter Function Residual | List of contributions\n"); + if(its==0) _printf0_("___________________________________________________________\n"); + _printf0_(setw(4)<J[its]=f; + + /*Retrieve objective functions independently*/ + for(int i=0;iResponsex(&f,EnumToStringx(responses[i])); + _printf0_(" "<(responses); + return 0; +} + +#else +void controltao_core(FemModel* femmodel){ + _error_("TAO not installed or PETSc version not supported"); +} +#endif //_HAVE_TAO_ Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 16518) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 16518) @@ -0,0 +1,60 @@ +/* + * cores.h: + */ + +#ifndef _ANALYSES_H_ +#define _ANALYSES_H_ + +/*forward declarations: */ +struct OptArgs; +class FemModel; +class Parameters; +template class Matrix; +template class Vector; + +#include "../shared/io/Comm/IssmComm.h" +#include "../shared/Numerics/types.h" + +/*cores: */ +void adjointstressbalance_core(FemModel* femmodel); +void adjointbalancethickness_core(FemModel* femmodel); +void gradient_core(FemModel* femmodel,int n=0,bool orthogonalize=false); +void stressbalance_core(FemModel* femmodel); +void hydrology_core(FemModel* femmodel); +void thermal_core(FemModel* femmodel); +void surfaceslope_core(FemModel* femmodel); +void bedslope_core(FemModel* femmodel); +void meshdeformation_core(FemModel* femmodel); +void control_core(FemModel* femmodel); +void controltao_core(FemModel* femmodel); +void masstransport_core(FemModel* femmodel); +void extrudefrombase_core(FemModel* femmodel); +void extrudefromtop_core(FemModel* femmodel); +void balancethickness_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 dummy_core(FemModel* femmodel); +void gia_core(FemModel* femmodel); +void damage_core(FemModel* femmodel); +IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs); + +//optimization +int GradJSearch(IssmDouble* search_vector,FemModel* femmodel,int step); + +//diverse +void ProcessArguments(int* solution,char** pbinname,char** poutbinname,char** ptoolkitsname,char** plockname,char** prootpath,int argc,char **argv); +void WriteLockFile(char* filename); +void ResetBoundaryConditions(FemModel* femmodel, int analysis_type); +void PrintBanner(void); + +//solution configuration +void AnalysisConfiguration(int** pcores,int* pnumcores, int solutiontype); +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); + +#endif Index: /issm/trunk-jpl/src/c/cores/dakota_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/dakota_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/dakota_core.cpp (revision 16518) @@ -0,0 +1,230 @@ +/*!\file: dakota_core.cpp + * \brief: wrapper to the Dakota capabilities. qmu fires up Dakota, and registers a Dakota Pluggin + * which will be in charge of running the solution sequences repeteadly, to garner statistics. + * + * This routine deals with running ISSM and Dakota in library mode. In library mode, Dakota does not + * run as an execuatble. Its capabilities are linked into the ISSM software. ISSM calls dakota routines + * directly from the dakota library. qmu.cpp is the code that is in charge of calling those routines. + * + * Dakota has its own way of running in parallel (for embarassingly parallel jobs). We do not want that, + * as ISSM knows exactly how to run "really parallel" jobs that use all CPUS. To bypass Dakota's parallelism, + * we overloaded the constructor for the parallel library (see the Dakota patch in the externalpackages/dakota + * directory). This overloaded constructor fires up Dakota serially on CPU 0 only! We take care of broadcasting + * to the other CPUS, hence ISSM is running in parallel, and Dakota serially on CPU0. + * + * Now, how does CPU 0 drive all other CPUS to carry out sensitivity analysese? By synchronizing its call to + * our ISSM cores (stressbalance_core, thermal_core, transient_core, etc ...) on CPU 0 with all other CPUS. + * This explains the structure of qmu.cpp, where cpu 0 runs Dakota, the Dakota pluggin fires up DakotaSpawnCore.cpp, + * while the other CPUS are waiting for a broadcast from CPU0, once they get it, they also fire up + * DakotaSpawnCore. In the end, DakotaSpawnCore is fired up on all CPUS, with CPU0 having Dakota inputs, that it will + * broacast to other CPUS. + * + * Now, how does dakota call the DakotaSpawnCore routine? The DakotaSpawnCore is embedded into the DakotaPlugin object + * which is derived from the Direct Interface Dakota objct. This is the only way to run Dakota in library + * mode (see their developper guide for more info). Dakota registers the DakotaPlugin object into its own + * database, and calls on the embedded DakotaSpawnCore from CPU0. + * + */ + +/*include files: {{{*/ +#ifdef HAVE_CONFIG_H +#include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../shared/shared.h" +#include "../classes/classes.h" + +#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. +#if _DAKOTA_MAJOR_ < 5 || (_DAKOTA_MAJOR_ == 5 && _DAKOTA_MINOR_ < 3) +#include +#include +#include +#include +#include +#else +#include +#include +#include +#include +#include +#endif +#include "./DakotaSpawnCore.h" +#endif +/*}}}*/ +/*DakotaPlugin class {{{*/ +#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. +#if _DAKOTA_MAJOR_ < 5 || (_DAKOTA_MAJOR_ == 5 && _DAKOTA_MINOR_ < 3) +#include +#include +#include +#include +#include +#include +#else +#include +#include +#include +#include +#include +#endif +namespace SIM { + class DakotaPlugin: public Dakota::DirectApplicInterface{ + public: + /*these fields are used by core solutions: */ + void *femmodel; + int counter; + /*Constructors/Destructors*/ + DakotaPlugin(const Dakota::ProblemDescDB& problem_db,void* in_femmodel):Dakota::DirectApplicInterface(problem_db){/*{{{*/ + femmodel = in_femmodel; + counter = 0; + }/*}}}*/ + ~DakotaPlugin(){/*{{{*/ + /* Virtual destructor handles referenceCount at Interface level. */ + }/*}}}*/ + protected: + /*execute the input filter portion of a direct evaluation invocation*/ + //int derived_map_if(const Dakota::String& if_name); + /*execute an analysis code portion of a direct evaluation invocation*/ + int derived_map_ac(const Dakota::String& driver){/*{{{*/ + + int i; + IssmDouble* variables=NULL; + char** variable_descriptors=NULL; + char* variable_descriptor=NULL; + IssmDouble* responses=NULL; + + /*increae counter: */ + counter++; + + /*Before launching analysis, we need to transfer the dakota inputs into Issm + *readable variables: */ + + /*First, the variables: */ + variables=xNew(numACV); + for(i=0;i(numACV); + for(i=0;i(strlen(label.c_str())+1); + memcpy(variable_descriptor,label.c_str(),(strlen(label.c_str())+1)*sizeof(char)); + + variable_descriptors[i]=variable_descriptor; + } + + /*Initialize responses: */ + responses=xNewZeroInit(numFns); + + /*run core solution: */ + DakotaSpawnCore(responses,numFns, variables,variable_descriptors,numACV,femmodel,counter); + + /*populate responses: */ + for(i=0;i(variables); + for(i=0;i(variable_descriptor); + } + xDelete(variable_descriptors); + xDelete(responses); + + return 0; + }/*}}}*/ + /*execute the output filter portion of a direct evaluation invocation*/ + //int derived_map_of(const Dakota::String& of_name); + /*add for issm: */ + int GetCounter(){/*{{{*/ + return counter; + }/*}}}*/ + private: + }; +} +#endif +/*}}}*/ + +void dakota_core(FemModel* femmodel){ + + #ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. + + int my_rank; + char *dakota_input_file = NULL; + char *dakota_output_file = NULL; + char *dakota_error_file = NULL; + Dakota::ModelLIter ml_iter; + + /*Recover dakota_input_file, dakota_output_file and dakota_error_file, in the parameters dataset in parallel */ + femmodel->parameters->FindParam(&dakota_input_file,QmuInNameEnum); + femmodel->parameters->FindParam(&dakota_output_file,QmuOutNameEnum); + femmodel->parameters->FindParam(&dakota_error_file,QmuErrNameEnum); + + /*recover my_rank:*/ + my_rank=IssmComm::GetRank(); + + if(my_rank==0){ + + // Instantiate/initialize the parallel library and problem description + // database objects. + char* dakotamode=xNew(strlen("serial")+1); + xMemCpy(dakotamode,"serial",strlen("serial")+1); + Dakota::ParallelLibrary parallel_lib(dakotamode); //use our own ISSM Dakota library mode constructor, which only fires up Dakota on CPU 0. + Dakota::ProblemDescDB problem_db(parallel_lib); + xDelete(dakotamode); + + // Manage input file parsing, output redirection, and restart processing + // without a CommandLineHandler. This version relies on parsing of an + // input file. + problem_db.manage_inputs(dakota_input_file); + // specify_outputs_restart() is only necessary if specifying non-defaults + parallel_lib.specify_outputs_restart(dakota_output_file,dakota_error_file,NULL,NULL); + + // Instantiate the Strategy object (which instantiates all Model and + // Iterator objects) using the parsed information in problem_db. + Dakota::Strategy selected_strategy(problem_db); + + // convenience function for iterating over models and performing any + // interface plug-ins + Dakota::ModelList& models = problem_db.model_list(); + + for (ml_iter = models.begin(); ml_iter != models.end(); ml_iter++) { + + Dakota::Interface& interface = ml_iter->interface(); + + //set DB nodes to the existing Model specification + problem_db.set_db_model_nodes(ml_iter->model_id()); + + // Serial case: plug in derived Interface object without an analysisComm + interface.assign_rep(new SIM::DakotaPlugin(problem_db,(void*)femmodel), false); + } + + // Execute the strategy + problem_db.lock(); // prevent run-time DB queries + selected_strategy.run_strategy(); + + //Warn other cpus that we are done running the dakota iterator, by setting the counter to -1: + DakotaSpawnCore(NULL,0, NULL,NULL,0,femmodel,-1); + + } + else{ + + for(;;){ + if(!DakotaSpawnCore(NULL,0, NULL,NULL,0,femmodel,0))break; //counter came in at -1 on cpu0, bail out. + } + } + + /*Free ressources:*/ + xDelete(dakota_input_file); + xDelete(dakota_error_file); + xDelete(dakota_output_file); + + #endif //#ifdef _HAVE_DAKOTA_ +} Index: /issm/trunk-jpl/src/c/cores/damage_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/damage_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/damage_core.cpp (revision 16518) @@ -0,0 +1,39 @@ +/* + * \brief: damgage_core.cpp: core for the damage solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void damage_core(FemModel* femmodel){ + + /*intermediary*/ + bool save_results; + bool dakota_analysis = false; + int solution_type; + + if(VerboseSolution()) _printf0_(" computing damage\n"); + + //first recover parameters common to all solutions + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + + if(dakota_analysis && solution_type!=TransientSolutionEnum){ + femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum); + ResetConstraintsx(femmodel); + } + + femmodel->SetCurrentConfiguration(DamageEvolutionAnalysisEnum); + solutionsequence_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + int outputs = DamageDEnum; + femmodel->RequestedOutputsx(&femmodel->results,&outputs,1); + } +} Index: /issm/trunk-jpl/src/c/cores/dummy_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/dummy_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/dummy_core.cpp (revision 16518) @@ -0,0 +1,11 @@ +/*!\file: dummy_core.cpp + * \brief: dummy core (nothing done) + */ + +class FemModel; + +void dummy_core(FemModel* femmodel){ + + /*We do not do anything*/ + +} Index: /issm/trunk-jpl/src/c/cores/extrudefrombase_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/extrudefrombase_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/extrudefrombase_core.cpp (revision 16518) @@ -0,0 +1,21 @@ +/*!\file: extrudefrombase_core.cpp + * \brief: core of the extrusion solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../solutionsequences/solutionsequences.h" +#include "../modules/modules.h" + +void extrudefrombase_core(FemModel* femmodel){ + + if(VerboseSolution()) _printf0_("extruding solution from base...\n"); + + /*Call on core computations: */ + femmodel->SetCurrentConfiguration(ExtrudeFromBaseAnalysisEnum); + femmodel->UpdateConstraintsExtrudeFromBasex(); + solutionsequence_linear(femmodel); + +} Index: /issm/trunk-jpl/src/c/cores/extrudefromtop_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/extrudefromtop_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/extrudefromtop_core.cpp (revision 16518) @@ -0,0 +1,21 @@ +/*!\file: extrudefromtop_core.cpp + * \brief: core of the extrusion solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../solutionsequences/solutionsequences.h" +#include "../modules/modules.h" + +void extrudefromtop_core(FemModel* femmodel){ + + if(VerboseSolution()) _printf0_("extruding solution from top...\n"); + + /*Call on core computations: */ + femmodel->SetCurrentConfiguration(ExtrudeFromTopAnalysisEnum); + femmodel->UpdateConstraintsExtrudeFromTopx(); + solutionsequence_linear(femmodel); + +} Index: /issm/trunk-jpl/src/c/cores/gia_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/gia_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/gia_core.cpp (revision 16518) @@ -0,0 +1,56 @@ +/*!\file: gia_core.cpp + * \brief: core of the GIA solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" +void gia_core(FemModel* femmodel){ + + Vector *wg = NULL; + Vector *dwdtg = NULL; + IssmDouble *x = NULL; + IssmDouble *y = NULL; + + /*parameters: */ + bool save_results; + int gsize; + int configuration_type; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum); + + if(VerboseSolution()) _printf0_(" computing GIA\n"); + + /*Call on core computations: */ + femmodel->SetCurrentConfiguration(GiaAnalysisEnum); + + /*Figure out size of g-set deflection vector and allocate solution vector: */ + gsize = femmodel->nodes->NumberOfDofs(configuration_type,GsetEnum); + wg = new Vector(gsize); + dwdtg = new Vector(gsize); + + /*first, recover x and y vectors from vertices: */ + VertexCoordinatesx(&x,&y,NULL,femmodel->vertices); //no need for z coordinate + + /*call the main module: */ + femmodel->Deflection(wg,dwdtg,x,y); + + /*assemble vector: */ + wg->Assemble(); + dwdtg->Assemble(); + + InputUpdateFromVectorx(femmodel,wg,GiaWEnum,VertexSIdEnum); + InputUpdateFromVectorx(femmodel,dwdtg,GiadWdtEnum,VertexSIdEnum); + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + const int outputs[2] = {GiaWEnum,GiadWdtEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2); + } + +} Index: /issm/trunk-jpl/src/c/cores/gradient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/gradient_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/gradient_core.cpp (revision 16518) @@ -0,0 +1,48 @@ +/*!\file: gradient_core.cpp + * \brief compute inverse method gradient direction. + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void gradient_core(FemModel* femmodel,int step,bool orthogonalize){ + + /*Intermediaries*/ + IssmDouble norm_inf; + IssmDouble *norm_list = NULL; + Vector* new_gradient = NULL; + Vector* gradient = NULL; + Vector* old_gradient = NULL; + + /*Compute gradient*/ + if(VerboseControl()) _printf0_(" compute cost function gradient\n"); + Gradjx(&gradient,&norm_list,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters); + + if (orthogonalize){ + if(VerboseControl()) _printf0_(" orthogonalization\n"); + ControlInputGetGradientx(&old_gradient,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters); + Orthx(&new_gradient,gradient,old_gradient); delete old_gradient; delete gradient; + } + else{ + new_gradient=gradient; + } + + /*Check that gradient is clean*/ + norm_inf=new_gradient->Norm(NORM_INF); + if(norm_inf<=0) _error_("||∂J/∂α||∞ = 0 gradient norm is zero"); + if(xIsNan(norm_inf))_error_("||∂J/∂α||∞ = NaN gradient norm is NaN"); + + /*plug back into inputs: */ + ControlInputSetGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,new_gradient); + delete new_gradient; + + /*Scale Gradients*/ + ControlInputScaleGradientx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,norm_list,step); + + /*Clean up and return*/ + xDelete(norm_list); +} Index: /issm/trunk-jpl/src/c/cores/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/hydrology_core.cpp (revision 16518) @@ -0,0 +1,105 @@ +/*!\file: hydrology_core.cpp + * \brief: core of the hydrology solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void hydrology_core(FemModel* femmodel){ + + int i; + + /*intermediary*/ + int step,nsteps; + int output_frequency,hydrology_model; + bool save_results; + bool modify_loads=true; + bool isefficientlayer; + IssmDouble starttime,final_time; + IssmDouble time,dt; + + /*first recover parameters common to all solutions*/ + femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + femmodel->parameters->FindParam(&final_time,TimesteppingFinalTimeEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum); + femmodel->parameters->FindParam(&hydrology_model,HydrologyModelEnum); + + /*first compute slopes: */ + if (hydrology_model==HydrologyshreveEnum){ + surfaceslope_core(femmodel); + bedslope_core(femmodel); + } + + /*Compute number of time steps: */ + if((dt==0)|| (final_time==0)){ + dt=0; + nsteps=1; + } + else nsteps=reCast((final_time-starttime)/dt); + + /*initialize: */ + step=0; + time=starttime; + + /*Loop through time: */ + for(i=0;iparameters->SetParam(time,TimeEnum); + femmodel->parameters->SetParam(step,StepEnum); + + if (hydrology_model==HydrologyshreveEnum){ + if(VerboseSolution()) _printf0_(" computing water column\n"); + femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum); + solutionsequence_nonlinear(femmodel,modify_loads); + + /*transfer water column thickness to old water column thickness: */ + + InputDuplicatex(femmodel,WatercolumnEnum,WaterColumnOldEnum); + + if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){ + if(VerboseSolution()) _printf0_(" saving results \n"); + int outputs[3] = {WatercolumnEnum,HydrologyWaterVxEnum,HydrologyWaterVyEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3); + + /*unload results*/ + if(VerboseSolution()) _printf0_(" saving temporary results\n"); + OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results); + } + } + + else if (hydrology_model==HydrologydcEnum){ + InputDuplicatex(femmodel,SedimentHeadEnum,SedimentHeadOldEnum); + femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if (isefficientlayer){ + InputDuplicatex(femmodel,EplHeadEnum,EplHeadOldEnum); + } + + if(VerboseSolution()) _printf0_(" computing water head\n"); + solutionsequence_hydro_nonlinear(femmodel); + if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){ + if(VerboseSolution()) _printf0_(" saving results \n"); + if(isefficientlayer){ + int outputs[4] = {SedimentHeadEnum,SedimentHeadResidualEnum,EplHeadEnum,HydrologydcMaskEplactiveEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],4); + } + else{ + int outputs[2] = {SedimentHeadEnum,SedimentHeadResidualEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2); + } + /*unload results*/ + if(VerboseSolution()) _printf0_(" saving temporary results\n"); + OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results); + } + } + + } +} Index: /issm/trunk-jpl/src/c/cores/masstransport_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/masstransport_core.cpp (revision 16518) @@ -0,0 +1,87 @@ +/*!\file: masstransport_core.cpp + * \brief: core of the masstransport solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void masstransport_core(FemModel* femmodel){ + + /*parameters: */ + int i; + int numoutputs,meshtype; + bool save_results; + bool issmbgradients,ispdd,isdelta18o,isFS,isfreesurface,dakota_analysis; + int solution_type; + char** requested_outputs = NULL; + + /*activate configuration*/ + femmodel->SetCurrentConfiguration(MasstransportAnalysisEnum); + + /*recover parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&issmbgradients,SurfaceforcingsIssmbgradientsEnum); + femmodel->parameters->FindParam(&ispdd,SurfaceforcingsIspddEnum); + femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); + femmodel->parameters->FindParam(&isfreesurface,MasstransportIsfreesurfaceEnum); + femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&meshtype,MeshTypeEnum); + femmodel->parameters->FindParam(&numoutputs,MasstransportNumRequestedOutputsEnum); + if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,MasstransportRequestedOutputsEnum); + + /*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/ + if(dakota_analysis && solution_type==MasstransportSolutionEnum){ + InputDuplicatex(femmodel,QmuSurfaceEnum,SurfaceEnum); + InputDuplicatex(femmodel,QmuThicknessEnum,ThicknessEnum); + InputDuplicatex(femmodel,QmuBedEnum,BedEnum); + } + + if(issmbgradients){ + if(VerboseSolution())_printf_(" call smb gradients module\n\n"); + SmbGradientsx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + if(ispdd){ + if(isdelta18o){ + if(VerboseSolution()) _printf0_(" call Delta18oParametrization module\n"); + Delta18oParameterizationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + if(VerboseSolution()) _printf0_(" call positive degree day module\n"); + PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + + if(isFS && isfreesurface){ + if(VerboseSolution()) _printf0_(" call free surface computational core\n"); + femmodel->SetCurrentConfiguration(FreeSurfaceBaseAnalysisEnum); + solutionsequence_linear(femmodel); + if(meshtype==Mesh2DverticalEnum){ + femmodel->parameters->SetParam(BedEnum,InputToExtrudeEnum); + extrudefrombase_core(femmodel); + } + femmodel->SetCurrentConfiguration(FreeSurfaceTopAnalysisEnum); + solutionsequence_linear(femmodel); + if(meshtype==Mesh2DverticalEnum){ + femmodel->parameters->SetParam(SurfaceEnum,InputToExtrudeEnum); + extrudefromtop_core(femmodel); + } + } + else{ + if(VerboseSolution()) _printf0_(" call computational core\n"); + solutionsequence_linear(femmodel); + } + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + } + + if(solution_type==MasstransportSolutionEnum)femmodel->RequestedDependentsx(); + + /*Free ressources:*/ + if(numoutputs){for(int i=0;i(string);} xDelete(requested_outputs);} +} Index: /issm/trunk-jpl/src/c/cores/meshdeformation_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/meshdeformation_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/meshdeformation_core.cpp (revision 16518) @@ -0,0 +1,30 @@ +/*!\file: meshdeformation_core.cpp + * \brief: core of the slope solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../solutionsequences/solutionsequences.h" +#include "../modules/modules.h" + +void meshdeformation_core(FemModel* femmodel){ + + /*parameters: */ + bool save_results; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + if(VerboseSolution()) _printf0_("computing mesh deformation (elasticity)...\n"); + + /*Call on core computations: */ + femmodel->SetCurrentConfiguration(MeshdeformationAnalysisEnum); + solutionsequence_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _printf0_("saving results:\n"); + } + +} Index: /issm/trunk-jpl/src/c/cores/objectivefunction.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/objectivefunction.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/objectivefunction.cpp (revision 16518) @@ -0,0 +1,77 @@ +/*!\file: objectivefunction + * \brief objective function that returns a misfit, for a certain parameter. + */ + +/*include files: {{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../solutionsequences/solutionsequences.h" +#include "../modules/modules.h" +/*}}}*/ + +IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs){ + + /*output: */ + IssmDouble J; + + /*parameters: */ + int solution_type,analysis_type; + bool isFS = false; + bool conserve_loads = true; + FemModel *femmodel = NULL; + + /*Recover finite element model: */ + femmodel=optargs->femmodel; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); + femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + + /*set analysis type to compute velocity: */ + if (solution_type==SteadystateSolutionEnum || solution_type==StressbalanceSolutionEnum){ + femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); + } + else if (solution_type==BalancethicknessSolutionEnum){ + femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); + } + else if (solution_type==BalancethicknessSoftSolutionEnum){ + femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); + } + else{ + _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet"); + } + + /*update parameter according to scalar: */ //false means: do not save control + InputControlUpdatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,search_scalar,false); + + /*Run stressbalance with updated inputs: */ + if (solution_type==SteadystateSolutionEnum){ + stressbalance_core(femmodel); //We need a 3D velocity!! (vz is required for the next thermal run) + } + else if (solution_type==StressbalanceSolutionEnum){ + solutionsequence_nonlinear(femmodel,conserve_loads); + } + else if (solution_type==BalancethicknessSolutionEnum){ + solutionsequence_linear(femmodel); + } + else if (solution_type==BalancethicknessSoftSolutionEnum){ + /*Don't do anything*/ + } + else{ + _error_("Solution " << EnumToStringx(solution_type) << " not implemented yet"); + } + + /*Compute misfit for this velocity field.*/ + femmodel->CostFunctionx(&J); + + /*Free ressources:*/ + return J; +} Index: /issm/trunk-jpl/src/c/cores/steadystate_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/steadystate_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/steadystate_core.cpp (revision 16518) @@ -0,0 +1,130 @@ +/*!\file: steadystate_core.cpp + * \brief: core of the steadystate solution + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +/*Local prototypes*/ +bool steadystateconvergence(Vector* tg,Vector* tg_old,Vector* ug,Vector* ug_old,IssmDouble reltol); + +void steadystate_core(FemModel* femmodel){ + + /*intermediary: */ + int i; + int step; + Vector* ug = NULL; + Vector* ug_old = NULL; + Vector* tg = NULL; + Vector* tg_old = NULL; + + /*parameters: */ + bool save_results,isenthalpy; + int maxiter; + IssmDouble reltol; + int numoutputs = 0; + char** requested_outputs = NULL; + + /* recover parameters:*/ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&maxiter,SteadystateMaxiterEnum); + femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); + femmodel->parameters->FindParam(&reltol,SteadystateReltolEnum); + femmodel->parameters->SetParam(false,SaveResultsEnum); + femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum); + if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SteadystateRequestedOutputsEnum); + + /*intialize counters: */ + step=1; + + for(;;){ + + if(VerboseSolution()) _printf0_(" computing temperature and velocity for step: " << step << "\n"); + #ifdef _HAVE_THERMAL_ + thermal_core(femmodel); + if(!isenthalpy)femmodel->SetCurrentConfiguration(ThermalAnalysisEnum);/*Could be MeltingAnalysis...*/ + GetSolutionFromInputsx(&tg,femmodel); + #else + _error_("ISSM was not compiled with thermal capabilities. Exiting"); + #endif + + if(VerboseSolution()) _printf0_(" computing new velocity\n"); + stressbalance_core(femmodel); + GetSolutionFromInputsx(&ug,femmodel); + + if(step>1){ + if(VerboseSolution()) _printf0_(" checking steadystate convergence\n"); + if(steadystateconvergence(tg,tg_old,ug,ug_old,reltol)) break; + } + if(step>maxiter){ + if(VerboseSolution()) _printf0_(" maximum number steadystate iterations " << maxiter << " reached\n"); + break; + } + + /*update results and increase counter*/ + delete tg_old;tg_old=tg; + delete ug_old;ug_old=ug; + step++; + } + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + } + + /*Free ressources:*/ + delete tg_old; + delete ug_old; + delete tg; + delete ug; + if(numoutputs){ for (i=0;i(string);} xDelete(requested_outputs);} +} +bool steadystateconvergence(Vector* tg,Vector* tg_old,Vector* ug,Vector* ug_old,IssmDouble reltol){ + + /*Output*/ + bool converged = true; + + /*Intermediary*/ + Vector* dug = NULL; + Vector* dtg = NULL; + IssmDouble ndt,nt; + IssmDouble ndu,nu; + + /*compute norm(du)/norm(u)*/ + dug=ug_old->Duplicate(); ug_old->Copy(dug); dug->AYPX(ug,-1.0); + ndu=dug->Norm(NORM_TWO); nu=ug_old->Norm(NORM_TWO); + if (xIsNan(ndu) || xIsNan(nu)) _error_("convergence criterion is NaN!"); + if((ndu/nu) "<Duplicate(); tg_old->Copy(dtg); dtg->AYPX(tg,-1.0); + ndt=dtg->Norm(NORM_TWO); nt=tg_old->Norm(NORM_TWO); + if (xIsNan(ndt) || xIsNan(nt)) _error_("convergence criterion is NaN!"); + if((ndt/nt) "<parameters->FindParam(&meshtype,MeshTypeEnum); + femmodel->parameters->FindParam(&isSIA,FlowequationIsSIAEnum); + femmodel->parameters->FindParam(&isSSA,FlowequationIsSSAEnum); + femmodel->parameters->FindParam(&isL1L2,FlowequationIsL1L2Enum); + femmodel->parameters->FindParam(&isHO,FlowequationIsHOEnum); + femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); + femmodel->parameters->FindParam(&newton,StressbalanceIsnewtonEnum); + femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&numoutputs,StressbalanceNumRequestedOutputsEnum); + if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,StressbalanceRequestedOutputsEnum); + + /*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/ + if(dakota_analysis && solution_type==StressbalanceSolutionEnum){ + InputDuplicatex(femmodel,QmuVxEnum,VxEnum); + InputDuplicatex(femmodel,QmuVyEnum,VyEnum); + InputDuplicatex(femmodel,QmuVzEnum,VzEnum); + if(isFS) InputDuplicatex(femmodel,QmuPressureEnum,PressureEnum); + } + + /*Compute slopes: */ + if(isSIA) surfaceslope_core(femmodel); + if(isFS){ + bedslope_core(femmodel); + femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); + ResetCoordinateSystemx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + + if(isSIA){ + if(VerboseSolution()) _printf0_(" computing SIA velocities\n"); + + /*Take the last velocity into account so that the velocity on the SSA domain is not zero*/ + if(isSSA || isL1L2 || isHO ) ResetBoundaryConditions(femmodel,StressbalanceSIAAnalysisEnum); + femmodel->SetCurrentConfiguration(StressbalanceSIAAnalysisEnum); + solutionsequence_linear(femmodel); + if(isSSA || isL1L2 || isHO) ResetBoundaryConditions(femmodel,StressbalanceAnalysisEnum); + } + + if ((isSSA || isHO || isL1L2) ^ isFS){ // ^ = xor + if(VerboseSolution()) _printf0_(" computing velocities\n"); + + femmodel->SetCurrentConfiguration(StressbalanceAnalysisEnum); + if(newton>0) + solutionsequence_newton(femmodel); + else + solutionsequence_nonlinear(femmodel,conserve_loads); + } + + if ((isSSA || isL1L2 || isHO) && isFS){ + if(VerboseSolution()) _printf0_(" computing coupling betweem lower order models and full-FS\n"); + solutionsequence_FScoupling_nonlinear(femmodel,conserve_loads); + } + + if (meshtype==Mesh3DEnum && (isSIA || isSSA || isL1L2 || isHO)){ + if(VerboseSolution()) _printf0_(" computing vertical velocities\n"); + femmodel->SetCurrentConfiguration(StressbalanceVerticalAnalysisEnum); + solutionsequence_linear(femmodel); + } + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + } + + if(solution_type==StressbalanceSolutionEnum)femmodel->RequestedDependentsx(); + + /*Free ressources:*/ + if(numoutputs){ for (i=0;i(string);} xDelete(requested_outputs);} +} Index: /issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/surfaceslope_core.cpp (revision 16518) @@ -0,0 +1,43 @@ +/*!\file: surfaceslope_core.cpp + * \brief: core of the slope solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../solutionsequences/solutionsequences.h" +#include "../modules/modules.h" + +void surfaceslope_core(FemModel* femmodel){ + + /*parameters: */ + bool save_results; + int meshtype; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&meshtype,MeshTypeEnum); + + if(VerboseSolution()) _printf0_("computing slope...\n"); + + /*Call on core computations: */ + femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum); + solutionsequence_linear(femmodel); + femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum); + solutionsequence_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _printf0_("saving results:\n"); + if(meshtype!=Mesh2DverticalEnum){ + int outputs[2] = {SurfaceSlopeXEnum,SurfaceSlopeYEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2); + + } + else{ + int outputs = SurfaceSlopeXEnum; + femmodel->RequestedOutputsx(&femmodel->results,&outputs,1); + } + } + +} Index: /issm/trunk-jpl/src/c/cores/thermal_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/thermal_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/thermal_core.cpp (revision 16518) @@ -0,0 +1,65 @@ +/*!\file: thermal_core.cpp + * \brief: core of the thermal solution + */ + +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void thermal_core(FemModel* femmodel){ + + /*intermediary*/ + bool save_results,isenthalpy; + bool dakota_analysis; + int solution_type,numoutputs; + char** requested_outputs = NULL; + + /*first recover parameters common to all solutions*/ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); + femmodel->parameters->FindParam(&numoutputs,ThermalNumRequestedOutputsEnum); + if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,ThermalRequestedOutputsEnum); + + if(dakota_analysis && solution_type!=TransientSolutionEnum){ + InputDuplicatex(femmodel,QmuVxMeshEnum,VxMeshEnum); + InputDuplicatex(femmodel,QmuVyMeshEnum,VyMeshEnum); + InputDuplicatex(femmodel,QmuVzMeshEnum,VzMeshEnum); + InputDuplicatex(femmodel,QmuTemperatureEnum,TemperatureEnum); + InputDuplicatex(femmodel,QmuMeltingEnum,BasalforcingsMeltingRateEnum); + InputDuplicatex(femmodel,QmuMaterialsRheologyBEnum,MaterialsRheologyBEnum); + femmodel->SetCurrentConfiguration(ThermalAnalysisEnum); + ResetConstraintsx(femmodel); + } + + if(isenthalpy){ + + if(VerboseSolution()) _printf0_(" computing enthalpy\n"); + femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum); + solutionsequence_nonlinear(femmodel,true); + + /*transfer enthalpy to enthalpy picard for the next step: */ + InputDuplicatex(femmodel,EnthalpyEnum,EnthalpyPicardEnum); + + /*Post process*/ + if(solution_type!=SteadystateSolutionEnum) PostprocessingEnthalpyx(femmodel); + } + else{ + if(VerboseSolution()) _printf0_(" computing temperatures\n"); + femmodel->SetCurrentConfiguration(ThermalAnalysisEnum); + solutionsequence_thermal_nonlinear(femmodel); + + if(VerboseSolution()) _printf0_(" computing melting\n"); + femmodel->SetCurrentConfiguration(MeltingAnalysisEnum); + solutionsequence_linear(femmodel); + } + + if(save_results){ + if(VerboseSolution()) _printf0_(" saving results\n"); + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + } +} Index: /issm/trunk-jpl/src/c/cores/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/cores/transient_core.cpp (revision 16518) @@ -0,0 +1,182 @@ +/*!\file: transient_3d_core.cpp + * \brief: core of the transient_3d solution + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include +#include "./cores.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" +#include "../solutionsequences/solutionsequences.h" + +void transient_core(FemModel* femmodel){ + + /*parameters: */ + int i; + IssmDouble starttime,finaltime,dt,yts; + bool isstressbalance,ismasstransport,isFS,isthermal,isgroundingline,isdelta18o,isgia; + bool save_results,dakota_analysis; + bool time_adapt=false; + int output_frequency; + int meshtype,groundingline_migration; + int numoutputs = 0; + char** requested_outputs = NULL; + + + /*intermediary: */ + int step; + IssmDouble time; + + //first recover parameters common to all solutions + femmodel->parameters->FindParam(&meshtype,MeshTypeEnum); + femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + femmodel->parameters->FindParam(&finaltime,TimesteppingFinalTimeEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&yts,ConstantsYtsEnum); + femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + femmodel->parameters->FindParam(&output_frequency,SettingsOutputFrequencyEnum); + femmodel->parameters->FindParam(&time_adapt,TimesteppingTimeAdaptEnum); + femmodel->parameters->FindParam(&isstressbalance,TransientIsstressbalanceEnum); + femmodel->parameters->FindParam(&ismasstransport,TransientIsmasstransportEnum); + femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum); + femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum); + femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum); + femmodel->parameters->FindParam(&isFS,FlowequationIsFSEnum); + if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum); + femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum); + if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,TransientRequestedOutputsEnum); + femmodel->parameters->FindParam(&isdelta18o,SurfaceforcingsIsdelta18oEnum); + + /*initialize: */ + step=0; + time=starttime; + + /*for qmu analysis, reinitialize velocity so that fake sensitivities do not show up as a result of a different restart of the convergence at each trial.*/ + if(dakota_analysis){ + if(isstressbalance){ + InputDuplicatex(femmodel,QmuVxEnum,VxEnum); + InputDuplicatex(femmodel,QmuVyEnum,VyEnum); + if(meshtype==Mesh3DEnum){ + InputDuplicatex(femmodel,QmuVzEnum,VzEnum); + if(isFS)InputDuplicatex(femmodel,QmuPressureEnum,PressureEnum); + } + } + if(ismasstransport || isgroundingline){ + InputDuplicatex(femmodel,QmuThicknessEnum,ThicknessEnum); + InputDuplicatex(femmodel,QmuSurfaceEnum,SurfaceEnum); + InputDuplicatex(femmodel,QmuBedEnum,BedEnum); + InputDuplicatex(femmodel,QmuMaskIceLevelsetEnum,MaskIceLevelsetEnum); + } + if(isgroundingline) InputDuplicatex(femmodel,QmuMaskGroundediceLevelsetEnum,MaskGroundediceLevelsetEnum); + if(isthermal && meshtype==Mesh3DEnum){ + //Update Vertex Position after updating Thickness and Bed + femmodel->SetCurrentConfiguration(MasstransportAnalysisEnum); + femmodel->UpdateVertexPositionsx(); + InputDuplicatex(femmodel,QmuVxMeshEnum,VxMeshEnum); + InputDuplicatex(femmodel,QmuVyMeshEnum,VyMeshEnum); + InputDuplicatex(femmodel,QmuVzMeshEnum,VzMeshEnum); + InputDuplicatex(femmodel,QmuTemperatureEnum,TemperatureEnum); + InputDuplicatex(femmodel,QmuMeltingEnum,BasalforcingsMeltingRateEnum); + InputDuplicatex(femmodel,QmuMaterialsRheologyBEnum,MaterialsRheologyBEnum); + //Reset Thermal Constraints + femmodel->SetCurrentConfiguration(ThermalAnalysisEnum); + ResetConstraintsx(femmodel); + } + } + + while(time < finaltime - (yts*DBL_EPSILON)){ //make sure we run up to finaltime. + + /*Increment*/ + if(time_adapt){ + femmodel->TimeAdaptx(&dt); + if(time+dt>finaltime) dt=finaltime-time; + femmodel->parameters->SetParam(dt,TimesteppingTimeStepEnum); + } + step+=1; + time+=dt; + femmodel->parameters->SetParam(time,TimeEnum); + femmodel->parameters->SetParam(step,StepEnum); + + if(VerboseSolution()) _printf0_("iteration " << step << "/" << floor((finaltime-time)/dt)+step << " time [yr]: " << time/yts << " (time step: " << dt/yts << ")\n"); + if(step%output_frequency==0 || time==finaltime) + save_results=true; + else + save_results=false; + femmodel->parameters->SetParam(save_results,SaveResultsEnum); + + if(isthermal && meshtype==Mesh3DEnum){ + if(VerboseSolution()) _printf0_(" computing temperatures\n"); + #ifdef _HAVE_THERMAL_ + thermal_core(femmodel); + #else + _error_("ISSM was not compiled with thermal capabilities. Exiting"); + #endif + } + + if(isstressbalance){ + if(VerboseSolution()) _printf0_(" computing new velocity\n"); + #ifdef _HAVE_STRESSBALANCE_ + stressbalance_core(femmodel); + #else + _error_("ISSM was not compiled with stressbalance capabilities. Exiting"); + #endif + } + + if(ismasstransport){ + if(VerboseSolution()) _printf0_(" computing new thickness\n"); + masstransport_core(femmodel); + if(VerboseSolution()) _printf0_(" updating vertices positions\n"); + femmodel->UpdateVertexPositionsx(); + } + + if(isgroundingline){ + if(VerboseSolution()) _printf0_(" computing new grounding line position\n"); + #ifdef _HAVE_GROUNDINGLINE_ + GroundinglineMigrationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + #else + _error_("ISSM was not compiled with grounding line migration capabilities. Exiting"); + #endif + if(save_results){ + int outputs[3] = {SurfaceEnum,BedEnum,MaskGroundediceLevelsetEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],3); + } + } + if(isgia){ + if(VerboseSolution()) _printf0_(" computing glacial isostatic adjustment\n"); + #ifdef _HAVE_GIA_ + gia_core(femmodel); + #else + _error_("ISSM was not compiled with gia capabilities. Exiting"); + #endif + + } + + /*unload results*/ + if(save_results){ + if(VerboseSolution()) _printf0_(" saving transient results\n"); + femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs); + if(isdelta18o){ + int outputs[2] = {SurfaceforcingsMonthlytemperaturesEnum,SurfaceforcingsPrecipitationEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],2); + } + if(isgroundingline && (groundingline_migration==SubelementMigrationEnum || groundingline_migration==SubelementMigration2Enum)){ + int outputs[1] = {MaskGroundediceLevelsetEnum}; + femmodel->RequestedOutputsx(&femmodel->results,&outputs[0],1); + } + if(VerboseSolution()) _printf0_(" saving temporary results\n"); + OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results); + } + } + + femmodel->RequestedDependentsx(); + + /*Free ressources:*/ + if(numoutputs){ for (i=0;i(string);} xDelete(requested_outputs);} +} Index: /issm/trunk-jpl/src/c/main/EnvironmentFinalize.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/EnvironmentFinalize.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/main/EnvironmentFinalize.cpp (revision 16518) @@ -0,0 +1,25 @@ +/*!\file: EnvironmentFinalize.cpp + * \brief: finalize Petsc, MPI, you name it + */ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif +#include "../toolkits/toolkits.h" +#include "../shared/shared.h" + +void EnvironmentFinalize(void){ + + int my_rank; + + /*Make sure we are all here*/ + ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); + + /*Print closing statement*/ + ISSM_MPI_Comm_rank(ISSM_MPI_COMM_WORLD,&my_rank); + if(!my_rank) printf("closing MPI\n"); + + /*Finalize: */ + ISSM_MPI_Finalize(); +} Index: /issm/trunk-jpl/src/c/main/EnvironmentInit.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/EnvironmentInit.cpp (revision 16518) +++ /issm/trunk-jpl/src/c/main/EnvironmentInit.cpp (revision 16518) @@ -0,0 +1,35 @@ +/*!\file: EnvironmentInit.cpp + * \brief: initialize Petsc, MPI, you name it + */ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif +#include +#include "../toolkits/toolkits.h" + +ISSM_MPI_Comm EnvironmentInit(int argc,char** argv){ + + /*Output*/ + ISSM_MPI_Comm comm = 0; + + /*Initialize MPI environment: */ + #if defined(_HAVE_MPI_) + ISSM_MPI_Init(&argc,&argv); + comm = ISSM_MPI_COMM_WORLD; + #else + comm = 1; //bogus number for comm, which does not exist anyway. + #endif + + /*Print Banner*/ + int my_rank = 0; + ISSM_MPI_Comm_rank(comm,&my_rank); + if(!my_rank) printf("\n"); + if(!my_rank) printf("%s version %s\n",PACKAGE_NAME,PACKAGE_VERSION); + if(!my_rank) printf("(website: %s contact: %s)\n",PACKAGE_URL,PACKAGE_BUGREPORT); + if(!my_rank) printf("\n"); + + /*Return communicator*/ + return comm; +} Index: /issm/trunk-jpl/src/c/main/issm.h =================================================================== --- /issm/trunk-jpl/src/c/main/issm.h (revision 16517) +++ /issm/trunk-jpl/src/c/main/issm.h (revision 16518) @@ -16,7 +16,11 @@ #include "../classes/classes.h" #include "../toolkits/toolkits.h" -#include "../analyses/analyses.h" +#include "../cores/cores.h" #include "../solutionsequences/solutionsequences.h" #include "../modules/modules.h" +/*Environment*/ +ISSM_MPI_Comm EnvironmentInit(int argc,char** argv); +void EnvironmentFinalize(void); + #endif //ifndef _ISSM_H_