Index: /issm/trunk-jpl/src/c/analyses/HydrologyArmapwAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyArmapwAnalysis.cpp (revision 27463) +++ /issm/trunk-jpl/src/c/analyses/HydrologyArmapwAnalysis.cpp (revision 27463) @@ -0,0 +1,164 @@ +#include "./HydrologyArmapwAnalysis.h" +#include "../toolkits/toolkits.h" +#include "../classes/classes.h" +#include "../shared/shared.h" +#include "../modules/modules.h" + +/*Model processing*/ +void HydrologyArmapwAnalysis::CreateConstraints(Constraints* constraints,IoModel* iomodel){/*{{{*/ + + return; + +}/*}}}*/ +void HydrologyArmapwAnalysis::CreateLoads(Loads* loads, IoModel* iomodel){/*{{{*/ + /*No loads*/ +}/*}}}*/ +void HydrologyArmapwAnalysis::CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr){/*{{{*/ + + return; + +}/*}}}*/ +int HydrologyArmapwAnalysis::DofsPerNode(int** doflist,int domaintype,int approximation){/*{{{*/ + return 0; +}/*}}}*/ +void HydrologyArmapwAnalysis::UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type){/*{{{*/ + + /*Fetch data needed: */ + int hydrology_model,frictionlaw; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want armapw?*/ + if(hydrology_model!=HydrologyarmapwEnum) return; + + /*Add input to elements*/ + iomodel->FetchDataToInput(inputs,elements,"md.mask.ice_levelset",MaskIceLevelsetEnum); + iomodel->FetchDataToInput(inputs,elements,"md.mask.ocean_levelset",MaskOceanLevelsetEnum); + iomodel->FetchDataToInput(inputs,elements,"md.basalforcings.groundedice_melting_rate",BasalforcingsGroundediceMeltingRateEnum); + iomodel->FetchDataToInput(inputs,elements,"md.hydrology.basin_id",HydrologyBasinsIdEnum); + iomodel->FetchDataToInput(inputs,elements,"md.initialization.watercolumn",WatercolumnEnum,0.); + +}/*}}}*/ +void HydrologyArmapwAnalysis::UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum){/*{{{*/ + + /*retrieve some parameters: */ + int hydrology_model; + int numoutputs; + char** requestedoutputs = NULL; + iomodel->FindConstant(&hydrology_model,"md.hydrology.model"); + + /*Now, do we really want Armapw?*/ + if(hydrology_model!=HydrologyarmapwEnum) return; + + parameters->AddObject(new IntParam(HydrologyModelEnum,hydrology_model)); + /*Requested outputs*/ + iomodel->FindConstant(&requestedoutputs,&numoutputs,"md.hydrology.requested_outputs"); + parameters->AddObject(new IntParam(HydrologyNumRequestedOutputsEnum,numoutputs)); + if(numoutputs)parameters->AddObject(new StringArrayParam(HydrologyRequestedOutputsEnum,requestedoutputs,numoutputs)); + iomodel->DeleteData(&requestedoutputs,numoutputs,"md.hydrology.requested_outputs"); + +}/*}}}*/ + +/*Finite Element Analysis*/ +void HydrologyArmapwAnalysis::Core(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +void HydrologyArmapwAnalysis::PreCore(FemModel* femmodel){/*{{{*/ + _error_("not implemented"); +}/*}}}*/ +ElementVector* HydrologyArmapwAnalysis::CreateDVector(Element* element){/*{{{*/ + /*Default, return NULL*/ + return NULL; +}/*}}}*/ +ElementMatrix* HydrologyArmapwAnalysis::CreateJacobianMatrix(Element* element){/*{{{*/ + _error_("Not implemented"); +}/*}}}*/ +ElementMatrix* HydrologyArmapwAnalysis::CreateKMatrix(Element* element){/*{{{*/ + _error_("Not implemented"); +}/*}}}*/ +ElementVector* HydrologyArmapwAnalysis::CreatePVector(Element* element){/*{{{*/ + _error_("Not implemented"); +}/*}}}*/ +void HydrologyArmapwAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ + _error_("Not implemented"); +}/*}}}*/ +void HydrologyArmapwAnalysis::GradientJ(Vector* gradient,Element* element,int control_type,int control_interp,int control_index){/*{{{*/ + _error_("Not implemented"); +}/*}}}*/ +void HydrologyArmapwAnalysis::InputUpdateFromSolution(IssmDouble* solution,Element* element){/*{{{*/ + _error_("Not implemented"); +}/*}}}*/ +void HydrologyArmapwAnalysis::UpdateConstraints(FemModel* femmodel){/*{{{*/ + _error_("Not implemented"); +}/*}}}*/ + +/*Additional methods*/ +void HydrologyArmapwAnalysis::UpdateSubglacialWaterPressure(FemModel* femmodel){/*{{{*/ + + /*Get time parameters*/ + IssmDouble time,dt,starttime,tstep_arma; + femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->parameters->FindParam(&dt,TimesteppingTimeStepEnum); + femmodel->parameters->FindParam(&starttime,TimesteppingStartTimeEnum); + femmodel->parameters->FindParam(&tstep_arma,HydrologyarmaTimestepEnum); + + /*Determine if this is a time step for the ARMA model*/ + bool isstepforarma = false; + + #ifndef _HAVE_AD_ + if((fmod(time,tstep_arma)parameters->FindParam(&numbasins,HydrologyNumBasinsEnum); + femmodel->parameters->FindParam(&numparams,HydrologyarmaNumParamsEnum); + femmodel->parameters->FindParam(&numbreaks,HydrologyarmaNumBreaksEnum); + femmodel->parameters->FindParam(&arorder,HydrologyarmaarOrderEnum); + femmodel->parameters->FindParam(&maorder,HydrologyarmamaOrderEnum); + IssmDouble* datebreaks = NULL; + IssmDouble* arlagcoefs = NULL; + IssmDouble* malagcoefs = NULL; + IssmDouble* monthlyfactors = NULL; + IssmDouble* polyparams = NULL; + + femmodel->parameters->FindParam(&datebreaks,&M,&N,HydrologyarmadatebreaksEnum); _assert_(M==numbasins); _assert_(N==max(numbreaks,1)); + femmodel->parameters->FindParam(&polyparams,&M,&N,HydrologyarmapolyparamsEnum); _assert_(M==numbasins); _assert_(N==(numbreaks+1)*numparams); + femmodel->parameters->FindParam(&arlagcoefs,&M,&N,HydrologyarmaarlagcoefsEnum); _assert_(M==numbasins); _assert_(N==arorder); + femmodel->parameters->FindParam(&malagcoefs,&M,&N,HydrologyarmamalagcoefsEnum); _assert_(M==numbasins); _assert_(N==maorder); + femmodel->parameters->FindParam(&monthlyfactors,&M,&N,HydrologyarmaMonthlyFactorsEnum); _assert_(M==numbasins); _assert_(N==12); + + femmodel->parameters->FindParam(&isstochastic,StochasticForcingIsStochasticForcingEnum); + if(isstochastic){ + int numstochasticfields; + int* stochasticfields; + femmodel->parameters->FindParam(&numstochasticfields,StochasticForcingNumFieldsEnum); + femmodel->parameters->FindParam(&stochasticfields,&N,StochasticForcingFieldsEnum); _assert_(N==numstochasticfields); + for(int i=0;i(stochasticfields); + } + + /*Loop over each element to compute Subglacial Water Pressure at vertices*/ + for(Object* &object:femmodel->elements->objects){ + Element* element = xDynamicCast(object); + /*Compute ARMA perturbation values*/ + element->ArmaProcess(isstepforarma,arorder,maorder,numparams,numbreaks,tstep_arma,polyparams,arlagcoefs,malagcoefs,datebreaks,ispwstochastic,HydrologyarmapwEnum); + /*Compute subglacial water pressure with the ARMA perturbation*/ + element->SubglacialWaterPressure(FrictionWaterPressureEnum); + /*Scale with monthly factors*/ + element->MonthlyFactorBasin(monthlyfactors,HydrologyarmapwEnum); + } + + /*Cleanup*/ + xDelete(arlagcoefs); + xDelete(malagcoefs); + xDelete(polyparams); + xDelete(datebreaks); + xDelete(monthlyfactors); +}/*}}}*/ + Index: /issm/trunk-jpl/src/c/analyses/HydrologyArmapwAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/HydrologyArmapwAnalysis.h (revision 27463) +++ /issm/trunk-jpl/src/c/analyses/HydrologyArmapwAnalysis.h (revision 27463) @@ -0,0 +1,37 @@ +/*! \file HydrologyArmapwAnalysis.h + * \brief: header file for generic external result object + */ + +#ifndef _HydrologyArmapwAnalysis_ +#define _HydrologyArmapwAnalysis_ + +/*Headers*/ +#include "./Analysis.h" + +class HydrologyArmapwAnalysis: public Analysis{ + + public: + /*Model processing*/ + void CreateConstraints(Constraints* constraints,IoModel* iomodel); + void CreateLoads(Loads* loads, IoModel* iomodel); + void CreateNodes(Nodes* nodes,IoModel* iomodel,bool isamr=false); + int DofsPerNode(int** doflist,int domaintype,int approximation); + void UpdateElements(Elements* elements,Inputs* inputs,IoModel* iomodel,int analysis_counter,int analysis_type); + void UpdateParameters(Parameters* parameters,IoModel* iomodel,int solution_enum,int analysis_enum); + + /*Finite element Analysis*/ + void Core(FemModel* femmodel); + void PreCore(FemModel* femmodel); + ElementVector* CreateDVector(Element* element); + ElementMatrix* CreateJacobianMatrix(Element* element); + ElementMatrix* CreateKMatrix(Element* element); + ElementVector* CreatePVector(Element* element); + void GetSolutionFromInputs(Vector* solution,Element* element); + void GradientJ(Vector* gradient,Element* element,int control_type,int control_interp,int control_index); + void InputUpdateFromSolution(IssmDouble* solution,Element* element); + void UpdateConstraints(FemModel* femmodel); + + /*Intermediaries*/ + void UpdateSubglacialWaterPressure(FemModel* femmodel); +}; +#endif