Index: /issm/trunk-jpl/src/c/Makefile.am =================================================================== --- /issm/trunk-jpl/src/c/Makefile.am (revision 14984) +++ /issm/trunk-jpl/src/c/Makefile.am (revision 14985) @@ -322,13 +322,13 @@ ./modules/InputConvergencex/InputConvergencex.cpp\ ./modules/InputConvergencex/InputConvergencex.h\ - ./solutions/convergence.cpp\ - ./solutions/ProcessArguments.cpp\ - ./solutions/ResetBoundaryConditions.cpp\ - ./solutions/AnalysisConfiguration.cpp\ - ./solutions/WrapperCorePointerFromSolutionEnum.cpp\ - ./solutions/CorePointerFromSolutionEnum.cpp\ - ./solutions/EnvironmentInit.cpp\ - ./solutions/EnvironmentFinalize.cpp\ - ./solutions/ad_core.cpp\ + ./analyses/convergence.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\ ./solvers/solver_linear.cpp\ ./solvers/solver_nonlinear.cpp\ @@ -353,15 +353,15 @@ ./modules/AverageOntoPartitionx/AverageOntoPartitionx.h\ ./modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp\ - ./solutions/dakota_core.cpp\ - ./solutions/DakotaSpawnCore.h\ - ./solutions/DakotaSpawnCore.cpp + ./analyses/dakota_core.cpp\ + ./analyses/DakotaSpawnCore.h\ + ./analyses/DakotaSpawnCore.cpp #}}} #Transient sources {{{ transient_sources = ./modules/ModelProcessorx/Transient/UpdateElementsTransient.cpp \ - ./solutions/transient_core.cpp + ./analyses/transient_core.cpp #}}} #Steadystate sources {{{ -steadystate_sources = ./solutions/steadystate_core.cpp\ - ./solutions/steadystateconvergence.cpp +steadystate_sources = ./analyses/steadystate_core.cpp\ + ./analyses/steadystateconvergence.cpp #}}} #Prognostic sources {{{ @@ -370,5 +370,5 @@ ./modules/ModelProcessorx/Prognostic/CreateConstraintsPrognostic.cpp\ ./modules/ModelProcessorx/Prognostic/CreateLoadsPrognostic.cpp\ - ./solutions/prognostic_core.cpp + ./analyses/prognostic_core.cpp #}}} #Thermal sources {{{ @@ -388,6 +388,6 @@ ./modules/ConstraintsStatex/ThermalIsPresent.cpp\ ./modules/ResetConstraintsx/ThermalConstraintsReset.cpp \ - ./solutions/thermal_core.cpp\ - ./solutions/enthalpy_core.cpp\ + ./analyses/thermal_core.cpp\ + ./analyses/enthalpy_core.cpp\ ./solvers/solver_thermal_nonlinear.cpp #}}} @@ -435,13 +435,13 @@ ./shared/Numerics/BrentSearch.cpp\ ./shared/Numerics/OptimalSearch.cpp \ - ./solutions/control_core.cpp\ - ./solutions/controltao_core.cpp\ - ./solutions/controlrestart.cpp\ - ./solutions/controlconvergence.cpp\ - ./solutions/objectivefunction.cpp\ - ./solutions/gradient_core.cpp\ - ./solutions/adjointdiagnostic_core.cpp\ - ./solutions/adjointbalancethickness_core.cpp\ - ./solutions/AdjointCorePointerFromSolutionEnum.cpp\ + ./analyses/control_core.cpp\ + ./analyses/controltao_core.cpp\ + ./analyses/controlrestart.cpp\ + ./analyses/controlconvergence.cpp\ + ./analyses/objectivefunction.cpp\ + ./analyses/gradient_core.cpp\ + ./analyses/adjointdiagnostic_core.cpp\ + ./analyses/adjointbalancethickness_core.cpp\ + ./analyses/AdjointCorePointerFromSolutionEnum.cpp\ ./solvers/solver_adjoint_linear.cpp @@ -463,6 +463,6 @@ ./modules/ModelProcessorx/HydrologyDCEfficient/CreateLoadsHydrologyDCEfficient.cpp \ ./modules/ModelProcessorx/HydrologyDCEfficient/CreateParametersHydrologyDCEfficient.cpp \ - ./solutions/hydrology_core.cpp\ - ./solvers/solver_hydro_nonlinear.cpp + ./analyses/hydrology_core.cpp\ + ./solvers/solver_hydro_nonlinear.cpp #}}} #Diagnostic sources {{{ @@ -478,7 +478,7 @@ ./modules/ModelProcessorx/DiagnosticHutter/CreateNodesDiagnosticHutter.cpp \ ./modules/ModelProcessorx/DiagnosticHutter/CreateConstraintsDiagnosticHutter.cpp \ - ./modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp \ - ./solutions/diagnostic_core.cpp\ - ./solvers/solver_stokescoupling_nonlinear.cpp + ./modules/ModelProcessorx/DiagnosticHutter/CreateLoadsDiagnosticHutter.cpp \ + ./analyses/diagnostic_core.cpp\ + ./solvers/solver_stokescoupling_nonlinear.cpp #}}} #Balanced sources {{{ @@ -486,7 +486,7 @@ ./modules/ModelProcessorx/Balancethickness/CreateNodesBalancethickness.cpp\ ./modules/ModelProcessorx/Balancethickness/CreateConstraintsBalancethickness.cpp\ - ./modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp\ - ./solutions/balancethickness_core.cpp \ - ./solutions/dummy_core.cpp + ./modules/ModelProcessorx/Balancethickness/CreateLoadsBalancethickness.cpp\ + ./analyses/balancethickness_core.cpp \ + ./analyses/dummy_core.cpp #}}} #Slope sources {{{ @@ -499,23 +499,23 @@ ./modules/ModelProcessorx/SurfaceSlope/CreateConstraintsSurfaceSlope.cpp\ ./modules/ModelProcessorx/SurfaceSlope/CreateLoadsSurfaceSlope.cpp\ - ./solutions/surfaceslope_core.cpp\ - ./solutions/bedslope_core.cpp + ./analyses/surfaceslope_core.cpp\ + ./analyses/bedslope_core.cpp #}}} #Gia sources {{{ -gia_sources = ./solutions/gia_core.cpp\ - ./modules/ModelProcessorx/Gia/UpdateElementsGia.cpp\ - ./modules/ModelProcessorx/Gia/CreateNodesGia.cpp \ - ./modules/ModelProcessorx/Gia/CreateConstraintsGia.cpp\ - ./modules/ModelProcessorx/Gia/CreateLoadsGia.cpp\ - ./modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp\ - ./modules/GiaDeflectionCorex/GiaDeflectionCorex.h\ - ./modules/GiaDeflectionCorex/distme.f\ - ./modules/GiaDeflectionCorex/freed.f\ - ./modules/GiaDeflectionCorex/ojrule.f\ - ./modules/GiaDeflectionCorex/pwise.f\ - ./modules/GiaDeflectionCorex/qwise.f\ - ./modules/GiaDeflectionCorex/stot.f\ - ./modules/GiaDeflectionCorex/what0.f\ - ./classes/GiaDeflectionCoreArgs.h +gia_sources = ./analyses/gia_core.cpp\ + ./modules/ModelProcessorx/Gia/UpdateElementsGia.cpp\ + ./modules/ModelProcessorx/Gia/CreateNodesGia.cpp \ + ./modules/ModelProcessorx/Gia/CreateConstraintsGia.cpp\ + ./modules/ModelProcessorx/Gia/CreateLoadsGia.cpp\ + ./modules/GiaDeflectionCorex/GiaDeflectionCorex.cpp\ + ./modules/GiaDeflectionCorex/GiaDeflectionCorex.h\ + ./modules/GiaDeflectionCorex/distme.f\ + ./modules/GiaDeflectionCorex/freed.f\ + ./modules/GiaDeflectionCorex/ojrule.f\ + ./modules/GiaDeflectionCorex/pwise.f\ + ./modules/GiaDeflectionCorex/qwise.f\ + ./modules/GiaDeflectionCorex/stot.f\ + ./modules/GiaDeflectionCorex/what0.f\ + ./classes/GiaDeflectionCoreArgs.h #}}} Index: /issm/trunk-jpl/src/c/analyses/AdjointCorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AdjointCorePointerFromSolutionEnum.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/AdjointCorePointerFromSolutionEnum.cpp (revision 14985) @@ -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 "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void AdjointCorePointerFromSolutionEnum(void (**padjointcore)(FemModel*),int solutiontype){ + + /*output: */ + void (*adjointcore)(FemModel*)=NULL; + + switch(solutiontype){ + + case DiagnosticSolutionEnum: + adjointcore=&adjointdiagnostic_core; + break; + case SteadystateSolutionEnum: + adjointcore=&adjointdiagnostic_core; + break; + case BalancethicknessSolutionEnum: + adjointcore=&adjointbalancethickness_core; + break; + case WeakBalancethicknessSolutionEnum: + 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/analyses/AnalysisConfiguration.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/AnalysisConfiguration.cpp (revision 14985) @@ -0,0 +1,132 @@ +/*!\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 "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void AnalysisConfiguration(int** panalyses,int* pnumanalyses, int solutiontype){ + + /*output: */ + int numanalyses; + int* analyses=NULL; + + /*Analyses lists*/ + switch(solutiontype){ + + case DiagnosticSolutionEnum: + numanalyses=5; + analyses=xNew(numanalyses); + analyses[0]=DiagnosticHorizAnalysisEnum; + analyses[1]=DiagnosticVertAnalysisEnum; + analyses[2]=DiagnosticHutterAnalysisEnum; + analyses[3]=SurfaceSlopeAnalysisEnum; + analyses[4]=BedSlopeAnalysisEnum; + break; + + case SteadystateSolutionEnum: + numanalyses=8; + analyses=xNew(numanalyses); + analyses[0]=DiagnosticHorizAnalysisEnum; + analyses[1]=DiagnosticVertAnalysisEnum; + analyses[2]=DiagnosticHutterAnalysisEnum; + analyses[3]=SurfaceSlopeAnalysisEnum; + analyses[4]=BedSlopeAnalysisEnum; + analyses[5]=EnthalpyAnalysisEnum; + analyses[6]=ThermalAnalysisEnum; + analyses[7]=MeltingAnalysisEnum; + break; + + case ThermalSolutionEnum: + numanalyses=2; + analyses=xNew(numanalyses); + analyses[0]=ThermalAnalysisEnum; + analyses[1]=MeltingAnalysisEnum; + break; + + case EnthalpySolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=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 PrognosticSolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=PrognosticAnalysisEnum; + break; + + case BalancethicknessSolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=BalancethicknessAnalysisEnum; + break; + + case WeakBalancethicknessSolutionEnum: + numanalyses=1; + analyses=xNew(numanalyses); + analyses[0]=BalancethicknessAnalysisEnum; + 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 TransientSolutionEnum: + numanalyses=10-1; + analyses=xNew(numanalyses); + analyses[0]=DiagnosticHorizAnalysisEnum; + analyses[1]=DiagnosticVertAnalysisEnum; + analyses[2]=DiagnosticHutterAnalysisEnum; + analyses[3]=SurfaceSlopeAnalysisEnum; + analyses[4]=BedSlopeAnalysisEnum; + analyses[5]=ThermalAnalysisEnum; + analyses[6]=MeltingAnalysisEnum; + analyses[7]=EnthalpyAnalysisEnum; + analyses[8]=PrognosticAnalysisEnum; + 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/analyses/CMakeLists.txt =================================================================== --- /issm/trunk-jpl/src/c/analyses/CMakeLists.txt (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/CMakeLists.txt (revision 14985) @@ -0,0 +1,61 @@ +# Subdirectories {{{ +# }}} +# Include Directory {{{ +include_directories(AFTER $ENV{ISSM_DIR}/src/c/solutions) +# }}} +# CORE_SOURCES {{{ +set(CORE_SOURCES $ENV{ISSM_DIR}/src/c/solutions/ad_core.cpp + $ENV{ISSM_DIR}/src/c/solutions/AnalysisConfiguration.cpp + $ENV{ISSM_DIR}/src/c/solutions/convergence.cpp +$ENV{ISSM_DIR}/src/c/solutions/CorePointerFromSolutionEnum.cpp + $ENV{ISSM_DIR}/src/c/solutions/EnvironmentFinalize.cpp + $ENV{ISSM_DIR}/src/c/solutions/EnvironmentInit.cpp + $ENV{ISSM_DIR}/src/c/solutions/ProcessArguments.cpp + $ENV{ISSM_DIR}/src/c/solutions/ResetBoundaryConditions.cpp +$ENV{ISSM_DIR}/src/c/solutions/WrapperCorePointerFromSolutionEnum.cpp PARENT_SCOPE) +# }}} +# DAKOTA_SOURCES {{{ +set(DAKOTA_SOURCES $ENV{ISSM_DIR}/src/c/solutions/dakota_core.cpp + $ENV{ISSM_DIR}/src/c/solutions/DakotaSpawnCore.cpp PARENT_SCOPE) +# }}} +# TRANSIENT_SOURCES {{{ +set(TRANSIENT_SOURCES $ENV{ISSM_DIR}/src/c/solutions/transient_core.cpp PARENT_SCOPE) +# }}} +# STEADYSTATE_SOURCES {{{ +set(STEADYSTATE_SOURCES $ENV{ISSM_DIR}/src/c/solutions/steadystate_core.cpp + $ENV{ISSM_DIR}/src/c/solutions/steadystateconvergence.cpp PARENT_SCOPE) +# }}} +# PROGNOSTIC_SOURCES {{{ +set(PROGNOSTIC_SOURCES $ENV{ISSM_DIR}/src/c/solutions/prognostic_core.cpp PARENT_SCOPE) +# }}} +# THERMAL_SOURCES {{{ +set(THERMAL_SOURCES $ENV{ISSM_DIR}/src/c/solutions/enthalpy_core.cpp + $ENV{ISSM_DIR}/src/c/solutions/thermal_core.cpp PARENT_SCOPE) +# }}} +# HYDROLOGY_SOURCES {{{ +set(HYDROLOGY_SOURCES $ENV{ISSM_DIR}/src/c/solutions/hydrology_core.cpp + $ENV{ISSM_DIR}/src/c/solutions/hydrology_core_step.cpp PARENT_SCOPE) +# }}} +# DIAGNOSTIC_SOURCES {{{ +set(DIAGNOSTIC_SOURCES $ENV{ISSM_DIR}/src/c/solutions/diagnostic_core.cpp PARENT_SCOPE) +# }}} +# BALANCED_SOURCES {{{ +set(BALANCED_SOURCES $ENV{ISSM_DIR}/src/c/solutions/balancethickness_core.cpp + $ENV{ISSM_DIR}/src/c/solutions/dummy_core.cpp PARENT_SCOPE) +# }}} +# SLOPE_SOURCES {{{ +set(SLOPE_SOURCES $ENV{ISSM_DIR}/src/c/solutions/bedslope_core.cpp + $ENV{ISSM_DIR}/src/c/solutions/surfaceslope_core.cpp PARENT_SCOPE) +# }}} +# LIBISSM_LA_SOURCES {{{ +set(LIBISSM_LA_SOURCES $ENV{ISSM_DIR}/src/c/solutions/issm.cpp PARENT_SCOPE) +# }}} +# ISSM_SOURCES {{{ +set(ISSM_SOURCES $ENV{ISSM_DIR}/src/c/solutions/issm.cpp PARENT_SCOPE) +# }}} +# KRIGING_SOURCES {{{ +set(KRIGING_SOURCES $ENV{ISSM_DIR}/src/c/solutions/kriging.cpp PARENT_SCOPE) +# }}} +# ISSMROSE_EXE_SOURCES {{{ +set(ISSMROSE_EXE_SOURCES $ENV{ISSM_DIR}/src/c/main/issm.cpp PARENT_SCOPE) +# }}} Index: /issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/CorePointerFromSolutionEnum.cpp (revision 14985) @@ -0,0 +1,120 @@ +/*!\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 "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype){ + + /*output: */ + void (*solutioncore)(FemModel*)=NULL; + + switch(solutiontype){ + + case DiagnosticSolutionEnum: + #ifdef _HAVE_DIAGNOSTIC_ + solutioncore=&diagnostic_core; + #else + _error_("ISSM was not compiled with diagnostic 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 EnthalpySolutionEnum: + #ifdef _HAVE_THERMAL_ + solutioncore=&enthalpy_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 WeakBalancethicknessSolutionEnum: + #ifdef _HAVE_BALANCED_ + solutioncore=&dummy_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 PrognosticSolutionEnum: + #ifdef _HAVE_PROGNOSTIC_ + solutioncore=&prognostic_core; + #else + _error_("ISSM was not compiled with prognostic capabilities. Exiting"); + #endif + break; + + case GiaSolutionEnum: + #ifdef _HAVE_GIA_ + solutioncore=&gia_core; + #else + _error_("ISSM was not compiled with gia 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/analyses/DakotaSpawnCore.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/DakotaSpawnCore.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/DakotaSpawnCore.cpp (revision 14985) @@ -0,0 +1,191 @@ +/*!\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 "./analyses.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.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; + char *string = NULL; + 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: */ + #ifdef _HAVE_MPI_ + MPI_Bcast(&counter,1,MPI_INT,0,IssmComm::GetComm()); + #endif + 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()) _pprintLine_("qmu iteration: " << counter); + + /* 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->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,d_variables,d_variables_descriptors,d_numvariables); + + /*Determine solution sequence: */ + if(VerboseQmu()) _pprintLine_("Starting " << EnumToStringx(solution_type) << " core:"); + WrapperCorePointerFromSolutionEnum(&solutioncore,femmodel->parameters,solution_type,nodakotacore); + + /*Run the core solution sequence: */ + solutioncore(femmodel); + + /*compute responses: */ + if(VerboseQmu()) _pprintLine_("compute dakota responses:"); + 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: */ + MPI_Bcast(&numvariables,1,MPI_INT,0,IssmComm::GetComm()); + + /*variables:*/ + if(my_rank!=0)variables=xNew(numvariables); + 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); + MPI_Bcast(string,string_length,MPI_CHAR,0,IssmComm::GetComm()); + if(my_rank!=0)variables_descriptors[i]=string; + } + + /*numresponses: */ + MPI_Bcast(&numresponses,1,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/analyses/DakotaSpawnCore.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/DakotaSpawnCore.h (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/DakotaSpawnCore.h (revision 14985) @@ -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/analyses/EnvironmentFinalize.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnvironmentFinalize.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/EnvironmentFinalize.cpp (revision 14985) @@ -0,0 +1,28 @@ +/*!\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){ + + #ifdef _HAVE_MPI_ + + /*Make sure we are all here*/ + MPI_Barrier(MPI_COMM_WORLD); + + /*Print closing statement*/ + int my_rank; + MPI_Comm_rank(MPI_COMM_WORLD,&my_rank); + if(!my_rank) printf("closing MPI\n"); + + /*Finalize: */ + MPI_Finalize(); + + #endif +} Index: /issm/trunk-jpl/src/c/analyses/EnvironmentInit.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnvironmentInit.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/EnvironmentInit.cpp (revision 14985) @@ -0,0 +1,37 @@ +/*!\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" + +COMM EnvironmentInit(int argc,char** argv){ + + /*Output*/ + COMM comm = 0; + + /*Initialize MPI environment: */ + #if defined(_HAVE_MPI_) + MPI_Init(&argc,&argv); + comm = MPI_COMM_WORLD; + #else + comm = 1; //bogus number for comm, which does not exist anyway. + #endif + + /*Print Banner*/ + int my_rank = 0; + #ifdef _HAVE_MPI_ + MPI_Comm_rank(comm,&my_rank); + #endif + if(!my_rank) printf("\n"); + if(!my_rank) printf("Ice Sheet System Model (%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/analyses/ProcessArguments.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ProcessArguments.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/ProcessArguments.cpp (revision 14985) @@ -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/analyses/ResetBoundaryConditions.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ResetBoundaryConditions.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/ResetBoundaryConditions.cpp (revision 14985) @@ -0,0 +1,31 @@ +/*!\file: ResetBoundaryConditions.cpp + * \brief: change boundary conditions of a model, using a solution vector from another analysis + */ + +#include "../classes/objects/objects.h" +#include "../modules/modules.h" +#include "../shared/io/io.h" + +void ResetBoundaryConditions(FemModel* femmodel, int analysis_type){ + + /*variables: */ + Vector* yg = NULL; + Nodes *nodes = NULL; + + if(VerboseSolution()) _pprintLine_(" updating boundary conditions..."); + + /*set current analysis: */ + femmodel->SetCurrentConfiguration(analysis_type); + + /*recover nodes: */ + nodes=femmodel->nodes; + + /*retrieve boundary conditions from element inputs :*/ + GetSolutionFromInputsx( &yg, femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials, femmodel->parameters); + + /*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/analyses/WrapperCorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/WrapperCorePointerFromSolutionEnum.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/WrapperCorePointerFromSolutionEnum.cpp (revision 14985) @@ -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 "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.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/analyses/ad_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/ad_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/ad_core.cpp (revision 14985) @@ -0,0 +1,307 @@ +/*!\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 "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/shared.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" +/*}}}*/ + +void ad_core(FemModel* femmodel){ + + /*diverse: */ + int i; + int dummy; + int num_dependents; + int num_independents; + bool isautodiff = false; + char *driver = NULL; + size_t tape_stats[11]; + + /*state variables: */ + IssmDouble *axp = NULL; + double *xp = NULL; + + /*AD mode on?: */ + femmodel->parameters->FindParam(&isautodiff,AutodiffIsautodiffEnum); + + if(isautodiff){ + + #ifdef _HAVE_ADOLC_ + + /*First, stop tracing: */ + trace_off(); + + /*preliminary checks: */ + femmodel->parameters->FindParam(&num_dependents,AutodiffNumDependentsEnum); + femmodel->parameters->FindParam(&num_independents,AutodiffNumIndependentsEnum); + if(!(num_dependents*num_independents)) return; + + if(VerboseAutodiff())_pprintLine_(" start ad core"); + + /*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: */ + anEDF_for_solverx_p->fos_forward=EDF_fos_forward_for_solverx; + + /*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(1,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: */ + anEDF_for_solverx_p->fov_forward=EDF_fov_forward_for_solverx; + // 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(1,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); + + if (aDepIndex<0 || aDepIndex>=num_dependents) _error_("index value for AutodiffFosReverseIndexEnum should be in [0,num_dependents-1]"); + + aWeightVector=xNewZeroInit(num_dependents); + aWeightVector[aDepIndex]=1.0; + + weightVectorTimesJac=xNew(num_independents); + + /*set the forward method function pointer: */ + anEDF_for_solverx_p->fos_reverse=EDF_fos_reverse_for_solverx; + + /*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(1,num_dependents,num_independents, aWeightVector, weightVectorTimesJac ); + + /*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: */ + anEDF_for_solverx_p->fov_reverse=EDF_fov_reverse_for_solverx; + + /*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(1,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); + + /*Print statistics:*/ + tapestats(1,tape_stats); //reading of tape statistics + if(VerboseAutodiff()){ + _pprintLine_(" ADOLC statistics: "); + _pprintLine_(" "<parameters->FindParam(&save_results,SaveResultsEnum); + + /*compute thickness */ + if(VerboseSolution()) _pprintLine_(" computing thickness"); + femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); + solver_linear(femmodel); + + /*Call SurfaceAreax, because some it might be needed by PVector*/ + SurfaceAreax(NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + /*compute adjoint*/ + if(VerboseSolution()) _pprintLine_(" computing adjoint"); + femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum,AdjointBalancethicknessAnalysisEnum); + solver_adjoint_linear(femmodel); + + /*Save results*/ + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointEnum); + } +} Index: /issm/trunk-jpl/src/c/analyses/adjointdiagnostic_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/adjointdiagnostic_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/adjointdiagnostic_core.cpp (revision 14985) @@ -0,0 +1,47 @@ +/*!\file: adjointdiagnostic_core.cpp + * \brief compute inverse method adjoint state + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void adjointdiagnostic_core(FemModel* femmodel){ + + /*parameters: */ + bool isstokes; + bool save_results; + bool conserve_loads = true; + + /*retrieve parameters:*/ + femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum); + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + /*Compute velocities*/ + if(VerboseSolution()) _pprintLine_(" computing velocities"); + femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); + solver_nonlinear(femmodel,conserve_loads); + + /*Call SurfaceAreax, because some it might be needed by PVector*/ + SurfaceAreax(NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + /*Compute adjoint*/ + if(VerboseSolution()) _pprintLine_(" computing adjoint"); + femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum,AdjointHorizAnalysisEnum); + solver_adjoint_linear(femmodel); + + /*Save results*/ + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointyEnum); + if (isstokes){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointzEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,AdjointpEnum); + } + } +} Index: /issm/trunk-jpl/src/c/analyses/analyses.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/analyses.h (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/analyses.h (revision 14985) @@ -0,0 +1,64 @@ +/* + * analyses.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/Comm.h" +#include "../shared/Numerics/types.h" + +/*cores: */ +void adjointdiagnostic_core(FemModel* femmodel); +void adjointbalancethickness_core(FemModel* femmodel); +void gradient_core(FemModel* femmodel,int n=0,bool orthogonalize=false); +void diagnostic_core(FemModel* femmodel); +void hydrology_core(FemModel* femmodel); +void thermal_core(FemModel* femmodel); +void enthalpy_core(FemModel* femmodel); +void surfaceslope_core(FemModel* femmodel); +void bedslope_core(FemModel* femmodel); +void control_core(FemModel* femmodel); +void controltao_core(FemModel* femmodel); +void prognostic_core(FemModel* femmodel); +void balancethickness_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); +IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs); + +//convergence: +void convergence(bool* pconverged, Matrix* K_ff,Vector* p_f,Vector* u_f,Vector* u_f_old,Parameters* parameters); +bool controlconvergence(IssmDouble J,IssmDouble tol_cm); +bool steadystateconvergence(FemModel* femmodel); + +//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 controlrestart(FemModel* femmodel,IssmDouble* J); +void ResetBoundaryConditions(FemModel* femmodel, int analysis_type); +COMM EnvironmentInit(int argc,char** argv); +void EnvironmentFinalize(void); +void PrintBanner(void); + +//solution configuration +void AnalysisConfiguration(int** panalyses,int* pnumanalyses, 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/analyses/balancethickness_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/balancethickness_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/balancethickness_core.cpp (revision 14985) @@ -0,0 +1,32 @@ +/*!\file: balancethickness_core.cpp + * \brief: core of the balancethickness solution + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.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()) _pprintLine_("call computational core:"); + solver_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + } + +} Index: /issm/trunk-jpl/src/c/analyses/bedslope_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/bedslope_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/bedslope_core.cpp (revision 14985) @@ -0,0 +1,35 @@ +/*!\file: bedslope_core.cpp + * \brief: core of the slope solution + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void bedslope_core(FemModel* femmodel){ + + /*parameters: */ + bool save_results; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + if(VerboseSolution()) _pprintLine_(" computing slope"); + + /*Call on core computations: */ + femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeXAnalysisEnum); + solver_linear(femmodel); + femmodel->SetCurrentConfiguration(BedSlopeAnalysisEnum,BedSlopeYAnalysisEnum); + solver_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeXEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedSlopeYEnum); + } + +} Index: /issm/trunk-jpl/src/c/analyses/control_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/control_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/control_core.cpp (revision 14985) @@ -0,0 +1,134 @@ +/*!\file: control_core.cpp + * \brief: core of the control solution + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void control_core(FemModel* femmodel){ + + int i,n; + + /*parameters: */ + int num_controls,num_responses; + int nsteps; + IssmDouble tol_cm; + bool cm_gradient; + int dim; + int solution_type; + bool isstokes; + bool dakota_analysis=false; + + int* control_type = NULL; + IssmDouble* responses=NULL; + int* step_responses=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(&num_responses,InversionNumCostFunctionsEnum); + femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + femmodel->parameters->FindParam(&responses,NULL,NULL,InversionCostFunctionsEnum); + 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(&cm_gradient,InversionGradientOnlyEnum); + femmodel->parameters->FindParam(&dim,MeshDimensionEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum); + 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()) _pprintLine_(" preparing initial solution"); + if(isstokes) solutioncore(femmodel); + + /*Initialize responses: */ + J=xNew(nsteps); + step_responses=xNew(num_responses); + + /*Initialize some of the BrentSearch arguments: */ + optargs.femmodel=femmodel; + optpars.xmin=0; optpars.xmax=1; + + /*Start looping: */ + for(n=0;n(responses[n*num_responses+i]); + femmodel->parameters->SetParam(step_responses,1,num_responses,StepResponsesEnum); + + /*In steady state inversion, compute new temperature field now*/ + if(solution_type==SteadystateSolutionEnum) solutioncore(femmodel); + + if(VerboseControl()) _pprintLine_(" compute adjoint state:"); + adjointcore(femmodel); + gradient_core(femmodel,n,search_scalar==0); + + /*Return gradient if asked: */ + if (cm_gradient){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GradientEnum); + goto cleanup_and_return; + } + + if(VerboseControl()) _pprintLine_(" optimizing along gradient direction"); + optpars.maxiter=reCast(maxiter[n]); optpars.cm_jump=cm_jump[n]; + BrentSearch(&search_scalar,J+n,&optpars,&objectivefunction,&optargs); + + if(VerboseControl()) _pprintLine_(" updating parameter using optimized search scalar"); //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()) _pprintLine_(" preparing final solution"); + femmodel->parameters->SetParam(true,SaveResultsEnum); + solutioncore(femmodel); + + /*some results not computed by steadystate_core or diagnostic_core: */ + if(!dakota_analysis){ //do not save this if we are running the control core from a qmu run! + for(i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]); + + #ifdef _HAVE_ADOLC_ + IssmPDouble* J_passive=xNew(nsteps); + for(int 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 + } + + cleanup_and_return: + /*Free ressources: */ + xDelete(control_type); + xDelete(step_responses); + xDelete(maxiter); + xDelete(responses); + xDelete(cm_jump); + xDelete(J); +} Index: /issm/trunk-jpl/src/c/analyses/controlconvergence.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/controlconvergence.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/controlconvergence.cpp (revision 14985) @@ -0,0 +1,28 @@ +/*!\file: controlconvergence.cpp + * \brief: determine convergence of control_core solution + */ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./analyses.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../Container/Container.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" + +bool controlconvergence(IssmDouble J, IssmDouble tol_cm){ + + bool converged=false; + + /*Has convergence been reached?*/ + if (!xIsNan(tol_cm) && Jparameters->FindParam(&num_controls,InversionNumControlParametersEnum); + femmodel->parameters->FindParam(&control_type,NULL,InversionControlParametersEnum); + femmodel->parameters->FindParam(&nsteps,InversionNstepsEnum); + femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + + /*only save if we are not running qmu analysis. We certainly don't want to save control results each time we + * run on control core!: */ + if(!dakota_analysis){ + /*we essentially want J and the parameter: */ + for(int i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_type[i]); + #ifdef _HAVE_ADOLC_ + IssmPDouble* J_passive=xNew(nsteps); + for(int 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 + //femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,InversionControlParametersEnum,EnumToStringx(control_type),1,0)); + + /*write to disk: */ + OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results); + } + + /*Clean up and return*/ + xDelete(control_type); +} Index: /issm/trunk-jpl/src/c/analyses/controltao_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/controltao_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/controltao_core.cpp (revision 14985) @@ -0,0 +1,190 @@ +/*!\file: control_core.cpp + * \brief: core of the control solution + */ +#include +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.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; + IssmDouble *dummy = NULL; + int *control_list = NULL; + Vector *X = 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()) _pprintLine_(" Initializing the Toolkit for Advanced Optimization (TAO)"); + 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()) _pprintLine_(" Starting optimization"); + TaoSolve(tao); + TaoView(tao,PETSC_VIEWER_STDOUT_WORLD); + TaoGetSolutionVector(tao,&X->pvector->vector); + SetControlInputsFromVectorx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,X); + for(int i=0;ielements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,control_list[i]); + } + femmodel->results->AddObject(new GenericExternalResult(femmodel->results->Size()+1,JEnum,user.J,maxiter+3,1,1,0)); + + /*Finalize*/ + if(VerboseControl()) _pprintLine_(" preparing final solution"); + 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,num_cost_functions; + AppCtx *user = (AppCtx *)userCtx; + FemModel *femmodel = user->femmodel; + int *cost_functions = NULL; + IssmDouble *cost_functionsd = NULL; + 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); + femmodel->parameters->FindParam(&num_cost_functions,InversionNumCostFunctionsEnum); + femmodel->parameters->FindParam(&cost_functionsd,NULL,NULL,InversionCostFunctionsEnum); + + /*Prepare objective function*/ + cost_functions=xNew(num_cost_functions); + for(int i=0;iparameters->SetParam(cost_functions,1,num_cost_functions,StepResponsesEnum); + + /*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*/ + xDelete(cost_functions); + xDelete(cost_functionsd); + return 0; +} +int IssmMonitor(TaoSolver tao, void *userCtx){ + + int i,its,num_responses; + IssmDouble f,gnorm,cnorm,xdiff; + AppCtx *user = (AppCtx *)userCtx; + FemModel *femmodel = user->femmodel; + Element *element = NULL; + int *responses = NULL; + + femmodel->parameters->FindParam(&num_responses,InversionNumCostFunctionsEnum); + femmodel->parameters->FindParam(&responses,NULL,NULL,StepResponsesEnum); + + TaoGetSolutionStatus(tao, &its, &f, &gnorm, &cnorm, &xdiff, NULL); + if(its==0) _pprintLine_("Iter Function Residual | List of contributions"); + if(its==0) _pprintLine_("-----------------------------------+-----------------------"); + _pprintString_(setw(4)<J[its]=f; + + /*Retrieve objective functions independently*/ + for(i=0;iResponsex(&f,EnumToStringx(responses[i]),false,i); + _pprintString_(" "<(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/analyses/convergence.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/convergence.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/convergence.cpp (revision 14985) @@ -0,0 +1,147 @@ +/*!\file: convergence.cpp + * \brief: figure out if convergence has been reached + */ + +#include "../classes/objects/objects.h" +#include "../modules/modules.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" + +void convergence(bool* pconverged, Matrix* Kff,Vector* pf,Vector* uf,Vector* old_uf,Parameters* parameters){ + + /*output*/ + bool converged=false; + + /*intermediary*/ + Vector* KU=NULL; + Vector* KUF=NULL; + Vector* KUold=NULL; + Vector* KUoldF=NULL; + Vector* duf=NULL; + IssmDouble ndu,nduinf,nu; + IssmDouble nKUF; + IssmDouble nKUoldF; + IssmDouble nF; + IssmDouble solver_residue,res; + + /*convergence options*/ + IssmDouble eps_res; + IssmDouble eps_rel; + IssmDouble eps_abs; + IssmDouble yts; + + if(VerboseModule()) _pprintLine_(" checking convergence"); + + /*If uf is NULL in input, f-set is nil, model is fully constrained, therefore converged from + * the get go: */ + if(uf->IsEmpty()){ + *pconverged=true; + return; + } + + /*get convergence options*/ + parameters->FindParam(&eps_res,DiagnosticRestolEnum); + parameters->FindParam(&eps_rel,DiagnosticReltolEnum); + parameters->FindParam(&eps_abs,DiagnosticAbstolEnum); + parameters->FindParam(&yts,ConstantsYtsEnum); + + /*Display solver caracteristics*/ + if (VerboseConvergence()){ + + //compute KUF = KU - F = K*U - F + KU=uf->Duplicate(); Kff->MatMult(uf,KU); + KUF=KU->Duplicate(); KU->Copy(KUF); KUF->AYPX(pf,-1.0); + + //compute norm(KUF), norm(F) and residue + nKUF=KUF->Norm(NORM_TWO); + nF=pf->Norm(NORM_TWO); + solver_residue=nKUF/nF; + _pprintLine_("\n" << " solver residue: norm(KU-F)/norm(F)=" << solver_residue); + + //clean up + delete KU; + delete KUF; + } + + /*Force equilibrium (Mandatory)*/ + + //compute K[n]U[n-1]F = K[n]U[n-1] - F + _assert_(uf); _assert_(Kff); + KUold=uf->Duplicate(); Kff->MatMult(old_uf,KUold); + KUoldF=KUold->Duplicate();KUold->Copy(KUoldF); KUoldF->AYPX(pf,-1.0); + nKUoldF=KUoldF->Norm(NORM_TWO); + nF=pf->Norm(NORM_TWO); + res=nKUoldF/nF; + if (xIsNan(res)){ + _pprintLine_("norm nf = " << nF << "f and norm kuold = " << nKUoldF << "f"); + _error_("mechanical equilibrium convergence criterion is NaN!"); + } + + //clean up + delete KUold; + delete KUoldF; + + //print + if(res "<(eps_rel) || (VerboseConvergence())){ + + //compute norm(du)/norm(u) + duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0); + ndu=duf->Norm(NORM_TWO); nu=old_uf->Norm(NORM_TWO); + + if (xIsNan(ndu) || xIsNan(nu)) _error_("convergence criterion is NaN!"); + + //clean up + delete duf; + + //print + if (!xIsNan(eps_rel)){ + if((ndu/nu) " << eps_rel*100 << " %"); + converged=false; + } + } + else _pprintLine_(setw(50) << left << " Convergence criterion: norm(du)/norm(u)" << ndu/nu*100 << " %"); + + } + + /*Absolute criterion (Optional) = max(du)*/ + if (!xIsNan(eps_abs) || (VerboseConvergence())){ + + //compute max(du) + duf=old_uf->Duplicate(); old_uf->Copy(duf); duf->AYPX(uf,-1.0); + ndu=duf->Norm(NORM_TWO); nduinf=duf->Norm(NORM_INF); + if (xIsNan(ndu) || xIsNan(nu)) _error_("convergence criterion is NaN!"); + + //clean up + delete duf; + + //print + if (!xIsNan(eps_abs)){ + if ((nduinf*yts) " << eps_abs << " m/yr"); + converged=false; + } + } + else _pprintLine_(setw(50) << left << " Convergence criterion: max(du)" << nduinf*yts << " m/yr"); + + } + + /*assign output*/ + *pconverged=converged; +} Index: /issm/trunk-jpl/src/c/analyses/dakota_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/dakota_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/dakota_core.cpp (revision 14985) @@ -0,0 +1,128 @@ +/*!\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 (diagnostic_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 "./analyses.h" +#include "../shared/io/io.h" +#include "../toolkits/toolkits.h" +#include "../shared/Enum/Enum.h" +#include "../classes/dakota/DakotaPlugin.h" +#include "../classes/FemModel.h" +#include "../Container/Parameters.h" + +#ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. +#include +#include +#include +#include +#include +#include "./DakotaSpawnCore.h" +#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; + int status = 0; + Dakota::ModelLIter ml_iter; + Parameters *parameters = NULL; + + /*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. + Dakota::ParallelLibrary parallel_lib("serial"); //use our own ISSM Dakota library mode constructor, which only fires up Dakota on CPU 0. + Dakota::ProblemDescDB problem_db(parallel_lib); + + // 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/analyses/diagnostic_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/diagnostic_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/diagnostic_core.cpp (revision 14985) @@ -0,0 +1,109 @@ +/*!\file: diagnostic_core.cpp + * \brief: core of the diagnostic solution + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void diagnostic_core(FemModel* femmodel){ + + /*parameters: */ + bool dakota_analysis = false; + int dim = -1; + bool ishutter = false; + bool ismacayealpattyn = false; + bool isl1l2 = false; + bool isstokes = false; + bool conserve_loads = true; + bool modify_loads = true; + bool save_results; + int newton; + int solution_type; + int numoutputs = 0; + int *requested_outputs = NULL; + + /* recover parameters:*/ + femmodel->parameters->FindParam(&dim,MeshDimensionEnum); + femmodel->parameters->FindParam(&ishutter,FlowequationIshutterEnum); + femmodel->parameters->FindParam(&ismacayealpattyn,FlowequationIsmacayealpattynEnum); + femmodel->parameters->FindParam(&isl1l2,FlowequationIsl1l2Enum); + femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum); + femmodel->parameters->FindParam(&newton,DiagnosticIsnewtonEnum); + femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&numoutputs,DiagnosticNumRequestedOutputsEnum); + if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,DiagnosticRequestedOutputsEnum); + + /*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==DiagnosticSolutionEnum){ + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVxEnum,VxEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVyEnum,VyEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVzEnum,VzEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuPressureEnum,PressureEnum); + } + + /*Compute slopes: */ + if(ishutter) surfaceslope_core(femmodel); + if(isstokes){ + bedslope_core(femmodel); + femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); + ResetCoordinateSystemx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + + if(ishutter){ + + if(VerboseSolution()) _pprintLine_(" computing hutter velocities"); + + //Take the last velocity into account so that the velocity on the MacAyeal domain is not zero + if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticHutterAnalysisEnum); + + femmodel->SetCurrentConfiguration(DiagnosticHutterAnalysisEnum); + solver_linear(femmodel); + + if (ismacayealpattyn) ResetBoundaryConditions(femmodel,DiagnosticHorizAnalysisEnum); + } + + if ((ismacayealpattyn || isl1l2) ^ isstokes){ // ^ = xor + + if(VerboseSolution()) _pprintLine_(" computing velocities"); + femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); + if(newton>0) + solver_newton(femmodel); + else + solver_nonlinear(femmodel,conserve_loads); + } + + if (ismacayealpattyn && isstokes){ + + if(VerboseSolution()) _pprintLine_(" computing coupling macayealpattyn and stokes velocities and pressure "); + solver_stokescoupling_nonlinear(femmodel,conserve_loads); + } + + if (dim==3 & (ishutter || ismacayealpattyn)){ + + if(VerboseSolution()) _pprintLine_(" computing vertical velocities"); + femmodel->SetCurrentConfiguration(DiagnosticVertAnalysisEnum); + solver_linear(femmodel); + } + + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); + if(dim==3) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum); + femmodel->RequestedOutputsx(requested_outputs,numoutputs); + } + + if(solution_type==DiagnosticSolutionEnum)femmodel->RequestedDependentsx(); + + /*Free ressources:*/ + xDelete(requested_outputs); +} Index: /issm/trunk-jpl/src/c/analyses/dummy_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/dummy_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/dummy_core.cpp (revision 14985) @@ -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/analyses/enthalpy_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/enthalpy_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/enthalpy_core.cpp (revision 14985) @@ -0,0 +1,34 @@ +/*!\file: enthalpy_core.cpp + * \brief: core of the enthalpy solution + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../shared/io/io.h" +#include "../solvers/solvers.h" + +void enthalpy_core(FemModel* femmodel){ + + /*intermediary*/ + bool save_results; + + //first recover parameters common to all solutions + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + if(VerboseSolution()) _pprintLine_(" computing enthalpy"); + femmodel->SetCurrentConfiguration(EnthalpyAnalysisEnum); + solver_nonlinear(femmodel,true); + + /*transfer enthalpy to enthalpy picard for the next step: */ + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,EnthalpyPicardEnum); + + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WaterfractionEnum); + } +} Index: /issm/trunk-jpl/src/c/analyses/gia_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/gia_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/gia_core.cpp (revision 14985) @@ -0,0 +1,58 @@ +/*!\file: gia_core.cpp + * \brief: core of the GIA solution + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" +void gia_core(FemModel* femmodel){ + + int i; + 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()) _pprintLine_(" computing GIA"); + + /*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(); + + InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,wg); + InputUpdateFromVectorx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,dwdtg,GiadWdtEnum,GiaAnalysisEnum); + + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GiaWEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GiadWdtEnum); + } + +} Index: /issm/trunk-jpl/src/c/analyses/gradient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/gradient_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/gradient_core.cpp (revision 14985) @@ -0,0 +1,49 @@ +/*!\file: gradient_core.cpp + * \brief compute inverse method gradient direction. + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.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()) _pprintLine_(" compute cost function gradient"); + Gradjx(&gradient,&norm_list,femmodel->elements,femmodel->nodes, femmodel->vertices,femmodel->loads, femmodel->materials,femmodel->parameters); + + if (orthogonalize){ + if(VerboseControl()) _pprintLine_(" orthogonalization"); + 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/analyses/hydrology_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/hydrology_core.cpp (revision 14985) @@ -0,0 +1,98 @@ +/*!\file: hydrology_core.cpp + * \brief: core of the hydrology solution + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.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()) _pprintLine_(" computing water column"); + femmodel->SetCurrentConfiguration(HydrologyShreveAnalysisEnum); + solver_nonlinear(femmodel,modify_loads); + + /*transfer water column thickness to old water column thickness: */ + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WatercolumnEnum,WaterColumnOldEnum); + + if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){ + if(VerboseSolution()) _pprintLine_(" saving results "); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WatercolumnEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,HydrologyWaterVxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,HydrologyWaterVyEnum); + + /*unload results*/ + if(VerboseSolution()) _pprintLine_(" saving temporary results"); + OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results); + } + } + + else if (hydrology_model==HydrologydcEnum){ + femmodel->parameters->FindParam(&isefficientlayer,HydrologydcIsefficientlayerEnum); + if(VerboseSolution()) _pprintLine_(" computing water head"); + solver_hydro_nonlinear(femmodel); + if(save_results && ((i+1)%output_frequency==0 || (i+1)==nsteps)){ + if(VerboseSolution()) _pprintLine_(" saving results "); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SedimentHeadEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SedimentHeadResidualEnum); + if(isefficientlayer){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EplHeadEnum); + } + /*unload results*/ + if(VerboseSolution()) _pprintLine_(" saving temporary results"); + OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results); + } + } + + } +} Index: /issm/trunk-jpl/src/c/analyses/objectivefunction.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/objectivefunction.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/objectivefunction.cpp (revision 14985) @@ -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 "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/Enum/Enum.h" +#include "../solvers/solvers.h" +#include "../modules/modules.h" +/*}}}*/ + +IssmDouble objectivefunction(IssmDouble search_scalar,OptArgs* optargs){ + + /*output: */ + IssmDouble J; + + /*parameters: */ + int solution_type,analysis_type; + bool isstokes = false; + bool conserve_loads = true; + FemModel *femmodel = NULL; + + /*Recover finite element model: */ + femmodel=optargs->femmodel; + + /*Recover parameters: */ + femmodel->parameters->FindParam(&isstokes,FlowequationIsstokesEnum); + femmodel->parameters->FindParam(&analysis_type,AnalysisTypeEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + + /*set analysis type to compute velocity: */ + if (solution_type==SteadystateSolutionEnum || solution_type==DiagnosticSolutionEnum){ + femmodel->SetCurrentConfiguration(DiagnosticHorizAnalysisEnum); + } + else if (solution_type==BalancethicknessSolutionEnum){ + femmodel->SetCurrentConfiguration(BalancethicknessAnalysisEnum); + } + else if (solution_type==WeakBalancethicknessSolutionEnum){ + 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 diagnostic with updated inputs: */ + if (solution_type==SteadystateSolutionEnum){ + diagnostic_core(femmodel); //We need a 3D velocity!! (vz is required for the next thermal run) + } + else if (solution_type==DiagnosticSolutionEnum){ + solver_nonlinear(femmodel,conserve_loads); + } + else if (solution_type==BalancethicknessSolutionEnum){ + solver_linear(femmodel); + } + else if (solution_type==WeakBalancethicknessSolutionEnum){ + /*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/analyses/prognostic_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/prognostic_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/prognostic_core.cpp (revision 14985) @@ -0,0 +1,59 @@ +/*!\file: prognostic_core.cpp + * \brief: core of the prognostic solution + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void prognostic_core(FemModel* femmodel){ + + /*parameters: */ + bool save_results; + bool issmbgradients,ispdd,isdelta18o; + int solution_type; + int *requested_outputs = NULL; + int numoutputs=0; + + /*activate formulation: */ + femmodel->SetCurrentConfiguration(PrognosticAnalysisEnum); + + /*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(&solution_type,SolutionTypeEnum); + femmodel->parameters->FindParam(&numoutputs,PrognosticNumRequestedOutputsEnum); + if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,PrognosticRequestedOutputsEnum); + + if(issmbgradients){ + _printf_(VerboseSolution()," call smb gradients module\n"); + SmbGradientsx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + if(ispdd){ + if(isdelta18o){ + if(VerboseSolution()) _pprintLine_(" call Delta18oParametrization module"); + Delta18oParameterizationx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + if(VerboseSolution()) _pprintLine_(" call positive degree day module"); + PositiveDegreeDayx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + if(VerboseSolution()) _pprintLine_(" call computational core"); + solver_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + femmodel->RequestedOutputsx(requested_outputs,numoutputs); + } + + if(solution_type==PrognosticSolutionEnum)femmodel->RequestedDependentsx(); + + /*Free ressources:*/ + xDelete(requested_outputs); +} Index: /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/steadystate_core.cpp (revision 14985) @@ -0,0 +1,93 @@ +/*!\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 "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void steadystate_core(FemModel* femmodel){ + + /*intermediary: */ + int step; + + /*parameters: */ + bool save_results,isenthalpy; + int maxiter; + int numoutputs = 0; + int *requested_outputs = NULL; + + /* recover parameters:*/ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&maxiter,SteadystateMaxiterEnum); + femmodel->parameters->FindParam(&numoutputs,SteadystateNumRequestedOutputsEnum); + femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); + femmodel->parameters->SetParam(false,SaveResultsEnum); + if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SteadystateRequestedOutputsEnum); + + /*intialize counters: */ + step=1; + + for(;;){ + + if(VerboseSolution()) _pprintLine_(" computing temperature and velocity for step: " << step); + #ifdef _HAVE_THERMAL_ + if(isenthalpy==0){ + thermal_core(femmodel); + } + else{ + enthalpy_core(femmodel); + } + #else + _error_("ISSM was not compiled with thermal capabilities. Exiting"); + #endif + + if(VerboseSolution()) _pprintLine_(" computing new velocity"); + diagnostic_core(femmodel); + + if (step>1){ + if(VerboseSolution()) _pprintLine_(" checking velocity, temperature and pressure convergence"); + if(steadystateconvergence(femmodel)) break; + } + if(step>maxiter){ + if(VerboseSolution()) _pprintLine_(" maximum number steadystate iterations " << maxiter << " reached"); + break; + } + + if(VerboseSolution()) _pprintLine_(" saving velocity, temperature and pressure to check for convergence at next step"); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,VxPicardEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,VyPicardEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,VzPicardEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressurePicardEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureOldEnum); + + //increase counter + step++; + } + + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VelEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum); + if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WaterfractionEnum); + if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum); + if(!isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum); + femmodel->RequestedOutputsx(requested_outputs,numoutputs); + } + + /*Free ressources:*/ + xDelete(requested_outputs); +} Index: /issm/trunk-jpl/src/c/analyses/steadystateconvergence.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/steadystateconvergence.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/steadystateconvergence.cpp (revision 14985) @@ -0,0 +1,40 @@ +/*!\file: steadystateconvergence.cpp + * \brief: determine convergence of steady state solution + */ + +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif +#include "./analyses.h" +#include "../classes/objects/objects.h" +#include "../Container/Container.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" + +bool steadystateconvergence(FemModel* femmodel){ + + /*output: */ + bool converged=false; + bool velocity_converged=false; + bool temperature_converged=false; + + /*intermediary: */ + int velocityenums[8]={VxEnum,VxPicardEnum,VyEnum,VyPicardEnum,VzEnum,VzPicardEnum,PressureEnum,PressurePicardEnum}; //pairs of enums (new and old) on which to carry out the converence tests + int temperatureenums[2]={TemperatureEnum,TemperatureOldEnum}; + int convergencecriterion[1]={RelativeEnum}; //criterions for convergence, RelativeEnum or AbsoluteEnum + IssmDouble convergencecriterionvalue[1]; //value of criterion to be respected + + /*retrieve parameters: */ + femmodel->parameters->FindParam(&convergencecriterionvalue[0],SteadystateReltolEnum); + + /*figure out convergence at the input level, because we don't have the solution vectors!: */ + velocity_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,&velocityenums[0],8,&convergencecriterion[0],&convergencecriterionvalue[0],1); + temperature_converged=InputConvergencex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,&temperatureenums[0],2,&convergencecriterion[0],&convergencecriterionvalue[0],1); + + if(velocity_converged && temperature_converged) converged=true; + + /*return: */ + return converged; +} Index: /issm/trunk-jpl/src/c/analyses/surfaceslope_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/surfaceslope_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/surfaceslope_core.cpp (revision 14985) @@ -0,0 +1,35 @@ +/*!\file: surfaceslope_core.cpp + * \brief: core of the slope solution + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../solvers/solvers.h" +#include "../modules/modules.h" + +void surfaceslope_core(FemModel* femmodel){ + + /*parameters: */ + bool save_results; + + /*Recover some parameters: */ + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + + if(VerboseSolution()) _pprintLine_("computing slope..."); + + /*Call on core computations: */ + femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum); + solver_linear(femmodel); + femmodel->SetCurrentConfiguration(SurfaceSlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum); + solver_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _pprintLine_("saving results:"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeXEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceSlopeYEnum); + } + +} Index: /issm/trunk-jpl/src/c/analyses/thermal_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/thermal_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/thermal_core.cpp (revision 14985) @@ -0,0 +1,41 @@ +/*!\file: thermal_core.cpp + * \brief: core of the thermal solution + */ + +#include "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void thermal_core(FemModel* femmodel){ + + /*intermediary*/ + bool save_results; + bool dakota_analysis = false; + + //first recover parameters common to all solutions + femmodel->parameters->FindParam(&save_results,SaveResultsEnum); + femmodel->parameters->FindParam(&dakota_analysis,QmuIsdakotaEnum); + + if(dakota_analysis){ + femmodel->SetCurrentConfiguration(ThermalAnalysisEnum); + ResetConstraintsx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + } + + if(VerboseSolution()) _pprintLine_(" computing temperatures"); + femmodel->SetCurrentConfiguration(ThermalAnalysisEnum); + solver_thermal_nonlinear(femmodel); + + if(VerboseSolution()) _pprintLine_(" computing melting"); + femmodel->SetCurrentConfiguration(MeltingAnalysisEnum); + solver_linear(femmodel); + + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum); + } +} Index: /issm/trunk-jpl/src/c/analyses/transient_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 14985) +++ /issm/trunk-jpl/src/c/analyses/transient_core.cpp (revision 14985) @@ -0,0 +1,170 @@ +/*!\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 "./analyses.h" +#include "../toolkits/toolkits.h" +#include "../classes/objects/objects.h" +#include "../shared/io/io.h" +#include "../shared/Enum/Enum.h" +#include "../modules/modules.h" +#include "../solvers/solvers.h" + +void transient_core(FemModel* femmodel){ + + /*parameters: */ + IssmDouble starttime,finaltime,dt,yts; + bool isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy,isdelta18o,isgia; + bool save_results,dakota_analysis; + bool time_adapt=false; + int output_frequency; + int dim,groundingline_migration; + int numoutputs = 0; + int *requested_outputs = NULL; + + /*intermediary: */ + int step; + IssmDouble time; + + //first recover parameters common to all solutions + femmodel->parameters->FindParam(&dim,MeshDimensionEnum); + 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(&isdiagnostic,TransientIsdiagnosticEnum); + femmodel->parameters->FindParam(&isprognostic,TransientIsprognosticEnum); + femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum); + femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum); + femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum); + femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum); + 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){ + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVxEnum,VxEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVyEnum,VyEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuVzEnum,VzEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuPressureEnum,PressureEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuThicknessEnum,ThicknessEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuSurfaceEnum,SurfaceEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuBedEnum,BedEnum); + if(isthermal && dim==3){ + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuTemperatureEnum,TemperatureEnum); + InputDuplicatex(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,QmuMeltingEnum,BasalforcingsMeltingRateEnum); + } + } + + 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()) _pprintLine_("iteration " << step << "/" << floor((finaltime-time)/dt)+step << " time [yr]: " << time/yts << " (time step: " << dt/yts << ")"); + if(step%output_frequency==0 || time==finaltime) + save_results=true; + else + save_results=false; + femmodel->parameters->SetParam(save_results,SaveResultsEnum); + + if(isthermal && dim==3){ + if(VerboseSolution()) _pprintLine_(" computing temperatures"); + #ifdef _HAVE_THERMAL_ + if(isenthalpy==0){ + thermal_core(femmodel); + } + else{ + enthalpy_core(femmodel); + } + #else + _error_("ISSM was not compiled with thermal capabilities. Exiting"); + #endif + } + + if(isdiagnostic){ + if(VerboseSolution()) _pprintLine_(" computing new velocity"); + #ifdef _HAVE_DIAGNOSTIC_ + diagnostic_core(femmodel); + #else + _error_("ISSM was not compiled with diagnostic capabilities. Exiting"); + #endif + } + + if(isprognostic){ + if(VerboseSolution()) _pprintLine_(" computing new thickness"); + prognostic_core(femmodel); + if(VerboseSolution()) _pprintLine_(" updating vertices positions"); + femmodel->UpdateVertexPositionsx(); + } + + if(isgroundingline){ + if(VerboseSolution()) _pprintLine_(" computing new grounding line position"); + #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(isgia){ + if(VerboseSolution()) _pprintLine_(" computing glacial isostatic adjustment"); + #ifdef _HAVE_GIA_ + gia_core(femmodel); + #else + _error_("ISSM was not compiled with gia capabilities. Exiting"); + #endif + + } + + /*unload results*/ + if(save_results){ + if(VerboseSolution()) _pprintLine_(" saving transient results"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsMassBalanceEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaskElementonfloatingiceEnum); + femmodel->RequestedOutputsx(requested_outputs,numoutputs); + + if(isdelta18o){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsMonthlytemperaturesEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsPrecipitationEnum); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalFrictionEnum); + } + if(isgroundingline && (groundingline_migration==SubelementMigrationEnum || groundingline_migration==SubelementMigration2Enum)){ + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,GLlevelsetEnum); + } + + if(VerboseSolution()) _pprintLine_(" saving temporary results"); + OutputResultsx(femmodel->elements, femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,femmodel->results); + } + } + + femmodel->RequestedDependentsx(); + + /*Free ressources:*/ + xDelete(requested_outputs); +} Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 14984) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 14985) @@ -11,5 +11,5 @@ #include #include "../Container/Container.h" -#include "../solutions/solutions.h" +#include "../analyses/analyses.h" #include "../shared/io/io.h" #include "./classes.h" Index: /issm/trunk-jpl/src/c/classes/dakota/DakotaPlugin.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/dakota/DakotaPlugin.cpp (revision 14984) +++ /issm/trunk-jpl/src/c/classes/dakota/DakotaPlugin.cpp (revision 14985) @@ -27,5 +27,5 @@ #include "../../shared/Numerics/types.h" #include "../../shared/MemOps/MemOps.h" -#include "../../solutions/DakotaSpawnCore.h" +#include "../../analyses/DakotaSpawnCore.h" #ifdef _HAVE_DAKOTA_ //only works if dakota library has been compiled in. Index: /issm/trunk-jpl/src/c/main/issm.h =================================================================== --- /issm/trunk-jpl/src/c/main/issm.h (revision 14984) +++ /issm/trunk-jpl/src/c/main/issm.h (revision 14985) @@ -18,5 +18,5 @@ #include "../shared/Enum/Enum.h" #include "../toolkits/toolkits.h" -#include "../solutions/solutions.h" +#include "../analyses/analyses.h" #include "../solvers/solvers.h" #include "../modules/modules.h" Index: /issm/trunk-jpl/src/c/solvers/solver_newton.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_newton.cpp (revision 14984) +++ /issm/trunk-jpl/src/c/solvers/solver_newton.cpp (revision 14985) @@ -3,4 +3,5 @@ */ +#include "./solvers.h" #include "../toolkits/toolkits.h" #include "../classes/objects/objects.h" @@ -8,6 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solutions/solutions.h" -#include "./solvers.h" +#include "../analyses/analyses.h" void solver_newton(FemModel* femmodel){ Index: /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp (revision 14984) +++ /issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp (revision 14985) @@ -3,4 +3,5 @@ */ +#include "./solvers.h" #include "../toolkits/toolkits.h" #include "../classes/objects/objects.h" @@ -8,6 +9,5 @@ #include "../shared/Enum/Enum.h" #include "../modules/modules.h" -#include "../solutions/solutions.h" -#include "./solvers.h" +#include "../analyses/analyses.h" void solver_nonlinear(FemModel* femmodel,bool conserve_loads){ Index: /issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp =================================================================== --- /issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp (revision 14984) +++ /issm/trunk-jpl/src/c/solvers/solver_stokescoupling_nonlinear.cpp (revision 14985) @@ -3,4 +3,5 @@ */ +#include "./solvers.h" #include "../toolkits/toolkits.h" #include "../classes/objects/objects.h" @@ -8,6 +9,5 @@ #include "../shared/io/io.h" #include "../modules/modules.h" -#include "../solutions/solutions.h" -#include "./solvers.h" +#include "../analyses/analyses.h" void solver_stokescoupling_nonlinear(FemModel* femmodel,bool conserve_loads){