Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 16811) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.cpp (revision 16812) @@ -188,5 +188,155 @@ }/*}}}*/ ElementVector* EnthalpyAnalysis::CreatePVector(Element* element){/*{{{*/ -_error_("not implemented yet"); + + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorVolume(element); + ElementVector* pe2=CreatePVectorSheet(element); + ElementVector* pe3=CreatePVectorShelf(element); + ElementVector* pe =new ElementVector(pe1,pe2,pe3); + + /*clean-up and return*/ + delete pe1; + delete pe2; + delete pe3; + return pe; +}/*}}}*/ +ElementVector* EnthalpyAnalysis::CreatePVectorVolume(Element* element){/*{{{*/ + + /*Intermediaries*/ + int stabilization; + IssmDouble Jdet,phi,dt; + IssmDouble enthalpy; + IssmDouble kappa,tau_parameter,diameter; + IssmDouble u,v,w; + IssmDouble scalar_def,scalar_transient; + IssmDouble* xyz_list = NULL; + + /*Fetch number of nodes and dof for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + int numvertices = element->GetNumberOfVertices(); + + /*Initialize Element vector*/ + ElementVector* pe = element->NewElementVector(); + IssmDouble* basis = xNew(numnodes); + IssmDouble* dbasis = xNew(3*numnodes); + IssmDouble* pressure = xNew(numvertices); + IssmDouble* enthalpypicard = xNew(numvertices); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinates(&xyz_list); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble thermalconductivity = element->GetMaterialParameter(MaterialsThermalconductivityEnum); + element->FindParam(&dt,TimesteppingTimeStepEnum); + element->FindParam(&stabilization,ThermalStabilizationEnum); + Input* vx_input=element->GetInput(VxEnum); _assert_(vx_input); + Input* vy_input=element->GetInput(VyEnum); _assert_(vy_input); + Input* vz_input=element->GetInput(VzEnum); _assert_(vz_input); + Input* enthalpy_input = NULL; + if(reCast(dt)){enthalpy_input = element->GetInput(EnthalpyEnum); _assert_(enthalpy_input);} + if(stabilization==2){ + element->GetInputListOnVertices(enthalpypicard,EnthalpyPicardEnum); + element->GetInputListOnVertices(pressure,PressureEnum); + diameter=element->MinEdgeLength(xyz_list); + } + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGauss(3); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminant(&Jdet,xyz_list,gauss); + element->NodalFunctions(basis,gauss); + element->ViscousHeating(&phi,xyz_list,gauss,vx_input,vy_input,vz_input); + + scalar_def=phi/rho_ice*Jdet*gauss->weight; + if(reCast(dt)) scalar_def=scalar_def*dt; + + /*TODO: add -beta*laplace T_m(p)*/ + for(int i=0;ivalues[i]+=scalar_def*basis[i]; + + /* Build transient now */ + if(reCast(dt)){ + enthalpy_input->GetInputValue(&enthalpy, gauss); + scalar_transient=enthalpy*Jdet*gauss->weight; + for(int i=0;ivalues[i]+=scalar_transient*basis[i]; + } + + if(stabilization==2){ + element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss); + + vx_input->GetInputValue(&u,gauss); + vy_input->GetInputValue(&v,gauss); + vz_input->GetInputValue(&w,gauss); + kappa = element->EnthalpyDiffusionParameterVolume(numvertices,enthalpypicard,pressure) / rho_ice; + tau_parameter = element->StabilizationParameter(u,v,w,diameter,kappa); + + for(int i=0;ivalues[i]+=tau_parameter*scalar_def*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]); + if(reCast(dt)){ + for(int i=0;ivalues[i]+=tau_parameter*scalar_transient*(u*dbasis[0*3+i]+v*dbasis[1*3+i]+w*dbasis[2*3+i]); + } + } + } + + /*Clean up and return*/ + xDelete(basis); + xDelete(dbasis); + xDelete(pressure); + xDelete(enthalpypicard); + xDelete(xyz_list); + delete gauss; + return pe; + +}/*}}}*/ +ElementVector* EnthalpyAnalysis::CreatePVectorSheet(Element* element){/*{{{*/ + return NULL; +}/*}}}*/ +ElementVector* EnthalpyAnalysis::CreatePVectorShelf(Element* element){/*{{{*/ + + IssmDouble t_pmp,dt,Jdet,scalar_ocean,pressure; + IssmDouble *xyz_list_base = NULL; + + /*Get basal element*/ + if(!element->IsOnBed() || !element->IsFloating()) return NULL; + + /*Fetch number of nodes for this finite element*/ + int numnodes = element->GetNumberOfNodes(); + + /*Initialize vectors*/ + ElementVector* pe = element->NewElementVector(); + IssmDouble* basis = xNew(numnodes); + + /*Retrieve all inputs and parameters*/ + element->GetVerticesCoordinatesBase(&xyz_list_base); + element->FindParam(&dt,TimesteppingTimeStepEnum); + Input* pressure_input=element->GetInput(PressureEnum); _assert_(pressure_input); + IssmDouble gravity = element->GetMaterialParameter(ConstantsGEnum); + IssmDouble rho_water = element->GetMaterialParameter(MaterialsRhoWaterEnum); + IssmDouble rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum); + IssmDouble heatcapacity = element->GetMaterialParameter(MaterialsHeatcapacityEnum); + IssmDouble mixed_layer_capacity= element->GetMaterialParameter(MaterialsMixedLayerCapacityEnum); + IssmDouble thermal_exchange_vel= element->GetMaterialParameter(MaterialsThermalExchangeVelocityEnum); + + /* Start looping on the number of gaussian points: */ + Gauss* gauss=element->NewGaussBase(2); + for(int ig=gauss->begin();igend();ig++){ + gauss->GaussPoint(ig); + + element->JacobianDeterminantBase(&Jdet,xyz_list_base,gauss); + element->NodalFunctions(basis,gauss); + + pressure_input->GetInputValue(&pressure,gauss); + t_pmp=element->TMeltingPoint(pressure); + + scalar_ocean=gauss->weight*Jdet*rho_water*mixed_layer_capacity*thermal_exchange_vel*(t_pmp)/(heatcapacity*rho_ice); + if(reCast(dt)) scalar_ocean=dt*scalar_ocean; + + for(int i=0;ivalues[i]+=scalar_ocean*basis[i]; + } + + /*Clean up and return*/ + delete gauss; + xDelete(basis); + xDelete(xyz_list_base); + return pe; }/*}}}*/ void EnthalpyAnalysis::GetSolutionFromInputs(Vector* solution,Element* element){/*{{{*/ Index: /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h =================================================================== --- /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h (revision 16811) +++ /issm/trunk-jpl/src/c/analyses/EnthalpyAnalysis.h (revision 16812) @@ -23,4 +23,7 @@ ElementMatrix* CreateKMatrix(Element* element); ElementVector* CreatePVector(Element* element); + ElementVector* CreatePVectorVolume(Element* element); + ElementVector* CreatePVectorSheet(Element* element); + ElementVector* CreatePVectorShelf(Element* element); void GetSolutionFromInputs(Vector* solution,Element* element); void InputUpdateFromSolution(IssmDouble* solution,Element* element); Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16811) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 16812) @@ -49,4 +49,6 @@ virtual void DeleteMaterials(void)=0; virtual void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure)=0; + virtual IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure)=0; + virtual IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure)=0; virtual void FindParam(int* pvalue,int paramenum)=0; virtual void FindParam(IssmDouble* pvalue,int paramenum)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16811) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.cpp (revision 16812) @@ -873,4 +873,14 @@ void Penta::EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){ matpar->EnthalpyToThermal(ptemperature,pwaterfraction,enthalpy,pressure); +} +/*}}}*/ +/*FUNCTION Penta::EnthalpyDiffusionParameter{{{*/ +IssmDouble Penta::EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){ + return matpar->GetEnthalpyDiffusionParameter(enthalpy,pressure); +} +/*}}}*/ +/*FUNCTION Penta::EnthalpyDiffusionParameterVolume{{{*/ +IssmDouble Penta::EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){ + return matpar->GetEnthalpyDiffusionParameterVolume(numvertices,enthalpy,pressure); } /*}}}*/ @@ -3781,5 +3791,5 @@ GetInputListOnVertices(&enthalpy[0],EnthalpyPicardEnum); GetInputListOnVertices(&pressure[0],PressureEnum); - kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy,pressure); _assert_(kappa>0.); + kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,enthalpy,pressure); _assert_(kappa>0.); D_scalar_conduct=gauss->weight*Jdet*kappa/rho_ice; @@ -4250,5 +4260,5 @@ GetInputListOnVertices(&enthalpypicard[0],EnthalpyPicardEnum); GetInputListOnVertices(&pressure_p[0],PressureEnum); - kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpypicard,pressure_p); + kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,enthalpypicard,pressure_p); kappa/=rho_ice; tau_parameter=StabilizationParameter(u,v,w,diameter,kappa); @@ -4808,5 +4818,5 @@ else{ enthalpy_input->GetInputDerivativeValue(&d1enthalpy[0],&xyz_list[0][0],gauss); - kappa=matpar->GetEnthalpyDiffusionParameterVolume(enthalpy, pressure); _assert_(kappa>0.); + kappa=matpar->GetEnthalpyDiffusionParameterVolume(NUMVERTICES,&enthalpy[0],&pressure[0]); _assert_(kappa>0.); for(i=0;i<3;i++) vec_heatflux[i]=-kappa*d1enthalpy[i]; } Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16811) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 16812) @@ -204,4 +204,6 @@ ElementVector* CreatePVectorFreeSurfaceBase(void); ElementVector* CreatePVectorL2ProjectionBase(void); + IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); + IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure); void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[6][3],int numpoints); Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16811) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 16812) @@ -84,4 +84,6 @@ void Delta18oParameterization(void){_error_("not implemented yet");}; void EnthalpyToThermal(IssmDouble* ptemperature,IssmDouble* pwaterfraction,IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented yet");}; + IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");}; + IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");}; void FindParam(int* pvalue,int paramenum); void FindParam(IssmDouble* pvalue,int paramenum); Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16811) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 16812) @@ -239,4 +239,6 @@ ElementVector* CreatePVectorL2Projection(void); ElementVector* CreatePVectorL2ProjectionBase(void); + IssmDouble EnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure){_error_("not implemented");}; + IssmDouble EnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){_error_("not implemented");}; IssmDouble GetArea(void); void GetAreaCoordinates(IssmDouble *area_coordinates,IssmDouble xyz_zero[3][3],IssmDouble xyz_list[3][3],int numpoints); Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16811) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.cpp (revision 16812) @@ -417,13 +417,14 @@ /*}}}*/ /*FUNCTION Matpar::GetEnthalpyDiffusionParameterVolume{{{*/ -IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy[6], IssmDouble pressure[6]){ - - IssmDouble lambda; // fraction of cold ice - IssmDouble kappa,kappa_c,kappa_t; //enthalpy conductivities - IssmDouble Hc,Ht; - IssmDouble PIE[6],dHpmp[6]; - int iv; - - for (iv=0; iv<6; iv++){ +IssmDouble Matpar::GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure){ + + int iv; + IssmDouble lambda; // fraction of cold ice + IssmDouble kappa,kappa_c,kappa_t; //enthalpy conductivities + IssmDouble Hc,Ht; + IssmDouble* PIE = xNew(numvertices); + IssmDouble* dHpmp = xNew(numvertices); + + for(iv=0; iv=0)); + for(iv=1; iv=0)); if(allequalsign){ @@ -447,14 +446,18 @@ kappa_t=GetEnthalpyDiffusionParameter(PureIceEnthalpy(0.)+1.,0.); Hc=0.; Ht=0.; - for(iv=0; iv<6;iv++){ + for(iv=0; iv0.); - lambda=Hc/(Hc+Ht); - kappa=1./(lambda/kappa_c + (1.-lambda)/kappa_t); - } + lambda = Hc/(Hc+Ht); + kappa = 1./(lambda/kappa_c + (1.-lambda)/kappa_t); + } + + /*Clean up and return*/ + xDelete(PIE); + xDelete(dHpmp); return kappa; } Index: /issm/trunk-jpl/src/c/classes/Materials/Matpar.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16811) +++ /issm/trunk-jpl/src/c/classes/Materials/Matpar.h (revision 16812) @@ -129,5 +129,5 @@ IssmDouble PureIceEnthalpy(IssmDouble pressure); IssmDouble GetEnthalpyDiffusionParameter(IssmDouble enthalpy,IssmDouble pressure); - IssmDouble GetEnthalpyDiffusionParameterVolume(IssmDouble enthalpy[6], IssmDouble pressure[6]); + IssmDouble GetEnthalpyDiffusionParameterVolume(int numvertices,IssmDouble* enthalpy,IssmDouble* pressure); IssmDouble GetLithosphereShearModulus(); IssmDouble GetLithosphereDensity();