Index: /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h =================================================================== --- /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8482) +++ /issm/trunk/src/c/EnumDefinitions/EnumDefinitions.h (revision 8483) @@ -447,5 +447,8 @@ IsDiagnosticEnum, IsPrognosticEnum, - IsThermalEnum + IsThermalEnum, + EnthalpySolutionEnum, + EnthalpyAnalysisEnum, + EnthalpyEnum }; Index: /issm/trunk/src/c/Makefile.am =================================================================== --- /issm/trunk/src/c/Makefile.am (revision 8482) +++ /issm/trunk/src/c/Makefile.am (revision 8483) @@ -463,4 +463,8 @@ ./modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\ ./modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp\ + ./modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp\ + ./modules/ModelProcessorx/Enthalpy/CreateNodesEnthalpy.cpp\ + ./modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp\ + ./modules/ModelProcessorx/Enthalpy/CreateLoadsEnthalpy.cpp\ ./modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp\ ./modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp\ @@ -1121,4 +1125,8 @@ ./modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\ ./modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp\ + ./modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp\ + ./modules/ModelProcessorx/Enthalpy/CreateNodesEnthalpy.cpp\ + ./modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp\ + ./modules/ModelProcessorx/Enthalpy/CreateLoadsEnthalpy.cpp\ ./modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp\ ./modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp\ @@ -1324,4 +1332,5 @@ ./solutions/thermal_core.cpp\ ./solutions/thermal_core_step.cpp\ + ./solutions/enthalpy_core.cpp\ ./solutions/WriteLockFile.cpp\ ./solutions/control_core.cpp\ Index: /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp =================================================================== --- /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8482) +++ /issm/trunk/src/c/modules/EnumToStringx/EnumToStringx.cpp (revision 8483) @@ -391,4 +391,7 @@ case IsPrognosticEnum : return "IsPrognostic"; case IsThermalEnum : return "IsThermal"; + case EnthalpySolutionEnum : return "EnthalpySolution"; + case EnthalpyAnalysisEnum : return "EnthalpyAnalysis"; + case EnthalpyEnum : return "Enthalpy"; default : return "unknown"; Index: /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 8482) +++ /issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp (revision 8483) @@ -74,4 +74,11 @@ break; + case EnthalpyAnalysisEnum: + CreateNodesEnthalpy(pnodes, iomodel,iomodel_handle); + CreateConstraintsEnthalpy(pconstraints,iomodel,iomodel_handle); + CreateLoadsEnthalpy(ploads,iomodel,iomodel_handle); + UpdateElementsEnthalpy(elements,iomodel,iomodel_handle,analysis_counter,analysis_type); + break; + case HydrologyAnalysisEnum: CreateNodesHydrology(pnodes, iomodel,iomodel_handle); Index: /issm/trunk/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp (revision 8482) +++ /issm/trunk/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp (revision 8483) @@ -75,4 +75,7 @@ numdofs=1; } + else if (analysis_type==EnthalpyAnalysisEnum){ + numdofs=1; + } else if (analysis_type==HydrologyAnalysisEnum){ numdofs=1; Index: /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp (revision 8483) +++ /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateConstraintsEnthalpy.cpp (revision 8483) @@ -0,0 +1,60 @@ +/* + * CreateConstraintsEnthalpy.c: + */ + +#include "../../../Container/Container.h" +#include "../../../io/io.h" +#include "../../../toolkits/toolkits.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../shared/shared.h" +#include "../ModelProcessorx.h" + +void CreateConstraintsEnthalpy(Constraints** pconstraints, IoModel* iomodel,FILE* iomodel_handle){ + + /*Intermediary*/ + int i; + int count; + + /*Output*/ + Constraints* constraints = NULL; + Spc* spc = NULL; + + /*Recover pointer: */ + constraints=*pconstraints; + + /*Create constraints if they do not exist yet*/ + if(!constraints) constraints = new Constraints(ConstraintsEnum); + + /*return if 2d mesh*/ + if (iomodel->dim==2) goto cleanup_and_return; + + /*Fetch data: */ + IoModelFetchData(&iomodel->spctemperature,NULL,NULL,iomodel_handle,"spctemperature"); + + /*Initialize counter*/ + count=0; + + /*Create spcs from x,y,z, as well as the spc values on those spcs: */ + for (i=0;inumberofvertices;i++){ + /*keep only this partition's nodes:*/ + if((iomodel->my_vertices[i])){ + + if ((int)iomodel->spctemperature[2*i]){ + + constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,iomodel->spctemperature[2*i+1],EnthalpyAnalysisEnum)); + count++; + + } + + } //if((my_nodes[i]==1)) + } + + /*Free data: */ + xfree((void**)&iomodel->spctemperature); + + cleanup_and_return: + + /*Assign output pointer: */ + *pconstraints=constraints; +} Index: /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateLoadsEnthalpy.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateLoadsEnthalpy.cpp (revision 8483) +++ /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateLoadsEnthalpy.cpp (revision 8483) @@ -0,0 +1,27 @@ +/*! \file CreateLoadsEnthalpy.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../shared/shared.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void CreateLoadsEnthalpy(Loads** ploads, IoModel* iomodel,FILE* iomodel_handle){ + + /*DataSet*/ + Loads* loads=NULL; + + /*Recover pointer: */ + loads=*ploads; + + /*Create loads if they do not exist yet*/ + if(!loads) loads = new Loads(LoadsEnum); + + + /*Assign output pointer: */ + *ploads=loads; +} Index: /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateNodesEnthalpy.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateNodesEnthalpy.cpp (revision 8483) +++ /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/CreateNodesEnthalpy.cpp (revision 8483) @@ -0,0 +1,61 @@ +/* + * CreateNodesEnthalpy.c: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../include/include.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../ModelProcessorx.h" + +void CreateNodesEnthalpy(Nodes** pnodes, IoModel* iomodel,FILE* iomodel_handle){ + + /*Intermediary*/ + int i; + bool continuous_galerkin=true; + + /*DataSets: */ + Nodes* nodes = NULL; + + /*Recover pointer: */ + nodes=*pnodes; + + /*Create nodes if they do not exist yet*/ + if(!nodes) nodes = new Nodes(NodesEnum); + + /*Continuous Galerkin partition of nodes: */ + NodesPartitioning(&iomodel->my_nodes,iomodel->my_elements,iomodel->my_vertices,iomodel,iomodel_handle,continuous_galerkin); + + /*Create nodes and vertices: */ + IoModelFetchData(&iomodel->nodeonbed,NULL,NULL,iomodel_handle,"nodeonbed"); + IoModelFetchData(&iomodel->nodeonsurface,NULL,NULL,iomodel_handle,"nodeonsurface"); + IoModelFetchData(&iomodel->nodeonicesheet,NULL,NULL,iomodel_handle,"nodeonicesheet"); + IoModelFetchData(&iomodel->nodeoniceshelf,NULL,NULL,iomodel_handle,"nodeoniceshelf"); + IoModelFetchData(&iomodel->vertices_type,NULL,NULL,iomodel_handle,"vertices_type"); + IoModelFetchData(&iomodel->nodeonwater,NULL,NULL,iomodel_handle,"nodeonwater"); + + for (i=0;inumberofvertices;i++){ + + if(iomodel->my_vertices[i]){ + + /*Add node to nodes dataset: */ + nodes->AddObject(new Node(iomodel->nodecounter+i+1,i,i+1,i,iomodel,EnthalpyAnalysisEnum)); + + } + } + + /*Clean fetched data: */ + xfree((void**)&iomodel->nodeonbed); + xfree((void**)&iomodel->nodeonsurface); + xfree((void**)&iomodel->nodeonicesheet); + xfree((void**)&iomodel->nodeonwater); + xfree((void**)&iomodel->nodeoniceshelf); + xfree((void**)&iomodel->vertices_type); + + /*Assign output pointer: */ + *pnodes=nodes; +} Index: /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp (revision 8483) +++ /issm/trunk/src/c/modules/ModelProcessorx/Enthalpy/UpdateElementsEnthalpy.cpp (revision 8483) @@ -0,0 +1,80 @@ +/* + * UpdateElementsEnthalpy: + */ + +#include "../../../Container/Container.h" +#include "../../../toolkits/toolkits.h" +#include "../../../io/io.h" +#include "../../../EnumDefinitions/EnumDefinitions.h" +#include "../../../objects/objects.h" +#include "../../../shared/shared.h" +#include "../../MeshPartitionx/MeshPartitionx.h" +#include "../../../include/include.h" +#include "../ModelProcessorx.h" + +void UpdateElementsEnthalpy(Elements* elements, IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type){ + + /*Intermediary*/ + int i; + int counter; + Element* element=NULL; + + /*Now, is the model 3d? otherwise, do nothing: */ + if (iomodel->dim==2)goto cleanup_and_return; + + /*Fetch data needed: */ + IoModelFetchData(&iomodel->elements,NULL,NULL,iomodel_handle,"elements"); + IoModelFetchData(&iomodel->thickness,NULL,NULL,iomodel_handle,"thickness"); + IoModelFetchData(&iomodel->surface,NULL,NULL,iomodel_handle,"surface"); + IoModelFetchData(&iomodel->bed,NULL,NULL,iomodel_handle,"bed"); + IoModelFetchData(&iomodel->drag_coefficient,NULL,NULL,iomodel_handle,"drag_coefficient"); + IoModelFetchData(&iomodel->drag_p,NULL,NULL,iomodel_handle,"drag_p"); + IoModelFetchData(&iomodel->drag_q,NULL,NULL,iomodel_handle,"drag_q"); + IoModelFetchData(&iomodel->elementoniceshelf,NULL,NULL,iomodel_handle,"elementoniceshelf"); + IoModelFetchData(&iomodel->elementonbed,NULL,NULL,iomodel_handle,"elementonbed"); + IoModelFetchData(&iomodel->elementonsurface,NULL,NULL,iomodel_handle,"elementonsurface"); + IoModelFetchData(&iomodel->elementonwater,NULL,NULL,iomodel_handle,"elementonwater"); + IoModelFetchData(&iomodel->elements_type,NULL,NULL,iomodel_handle,"elements_type"); + IoModelFetchData(&iomodel->rheology_B,NULL,NULL,iomodel_handle,"rheology_B"); + IoModelFetchData(&iomodel->rheology_n,NULL,NULL,iomodel_handle,"rheology_n"); + IoModelFetchData(&iomodel->pressure,NULL,NULL,iomodel_handle,"pressure"); + IoModelFetchData(&iomodel->temperature,NULL,NULL,iomodel_handle,"temperature"); + IoModelFetchData(&iomodel->geothermalflux,NULL,NULL,iomodel_handle,"geothermalflux"); + IoModelFetchData(&iomodel->vx,NULL,NULL,iomodel_handle,"vx"); + IoModelFetchData(&iomodel->vy,NULL,NULL,iomodel_handle,"vy"); + IoModelFetchData(&iomodel->vz,NULL,NULL,iomodel_handle,"vz"); + + /*Update elements: */ + counter=0; + for (i=0;inumberofelements;i++){ + if(iomodel->my_elements[i]){ + element=(Element*)elements->GetObjectByOffset(counter); + element->Update(i,iomodel,analysis_counter,analysis_type); //we need i to index into elements. + counter++; + } + } + + cleanup_and_return: + + /*Free data: */ + xfree((void**)&iomodel->elements); + xfree((void**)&iomodel->thickness); + xfree((void**)&iomodel->surface); + xfree((void**)&iomodel->bed); + xfree((void**)&iomodel->drag_coefficient); + xfree((void**)&iomodel->drag_p); + xfree((void**)&iomodel->drag_q); + xfree((void**)&iomodel->elementoniceshelf); + xfree((void**)&iomodel->elementonbed); + xfree((void**)&iomodel->elementonsurface); + xfree((void**)&iomodel->elementonwater); + xfree((void**)&iomodel->elements_type); + xfree((void**)&iomodel->rheology_B); + xfree((void**)&iomodel->rheology_n); + xfree((void**)&iomodel->pressure); + xfree((void**)&iomodel->temperature); + xfree((void**)&iomodel->geothermalflux); + xfree((void**)&iomodel->vx); + xfree((void**)&iomodel->vy); + xfree((void**)&iomodel->vz); +} Index: /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h =================================================================== --- /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 8482) +++ /issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h (revision 8483) @@ -61,4 +61,10 @@ void UpdateElementsThermal(Elements* elements,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type); +/*enthalpy:*/ +void CreateNodesEnthalpy(Nodes** pnodes,IoModel* iomodel,FILE* iomodel_handle); +void CreateConstraintsEnthalpy(Constraints** pconstraints,IoModel* iomodel,FILE* iomodel_handle); +void CreateLoadsEnthalpy(Loads** ploads, IoModel* iomodel, FILE* iomodel_handle); +void UpdateElementsEnthalpy(Elements* elements,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type); + /*hydrology:*/ void CreateNodesHydrology(Nodes** pnodes,IoModel* iomodel,FILE* iomodel_handle); @@ -66,5 +72,4 @@ void CreateLoadsHydrology(Loads** ploads, IoModel* iomodel, FILE* iomodel_handle); void UpdateElementsHydrology(Elements* elements,IoModel* iomodel,FILE* iomodel_handle,int analysis_counter,int analysis_type); - /*melting:*/ Index: /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp =================================================================== --- /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8482) +++ /issm/trunk/src/c/modules/StringToEnumx/StringToEnumx.cpp (revision 8483) @@ -389,4 +389,7 @@ else if (strcmp(name,"IsPrognostic")==0) return IsPrognosticEnum; else if (strcmp(name,"IsThermal")==0) return IsThermalEnum; + else if (strcmp(name,"EnthalpySolution")==0) return EnthalpySolutionEnum; + else if (strcmp(name,"EnthalpyAnalysis")==0) return EnthalpyAnalysisEnum; + else if (strcmp(name,"Enthalpy")==0) return EnthalpyEnum; 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 8482) +++ /issm/trunk/src/c/objects/Elements/Penta.cpp (revision 8483) @@ -1803,4 +1803,245 @@ } + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixEnthalpy {{{1*/ +ElementMatrix* Penta::CreateKMatrixEnthalpy(void){ + + /*compute all stiffness matrices for this element*/ + ElementMatrix* Ke1=CreateKMatrixEnthalpyVolume(); + ElementMatrix* Ke2=CreateKMatrixEnthalpyShelf(); + ElementMatrix* Ke =new ElementMatrix(Ke1,Ke2); + + /*clean-up and return*/ + delete Ke1; + delete Ke2; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixEnthalpyVolume {{{1*/ +ElementMatrix* Penta::CreateKMatrixEnthalpyVolume(void){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + int artdiff; + int i,j,ig,found=0; + double Jdet,u,v,w,um,vm,wm,epsvel; + double gravity,rho_ice,rho_water; + double heatcapacity,thermalconductivity,dt; + double latentheat,kappa=1; + double tau_parameter,diameter; + double xyz_list[NUMVERTICES][3]; + double B[3][numdof]; + double Bprime[3][numdof]; + double B_conduct[3][numdof]; + double B_advec[3][numdof]; + double B_artdiff[2][numdof]; + double Bprime_advec[3][numdof]; + double L[numdof]; + double dh1dh6[3][6]; + double D_scalar_conduct,D_scalar_advec; + double D_scalar_trans,D_scalar_artdiff; + double D[3][3]; + double K[2][2]={0.0}; + double Ke_gaussian_conduct[numdof][numdof]; + double Ke_gaussian_advec[numdof][numdof]; + double Ke_gaussian_artdiff[numdof][numdof]; + double Ke_gaussian_transient[numdof][numdof]; + Tria* tria=NULL; + GaussPenta *gauss=NULL; + + /*Initialize Element matrix*/ + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + gravity=matpar->GetG(); + heatcapacity=matpar->GetHeatCapacity(); + latentheat=matpar->GetLatentHeat(); + thermalconductivity=matpar->GetThermalConductivity(); + this->parameters->FindParam(&dt,DtEnum); + 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); + if (artdiff==2) diameter=MinEdgeLength(xyz_list); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(2,2); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + + /*Conduction: */ + /*Need to change that depending on enthalpy value -> cold or temperate ice: */ + + GetBConduct(&B_conduct[0][0],&xyz_list[0][0],gauss); + + //kappa=GetEnthalpyDiffusivityParameter(rho_ice,heatcapacity,thermalconductivity,latentheat,enthalpy); + D_scalar_conduct=gauss->weight*Jdet*kappa; + if(dt) D_scalar_conduct=D_scalar_conduct*dt; + + D[0][0]=D_scalar_conduct; D[0][1]=0; D[0][2]=0; + D[1][0]=0; D[1][1]=D_scalar_conduct; D[1][2]=0; + D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_conduct; + + TripleMultiply(&B_conduct[0][0],3,numdof,1, + &D[0][0],3,3,0, + &B_conduct[0][0],3,numdof,0, + &Ke_gaussian_conduct[0][0],0); + + /*Advection: */ + + GetBAdvec(&B_advec[0][0],&xyz_list[0][0],gauss); + GetBprimeAdvec(&Bprime_advec[0][0],&xyz_list[0][0],gauss); + + vx_input->GetParameterValue(&u, gauss); + vy_input->GetParameterValue(&v, gauss); + vz_input->GetParameterValue(&w, gauss); + vxm_input->GetParameterValue(&um,gauss); + vym_input->GetParameterValue(&vm,gauss); + vzm_input->GetParameterValue(&wm,gauss); + + D_scalar_advec=gauss->weight*Jdet; + if(dt) D_scalar_advec=D_scalar_advec*dt; + + D[0][0]=D_scalar_advec*(u-um);D[0][1]=0; D[0][2]=0; + D[1][0]=0; D[1][1]=D_scalar_advec*(v-vm);D[1][2]=0; + D[2][0]=0; D[2][1]=0; D[2][2]=D_scalar_advec*(w-wm); + + TripleMultiply(&B_advec[0][0],3,numdof,1, + &D[0][0],3,3,0, + &Bprime_advec[0][0],3,numdof,0, + &Ke_gaussian_advec[0][0],0); + + /*Transient: */ + + if(dt){ + GetNodalFunctionsP1(&L[0], gauss); + D_scalar_trans=gauss->weight*Jdet; + D_scalar_trans=D_scalar_trans; + + TripleMultiply(&L[0],numdof,1,0, + &D_scalar_trans,1,1,0, + &L[0],1,numdof,0, + &Ke_gaussian_transient[0][0],0); + } + else{ + for(i=0;iweight*Jdet/(pow(u-um,2)+pow(v-vm,2)+epsvel); + if(dt) D_scalar_artdiff=D_scalar_artdiff*dt; + K[0][0]=D_scalar_artdiff*pow(u,2); K[0][1]=D_scalar_artdiff*fabs(u)*fabs(v); + K[1][0]=D_scalar_artdiff*fabs(u)*fabs(v);K[1][1]=D_scalar_artdiff*pow(v,2); + + GetBArtdiff(&B_artdiff[0][0],&xyz_list[0][0],gauss); + + TripleMultiply(&B_artdiff[0][0],2,numdof,1, + &K[0][0],2,2,0, + &B_artdiff[0][0],2,numdof,0, + &Ke_gaussian_artdiff[0][0],0); + } + else if(artdiff==2){ + + GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss); + + tau_parameter=GetStabilizationParameter(u-um,v-vm,w-wm,diameter,rho_ice,heatcapacity,thermalconductivity); + + for(i=0;ivalues[i*numdof+j]+=Ke_gaussian_conduct[i][j]+Ke_gaussian_advec[i][j]+Ke_gaussian_transient[i][j]+Ke_gaussian_artdiff[i][j]; + } + + /*Clean up and return*/ + delete gauss; + return Ke; +} +/*}}}*/ +/*FUNCTION Penta::CreateKMatrixEnthalpyShelf {{{1*/ +ElementMatrix* Penta::CreateKMatrixEnthalpyShelf(void){ + + /*Constants*/ + const int numdof=NDOF1*NUMVERTICES; + + /*Intermediaries */ + int i,j,ig; + double mixed_layer_capacity,thermal_exchange_velocity; + double rho_ice,rho_water,heatcapacity; + double Jdet2d,dt; + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[NUMVERTICES2D][3]; + double l1l6[NUMVERTICES]; + double D_scalar; + double Ke_gaussian[numdof][numdof]={0.0}; + GaussPenta *gauss=NULL; + + /*Initialize Element matrix and return if necessary*/ + if (!IsOnBed() || !IsOnShelf()) return NULL; + ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum); + + /*Retrieve all inputs and parameters*/ + this->parameters->FindParam(&dt,DtEnum); + mixed_layer_capacity=matpar->GetMixedLayerCapacity(); + thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); + rho_water=matpar->GetRhoWater(); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + for(i=0;ibegin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); + GetNodalFunctionsP1(&l1l6[0], gauss); + + D_scalar=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity/(rho_ice*heatcapacity); + if(dt) D_scalar=dt*D_scalar; + + TripleMultiply(&l1l6[0],numdof,1,0, + &D_scalar,1,1,0, + &l1l6[0],1,numdof,0, + &Ke_gaussian[0][0],0); + + for(i=0;ivalues[i*numdof+j]+=Ke_gaussian[i][j]; + } + /*Clean up and return*/ delete gauss; @@ -3057,4 +3298,234 @@ /*Clean up and return*/ delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorEnthalpy {{{1*/ +ElementVector* Penta::CreatePVectorEnthalpy(void){ + + /*compute all load vectors for this element*/ + ElementVector* pe1=CreatePVectorEnthalpyVolume(); + ElementVector* pe2=CreatePVectorEnthalpySheet(); + ElementVector* pe3=CreatePVectorEnthalpyShelf(); + ElementVector* pe =new ElementVector(pe1,pe2,pe3); + + /*clean-up and return*/ + delete pe1; + delete pe2; + delete pe3; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorEnthalpyVolume {{{1*/ +ElementVector* Penta::CreatePVectorEnthalpyVolume(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF1; + + /*Intermediaries*/ + int i,j,ig,found=0; + int friction_type,artdiff; + double Jdet,phi,dt; + double rho_ice,heatcapacity; + double thermalconductivity; + double viscosity,temperature; + double tau_parameter,diameter; + double u,v,w; + double scalar_def,scalar_transient; + double temperature_list[NUMVERTICES]; + double xyz_list[NUMVERTICES][3]; + double L[numdof]; + double dh1dh6[3][6]; + double epsilon[6]; + GaussPenta *gauss=NULL; + + /*Initialize Element vector*/ + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + thermalconductivity=matpar->GetThermalConductivity(); + this->parameters->FindParam(&dt,DtEnum); + this->parameters->FindParam(&artdiff,ArtDiffEnum); + 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* temperature_input=NULL; + if (dt) temperature_input=inputs->GetInput(TemperatureEnum); _assert_(inputs); + if (artdiff==2) diameter=MinEdgeLength(xyz_list); + + /* Start looping on the number of gaussian points: */ + gauss=new GaussPenta(2,3); + for (ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetJacobianDeterminant(&Jdet, &xyz_list[0][0],gauss); + GetNodalFunctionsP1(&L[0], gauss); + + this->GetStrainRate3d(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input); + matice->GetViscosity3dStokes(&viscosity,&epsilon[0]); + GetPhi(&phi, &epsilon[0], viscosity); + + scalar_def=phi/(rho_ice)*Jdet*gauss->weight; + if(dt) scalar_def=scalar_def*dt; + + for(i=0;ivalues[i]+=scalar_def*L[i]; + + /* Build transient now */ + if(dt){ + temperature_input->GetParameterValue(&temperature, gauss); + scalar_transient=temperature*Jdet*gauss->weight; + for(i=0;ivalues[i]+=scalar_transient*L[i]; + } + + if(artdiff==2){ + GetNodalFunctionsP1Derivatives(&dh1dh6[0][0],&xyz_list[0][0], gauss); + + vx_input->GetParameterValue(&u, gauss); + vy_input->GetParameterValue(&v, gauss); + vz_input->GetParameterValue(&w, gauss); + + tau_parameter=GetStabilizationParameter(u,v,w,diameter,rho_ice,heatcapacity,thermalconductivity); + + for(i=0;ivalues[i]+=tau_parameter*scalar_def*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]); + if(dt){ + for(i=0;ivalues[i]+=tau_parameter*scalar_transient*(u*dh1dh6[0][i]+v*dh1dh6[1][i]+w*dh1dh6[2][i]); + } + } + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorEnthalpyShelf {{{1*/ +ElementVector* Penta::CreatePVectorEnthalpyShelf(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF1; + + /*Intermediaries */ + int i,j,ig; + double Jdet2d; + 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]; + double l1l6[NUMVERTICES]; + GaussPenta* gauss=NULL; + + /*Initialize Element vector*/ + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + + /* Ice/ocean heat exchange flux on ice shelf base */ + if (!IsOnBed() || !IsOnShelf()) return NULL; + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + for(i=0;iGetMixedLayerCapacity(); + thermal_exchange_velocity=matpar->GetThermalExchangeVelocity(); + rho_water=matpar->GetRhoWater(); + 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); + + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ + gauss=new GaussPenta(0,1,2,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); + GetNodalFunctionsP1(&l1l6[0], gauss); + + pressure_input->GetParameterValue(&pressure,gauss); + h_pmp=meltingpoint-beta*pressure; + + scalar_ocean=gauss->weight*Jdet2d*rho_water*mixed_layer_capacity*thermal_exchange_velocity*(h_pmp)/(rho_ice*heatcapacity); + if(dt) scalar_ocean=dt*scalar_ocean; + + for(i=0;ivalues[i]+=scalar_ocean*l1l6[i]; + } + + /*Clean up and return*/ + delete gauss; + return pe; +} +/*}}}*/ +/*FUNCTION Penta::CreatePVectorEnthalpySheet {{{1*/ +ElementVector* Penta::CreatePVectorEnthalpySheet(void){ + + /*Constants*/ + const int numdof=NUMVERTICES*NDOF1; + + /*Intermediaries */ + int i,j,ig; + int analysis_type,drag_type; + double xyz_list[NUMVERTICES][3]; + double xyz_list_tria[NUMVERTICES2D][3]={0.0}; + double Jdet2d,dt; + double rho_ice,heatcapacity,geothermalflux_value; + double basalfriction,alpha2,vx,vy; + double scalar; + double l1l6[NUMVERTICES]; + Friction* friction=NULL; + GaussPenta* gauss=NULL; + + /* Geothermal flux on ice sheet base and basal friction */ + if (!IsOnBed() || IsOnShelf()) return NULL; + + /*Initialize Element vector*/ + ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters); + + /*Retrieve all inputs and parameters*/ + GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES); + for(i=0;iFindParam(&analysis_type,AnalysisTypeEnum); + rho_ice=matpar->GetRhoIce(); + heatcapacity=matpar->GetHeatCapacity(); + this->parameters->FindParam(&dt,DtEnum); + 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* geothermalflux_input=inputs->GetInput(GeothermalFluxEnum); _assert_(geothermalflux_input); + + /*Build frictoin element, needed later: */ + inputs->GetParameterValue(&drag_type,DragTypeEnum); + if (drag_type!=2)_error_(" non-viscous friction not supported yet!"); + friction=new Friction("3d",inputs,matpar,analysis_type); + + /* Start looping on the number of gauss 2d (nodes on the bedrock) */ + gauss=new GaussPenta(0,1,2,2); + for(ig=gauss->begin();igend();ig++){ + + gauss->GaussPoint(ig); + + GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss); + GetNodalFunctionsP1(&l1l6[0], gauss); + + geothermalflux_input->GetParameterValue(&geothermalflux_value,gauss); + friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum); + vx_input->GetParameterValue(&vx,gauss); + vy_input->GetParameterValue(&vy,gauss); + basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0)); + + scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(rho_ice); + if(dt) scalar=dt*scalar; + + for(i=0;ivalues[i]+=scalar*l1l6[i]; + } + + /*Clean up and return*/ + delete gauss; + delete friction; return pe; } @@ -6634,4 +7105,12 @@ break; + case EnthalpyAnalysisEnum: + /*Initialize mesh velocity*/ + for(i=0;i<6;i++)nodeinputs[i]=0; + this->inputs->AddInput(new PentaVertexInput(VxMeshEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VyMeshEnum,nodeinputs)); + this->inputs->AddInput(new PentaVertexInput(VzMeshEnum,nodeinputs)); + break; + default: /*No update for other solution types*/ Index: /issm/trunk/src/c/objects/Elements/Penta.h =================================================================== --- /issm/trunk/src/c/objects/Elements/Penta.h (revision 8482) +++ /issm/trunk/src/c/objects/Elements/Penta.h (revision 8483) @@ -167,4 +167,7 @@ ElementMatrix* CreateKMatrixDiagnosticVertVolume(void); ElementMatrix* CreateKMatrixDiagnosticVertSurface(void); + ElementMatrix* CreateKMatrixEnthalpy(void); + ElementMatrix* CreateKMatrixEnthalpyVolume(void); + ElementMatrix* CreateKMatrixEnthalpyShelf(void); ElementMatrix* CreateKMatrixMelting(void); ElementMatrix* CreateKMatrixPrognostic(void); @@ -194,4 +197,8 @@ ElementVector* CreatePVectorDiagnosticStokesViscous(void); ElementVector* CreatePVectorDiagnosticStokesShelf(void); + ElementVector* CreatePVectorEnthalpy(void); + ElementVector* CreatePVectorEnthalpyVolume(void); + ElementVector* CreatePVectorEnthalpyShelf(void); + ElementVector* CreatePVectorEnthalpySheet(void); ElementVector* CreatePVectorAdjointStokes(void); ElementVector* CreatePVectorDiagnosticVert(void); Index: /issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp =================================================================== --- /issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp (revision 8482) +++ /issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp (revision 8483) @@ -34,4 +34,7 @@ solutioncore=&thermal_core; break; + case EnthalpySolutionEnum: + solutioncore=&enthalpy_core; + break; case PrognosticSolutionEnum: solutioncore=&prognostic_core; Index: /issm/trunk/src/c/solutions/SolutionConfiguration.cpp =================================================================== --- /issm/trunk/src/c/solutions/SolutionConfiguration.cpp (revision 8482) +++ /issm/trunk/src/c/solutions/SolutionConfiguration.cpp (revision 8483) @@ -58,4 +58,10 @@ analyses[0]=ThermalAnalysisEnum; analyses[1]=MeltingAnalysisEnum; + break; + + case EnthalpySolutionEnum: + numanalyses=1; + analyses=(int*)xmalloc(numanalyses*sizeof(int)); + analyses[0]=EnthalpyAnalysisEnum; break; Index: /issm/trunk/src/c/solutions/enthalpy_core.cpp =================================================================== --- /issm/trunk/src/c/solutions/enthalpy_core.cpp (revision 8483) +++ /issm/trunk/src/c/solutions/enthalpy_core.cpp (revision 8483) @@ -0,0 +1,58 @@ +/*!\file: enthalpy_core.cpp + * \brief: core of the enthalpy solution + */ + +#include "../toolkits/toolkits.h" +#include "../objects/objects.h" +#include "../shared/shared.h" +#include "../EnumDefinitions/EnumDefinitions.h" +#include "./solutions.h" +#include "../modules/modules.h" +#include "../include/include.h" +#include "../solvers/solvers.h" + +void enthalpy_core(FemModel* femmodel){ + + int i; + + /*intermediary*/ + double time; + int nsteps; + double ndt; + double dt; + double melting_offset; + bool control_analysis; + int solution_type; + + //first recover parameters common to all solutions + femmodel->parameters->FindParam(&ndt,NdtEnum); + femmodel->parameters->FindParam(&dt,DtEnum); + femmodel->parameters->FindParam(&control_analysis,ControlAnalysisEnum); + femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); + + /*Compute number of time steps: */ + if((dt==0)|| (ndt==0)){ + dt=0; + nsteps=1; + } + else nsteps=(int)(ndt/dt); + + /*Loop through time: */ + for(i=0;iSetCurrentConfiguration(EnthalpyAnalysisEnum); + solver_linear(femmodel); + + if(solution_type==EnthalpySolutionEnum && !control_analysis){ + _printf_(VerboseSolution()," saving results\n"); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,i+1,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,i+1,time); + InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalMeltingRateEnum,i+1,time); + } + + } +} Index: /issm/trunk/src/c/solutions/solutions.h =================================================================== --- /issm/trunk/src/c/solutions/solutions.h (revision 8482) +++ /issm/trunk/src/c/solutions/solutions.h (revision 8483) @@ -21,4 +21,5 @@ void thermal_core(FemModel* femmodel); void thermal_core_step(FemModel* femmodel,int step, double time); +void enthalpy_core(FemModel* femmodel); void surfaceslope_core(FemModel* femmodel); void bedslope_core(FemModel* femmodel);