Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8495) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8496) @@ -451,4 +451,5 @@ EnthalpyAnalysisEnum, EnthalpyEnum, + EnthalpyPicardEnum, WaterFractionEnum, ReferenceTemperatureEnum Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8495) +++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8496) @@ -394,4 +394,5 @@ case EnthalpyAnalysisEnum : return "EnthalpyAnalysis"; case EnthalpyEnum : return "Enthalpy"; + case EnthalpyPicardEnum : return "EnthalpyPicard"; case WaterFractionEnum : return "WaterFraction"; case ReferenceTemperatureEnum : return "ReferenceTemperature"; Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8495) +++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8496) @@ -392,4 +392,5 @@ else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; + else if (strcmp(name,"EnthalpyPicard")==0) return EnthalpyPicardEnum; else if (strcmp(name,"WaterFraction")==0) return WaterFractionEnum; else if (strcmp(name,"ReferenceTemperature")==0) return ReferenceTemperatureEnum; Index: /issm/trunk/src/c/objects/Elements/Penta.cpp =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8495) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8496) @@ -612,4 +612,7 @@ Ke=CreateKMatrixThermal(); break; + case EnthalpyAnalysisEnum: + Ke=CreateKMatrixEnthalpy(); + break; case MeltingAnalysisEnum: Ke=CreateKMatrixMelting(); @@ -1834,5 +1837,6 @@ double gravity,rho_ice,rho_water; double heatcapacity,thermalconductivity,dt; - double latentheat,kappa=1; + double pressure,enthalpy; + double latentheat,kappa; double tau_parameter,diameter; double xyz_list[NUMVERTICES][3]; @@ -1870,10 +1874,12 @@ this->parameters->FindParam(&artdiff,ArtDiffEnum); this->parameters->FindParam(&epsvel,EpsVelEnum); - Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); - Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); - Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); - Input* vxm_input=inputs->GetInput(VxMeshEnum); _assert_(vxm_input); - Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input); - Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input); + Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); + Input* enthalpy_input=inputs->GetInput(EnthalpyEnum); _assert_(enthalpy_input); + Input* vx_input=inputs->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=inputs->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=inputs->GetInput(VzEnum); _assert_(vz_input); + Input* vxm_input=inputs->GetInput(VxMeshEnum); _assert_(vxm_input); + Input* vym_input=inputs->GetInput(VyMeshEnum); _assert_(vym_input); + Input* vzm_input=inputs->GetInput(VzMeshEnum); _assert_(vzm_input); if (artdiff==2) diameter=MinEdgeLength(xyz_list); @@ -1891,5 +1897,7 @@ GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); - //kappa=GetEnthalpyDiffusivityParameter(rho_ice,heatcapacity,thermalconductivity,latentheat,enthalpy); + enthalpy_input->GetParameterValue(&enthalpy, gauss); + pressure_input->GetParameterValue(&pressure, gauss); + kappa=matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); D_scalar_conduct=gauss->weight*Jdet*kappa; if(dt) D_scalar_conduct=D_scalar_conduct*dt; @@ -2377,4 +2385,7 @@ pe=CreatePVectorThermal(); break; + case EnthalpyAnalysisEnum: + pe=CreatePVectorEnthalpy(); + break; case MeltingAnalysisEnum: pe=CreatePVectorMelting(); @@ -3412,7 +3423,7 @@ int i,j,ig; double Jdet2d; + double heatcapacity,h_pmp; double mixed_layer_capacity,thermal_exchange_velocity; double rho_ice,rho_water,pressure,dt,scalar_ocean; - double meltingpoint,beta,heatcapacity,h_pmp; double xyz_list[NUMVERTICES][3]; double xyz_list_tria[NUMVERTICES2D][3]; @@ -3434,6 +3445,4 @@ rho_ice=matpar->GetRhoIce(); heatcapacity=matpar->GetHeatCapacity(); - beta=matpar->GetBeta(); - meltingpoint=matpar->GetMeltingPoint(); this->parameters->FindParam(&dt,DtEnum); Input* pressure_input=inputs->GetInput(PressureEnum); _assert_(pressure_input); @@ -3449,5 +3458,5 @@ pressure_input->GetParameterValue(&pressure,gauss); - h_pmp=meltingpoint-beta*pressure; + h_pmp=matpar->PureIceEnthalpy(pressure); scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(h_pmp)/(rho_ice*heatcapacity); @@ -4177,4 +4186,7 @@ GetSolutionFromInputsThermal(solution); } + else if(analysis_type==EnthalpyAnalysisEnum){ + GetSolutionFromInputsEnthalpy(solution); + } else{ _error_("analysis: %i (%s) not supported yet",analysis_type,EnumToStringx(analysis_type)); @@ -4359,4 +4371,35 @@ t_input->GetParameterValue(&temp,gauss); values[i]=temp; + } + + /*Add value to global vector*/ + VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES); + + /*Free ressources:*/ + delete gauss; + xfree((void**)&doflist); +} +/*}}}*/ +/*FUNCTION Penta::GetSolutionFromInputsEnthalpy{{{1*/ +void Penta::GetSolutionFromInputsEnthalpy(Vec solution){ + + const int numdof=NDOF1*NUMVERTICES; + + int i; + int* doflist=NULL; + double values[numdof]; + double enthalpy; + GaussPenta *gauss=NULL; + + /*Get dof list: */ + GetDofList(&doflist,NoneApproximationEnum,GsetEnum); + Input* h_input=inputs->GetInput(EnthalpyEnum); _assert_(h_input); + + gauss=new GaussPenta(); + for(i=0;iGaussVertex(i); + h_input->GetParameterValue(&enthalpy,gauss); + values[i]=enthalpy; } @@ -5042,4 +5085,7 @@ InputUpdateFromSolutionThermal( solution); } + else if (analysis_type==EnthalpyAnalysisEnum){ + InputUpdateFromSolutionEnthalpy( solution); + } else if (analysis_type==MeltingAnalysisEnum){ InputUpdateFromSolutionOneDof(solution,BasalMeltingRateEnum); @@ -6026,4 +6072,5 @@ double xyz_list[NUMVERTICES][3]; double values[numdof]; + double pressure[NUMVERTICES]; double temperatures[numdof]; double waterfraction[numdof]; @@ -6045,20 +6092,13 @@ /*Get all inputs and parameters*/ GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + GetParameterListOnVertices(&pressure[0],PressureEnum); Input* surface_input=inputs->GetInput(SurfaceEnum); _assert_(surface_input); + this->inputs->GetParameterValue(&converged,ConvergedEnum); if(converged){ /*Convert enthalpy into temperature and water fraction*/ -// for(i=0;iEnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]); + this->inputs->AddInput(new PentaVertexInput(EnthalpyEnum,values)); this->inputs->AddInput(new PentaVertexInput(WaterFractionEnum,waterfraction)); @@ -6073,5 +6113,5 @@ break; case PatersonEnum: - B_average=Paterson((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0); + B_average=Paterson((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[5])/6.0); for(i=0;imatice->inputs->AddInput(new PentaVertexInput(RheologyBEnum,B)); @@ -6079,5 +6119,5 @@ case ArrheniusEnum: surface_input->GetParameterAverage(&s_average); - B_average=Arrhenius((values[0]+values[1]+values[2]+values[3]+values[4]+values[5])/6.0, + B_average=Arrhenius((temperatures[0]+temperatures[1]+temperatures[2]+temperatures[3]+temperatures[4]+temperatures[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()); @@ -6090,7 +6130,7 @@ } } -// else{ -// this->inputs->AddInput(new PentaVertexInput(EnthalpyPicardEnum,values)); -// } + else{ + this->inputs->AddInput(new PentaVertexInput(EnthalpyPicardEnum,values)); + } /*Free ressources:*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 8495) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 8496) @@ -225,4 +225,5 @@ void GetSolutionFromInputsDiagnosticVert(Vec solutiong); void GetSolutionFromInputsThermal(Vec solutiong); + void GetSolutionFromInputsEnthalpy(Vec solutiong); double GetStabilizationParameter(double u, double v, double w, double diameter, double rho_ice, double heatcapacity, double thermalconductivity); void GetStrainRate3dPattyn(double* epsilon,double* xyz_list, GaussPenta* gauss, Input* vx_input, Input* vy_input); Index: /issm/trunk/src/c/objects/Materials/Matpar.cpp =================================================================== --- /issm/trunk/src/c/objects/Materials/Matpar.cpp (revision 8495) +++ /issm/trunk/src/c/objects/Materials/Matpar.cpp (revision 8496) @@ -32,4 +32,5 @@ double matpar_beta; double matpar_meltingpoint; + double matpar_referencetemperature; double matpar_mixed_layer_capacity; double matpar_thermal_exchange_velocity; @@ -46,4 +47,5 @@ matpar_beta=iomodel->beta; matpar_meltingpoint=iomodel->meltingpoint; + matpar_referencetemperature=iomodel->referencetemperature; matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity; matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity; @@ -59,4 +61,5 @@ this->beta=matpar_beta; this->meltingpoint=matpar_meltingpoint; + this->referencetemperature=matpar_referencetemperature; this->mixed_layer_capacity=matpar_mixed_layer_capacity; this->thermal_exchange_velocity=matpar_thermal_exchange_velocity; @@ -86,4 +89,5 @@ printf(" beta: %g\n",beta); printf(" meltingpoint: %g\n",meltingpoint); + printf(" referencetemperature: %g\n",referencetemperature); printf(" mixed_layer_capacity: %g\n",mixed_layer_capacity); printf(" thermal_exchange_velocity: %g\n",thermal_exchange_velocity); @@ -106,4 +110,5 @@ printf(" beta: %g\n",beta); printf(" meltingpoint: %g\n",meltingpoint); + printf(" referencetemperature: %g\n",referencetemperature); printf(" mixed_layer_capacity: %g\n",mixed_layer_capacity); printf(" thermal_exchange_velocity: %g\n",thermal_exchange_velocity); @@ -147,4 +152,5 @@ memcpy(marshalled_dataset,&beta,sizeof(beta));marshalled_dataset+=sizeof(beta); memcpy(marshalled_dataset,&meltingpoint,sizeof(meltingpoint));marshalled_dataset+=sizeof(meltingpoint); + memcpy(marshalled_dataset,&referencetemperature,sizeof(referencetemperature));marshalled_dataset+=sizeof(referencetemperature); memcpy(marshalled_dataset,&mixed_layer_capacity,sizeof(mixed_layer_capacity));marshalled_dataset+=sizeof(mixed_layer_capacity); memcpy(marshalled_dataset,&thermal_exchange_velocity,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity); @@ -168,4 +174,5 @@ sizeof(beta)+ sizeof(meltingpoint)+ + sizeof(referencetemperature)+ sizeof(mixed_layer_capacity)+ sizeof(thermal_exchange_velocity)+ @@ -195,4 +202,5 @@ memcpy(&beta,marshalled_dataset,sizeof(beta));marshalled_dataset+=sizeof(beta); memcpy(&meltingpoint,marshalled_dataset,sizeof(meltingpoint));marshalled_dataset+=sizeof(meltingpoint); + memcpy(&referencetemperature,marshalled_dataset,sizeof(referencetemperature));marshalled_dataset+=sizeof(referencetemperature); memcpy(&mixed_layer_capacity,marshalled_dataset,sizeof(mixed_layer_capacity));marshalled_dataset+=sizeof(mixed_layer_capacity); memcpy(&thermal_exchange_velocity,marshalled_dataset,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity); @@ -282,4 +290,8 @@ break; + case ReferenceTemperatureEnum: + this->referencetemperature=constant; + break; + case MixedLayerCapacityEnum: this->mixed_layer_capacity=constant; @@ -347,4 +359,9 @@ double Matpar::GetMeltingPoint(){ return meltingpoint; +} +/*}}}1*/ +/*FUNCTION Matpar::GetReferenceTemperature {{{1*/ +double Matpar::GetReferenceTemperature(){ + return referencetemperature; } /*}}}1*/ @@ -390,2 +407,49 @@ } /*}}}1*/ +/*FUNCTION Matpar::PureIceEnthalpy{{{1*/ +double Matpar::PureIceEnthalpy(double pressure){ + return heatcapacity*(TMeltingPoint(pressure)-referencetemperature); +} +/*}}}1*/ +/*FUNCTION Matpar::GetEnthalpyDiffusionParameter{{{1*/ +double Matpar::GetEnthalpyDiffusionParameter(double enthalpy,double pressure){ + return thermalconductivity/(rho_ice*heatcapacity); +} +/*}}}1*/ +/*FUNCTION Matpar::EnthalpyToThermal {{{1*/ +void Matpar::EnthalpyToThermal(double* ptemperature,double* pwaterfraction,double enthalpy,double pressure){ + + /*Ouput*/ + double temperature,waterfraction; + + if(enthalpy