Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8486) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8487) @@ -450,5 +450,7 @@ EnthalpySolutionEnum, EnthalpyAnalysisEnum, - EnthalpyEnum + EnthalpyEnum, + WaterFractionEnum, + ReferenceTemperatureEnum }; Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8486) +++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8487) @@ -394,4 +394,6 @@ case EnthalpyAnalysisEnum : return "EnthalpyAnalysis"; case EnthalpyEnum : return "Enthalpy"; + case WaterFractionEnum : return "WaterFraction"; + case ReferenceTemperatureEnum : return "ReferenceTemperature"; default : return "unknown"; Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 8486) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateParameters.cpp (revision 8487) @@ -62,4 +62,5 @@ parameters->AddObject(new DoubleParam(BetaEnum,iomodel->beta)); parameters->AddObject(new DoubleParam(MeltingPointEnum,iomodel->meltingpoint)); + parameters->AddObject(new DoubleParam(ReferenceTemperatureEnum,iomodel->referencetemperature)); parameters->AddObject(new DoubleParam(LatentHeatEnum,iomodel->latentheat)); parameters->AddObject(new DoubleParam(HeatCapacityEnum,iomodel->heatcapacity)); Index: /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp (revision 8486) +++ /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp (revision 8487) @@ -40,4 +40,5 @@ IoModelFetchData(&iomodel->pressure,NULL,NULL,iomodel_handle,"pressure"); IoModelFetchData(&iomodel->temperature,NULL,NULL,iomodel_handle,"temperature"); + IoModelFetchData(&iomodel->waterfraction,NULL,NULL,iomodel_handle,"waterfraction"); IoModelFetchData(&iomodel->geothermalflux,NULL,NULL,iomodel_handle,"geothermalflux"); IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); @@ -74,4 +75,5 @@ xfree((void**)&iomodel->pressure); xfree((void**)&iomodel->temperature); + xfree((void**)&iomodel->waterfraction); xfree((void**)&iomodel->geothermalflux); xfree((void**)&iomodel->vx); Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8486) +++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8487) @@ -392,4 +392,6 @@ else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; + else if (strcmp(name,"WaterFraction")==0) return WaterFractionEnum; + else if (strcmp(name,"ReferenceTemperature")==0) return ReferenceTemperatureEnum; else _error_("Enum %s not found",name); Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8486) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8487) @@ -4866,4 +4866,8 @@ for(i=0;i<6;i++)nodeinputs[i]=iomodel->temperature[penta_vertex_ids[i]-1]; this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,nodeinputs)); + } + if (iomodel->waterfraction) { + for(i=0;i<6;i++)nodeinputs[i]=iomodel->waterfraction[penta_vertex_ids[i]-1]; + this->inputs->AddInput(new PentaVertexInput(WaterFractionEnum,nodeinputs)); } if (iomodel->dhdt) { @@ -6010,4 +6014,85 @@ this->inputs->AddInput(new PentaVertexInput(TemperaturePicardEnum,values)); } + + /*Free ressources:*/ + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Penta::InputUpdateFromSolutionEnthalpy {{{1*/ +void Penta::InputUpdateFromSolutionEnthalpy(double* solution){ + + const int numdof=NDOF1*NUMVERTICES; + + bool converged; + int i,rheology_law; + double xyz_list[NUMVERTICES][3]; + double values[numdof]; + double temperatures[numdof]; + double waterfraction[numdof]; + double B[numdof]; + double B_average,s_average; + int* doflist=NULL; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + + /*Use the dof list to index into the solution vector: */ + for(i=0;iGetInput(SurfaceEnum); _assert_(surface_input); + + this->inputs->GetParameterValue(&converged,ConvergedEnum); + if(converged){ + /*Convert enthalpy into temperature and water fraction*/ +// for(i=0;iinputs->AddInput(new PentaVertexInput(EnthalpyEnum,values)); + this->inputs->AddInput(new PentaVertexInput(WaterFractionEnum,waterfraction)); + this->inputs->AddInput(new PentaVertexInput(TemperatureEnum,temperatures)); + + /*Update Rheology only if converged (we must make sure that the temperature is below melting point + * otherwise the rheology could be negative*/ + this->parameters->FindParam(&rheology_law,RheologyLawEnum); + switch(rheology_law){ + case NoneEnum: + /*Do nothing: B is not temperature dependent*/ + break; + case PatersonEnum: + B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0); + for(i=0;imatice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B)); + break; + case ArrheniusEnum: + surface_input->GetParameterAverage(&s_average); + B_average=Arrhenius((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0, + s_average-((xyz_list[0][2]+xyz_list[1][2]+xyz_list[2][2]+xyz_list[3][2]+xyz_list[4][2]+xyz_list[5][2])/6.0), + matice->GetN()); + for(i=0;imatice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B)); + break; + default: + _error_("Rheology law %s not supported yet",EnumToStringx(rheology_law)); + + } + } +// else{ +// this->inputs->AddInput(new PentaVertexInput(EnthalpyPicardEnum,values)); +// } /*Free ressources:*/ @@ -7111,4 +7196,16 @@ this->inputs->AddInput(new PentaVertexInput(VyMeshEnum,nodeinputs)); this->inputs->AddInput(new PentaVertexInput(VzMeshEnum,nodeinputs)); + if (iomodel->temperature && iomodel->waterfraction) { + for(i=0;i<6;i++){ + if(iomodel->temperature[penta_vertex_ids[i]-1] < iomodel->meltingpoint-iomodel->beta*iomodel->pressure[penta_vertex_ids[i]-1]){ + nodeinputs[i]=iomodel->heatcapacity*(iomodel->temperature[penta_vertex_ids[i]-1]-iomodel->referencetemperature); + } + else nodeinputs[i]=iomodel->heatcapacity* + (iomodel->meltingpoint-iomodel->beta*iomodel->pressure[penta_vertex_ids[i]-1]-iomodel->referencetemperature) + +iomodel->latentheat*iomodel->waterfraction[penta_vertex_ids[i]-1]; + } + this->inputs->AddInput(new PentaVertexInput(EnthalpyEnum,nodeinputs)); + } + else _error_("temperature and waterfraction required for the enthalpy solution"); break; Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 8486) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 8487) @@ -245,4 +245,5 @@ void InputUpdateFromSolutionPrognostic(double* solutiong); void InputUpdateFromSolutionThermal( double* solutiong); + void InputUpdateFromSolutionEnthalpy( double* solutiong); void InputUpdateFromSolutionOneDof(double* solutiong,int enum_type); void InputUpdateFromSolutionOneDofCollapsed(double* solutiong,int enum_type); Index: /issm/trunk/src/c/objects/IoModel.cpp =================================================================== --- /issm/trunk/src/c/objects/IoModel.cpp (revision 8486) +++ /issm/trunk/src/c/objects/IoModel.cpp (revision 8487) @@ -69,4 +69,5 @@ xfree((void**)&this->pressure); xfree((void**)&this->temperature); + xfree((void**)&this->waterfraction); xfree((void**)&this->drag_coefficient); xfree((void**)&this->drag_p); @@ -213,4 +214,5 @@ IoModelFetchData(&this->hydro_kn,iomodel_handle,"hydro_kn"); IoModelFetchData(&this->meltingpoint,iomodel_handle,"meltingpoint"); + IoModelFetchData(&this->referencetemperature,iomodel_handle,"referencetemperature"); IoModelFetchData(&this->latentheat,iomodel_handle,"latentheat"); IoModelFetchData(&this->heatcapacity,iomodel_handle,"heatcapacity"); @@ -286,4 +288,5 @@ this->pressure=NULL; this->temperature=NULL; + this->waterfraction=NULL; this->gl_melting_rate=0; this->basal_melting_rate=NULL; @@ -383,4 +386,5 @@ this->beta=0; this->meltingpoint=0; + this->referencetemperature=0; this->latentheat=0; this->heatcapacity=0; Index: /issm/trunk/src/c/objects/IoModel.h =================================================================== --- /issm/trunk/src/c/objects/IoModel.h (revision 8486) +++ /issm/trunk/src/c/objects/IoModel.h (revision 8487) @@ -56,4 +56,5 @@ double* pressure; double* temperature; + double* waterfraction; /*i/o: */ @@ -172,4 +173,5 @@ double hydro_kn; double meltingpoint; + double referencetemperature; double latentheat; double heatcapacity,thermalconductivity;