Changeset 11347


Ignore:
Timestamp:
02/07/12 14:30:43 (13 years ago)
Author:
seroussi
Message:

added possibility to used enthalpy in transient solutions

Location:
issm/trunk-jpl/src
Files:
1 added
14 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/EnumDefinitions/EnumDefinitions.h

    r11322 r11347  
    161161        ThermalSpctemperatureEnum,
    162162        ThermalStabilizationEnum,
     163        ThermalIsenthalpyEnum,
    163164        ThicknessEnum,
    164165        TimesteppingCflCoefficientEnum,

  • issm/trunk-jpl/src/c/modules/EnumToStringx/EnumToStringx.cpp

    r11322 r11347  
    165165                case ThermalSpctemperatureEnum : return "ThermalSpctemperature";
    166166                case ThermalStabilizationEnum : return "ThermalStabilization";
     167                case ThermalIsenthalpyEnum : return "ThermalIsenthalpy";
    167168                case ThicknessEnum : return "Thickness";
    168169                case TimesteppingCflCoefficientEnum : return "TimesteppingCflCoefficient";

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateParameters.cpp

    r11322 r11347  
    7979        parameters->AddObject(iomodel->CopyConstantObject(TransientIsthermalEnum));
    8080        parameters->AddObject(iomodel->CopyConstantObject(TransientIsgroundinglineEnum));
     81        parameters->AddObject(iomodel->CopyConstantObject(ThermalIsenthalpyEnum));
    8182        parameters->AddObject(iomodel->CopyConstantObject(MaterialsRheologyLawEnum));
    8283        parameters->AddObject(iomodel->CopyConstantObject(AutodiffIsautodiffEnum));

  • issm/trunk-jpl/src/c/modules/ModelProcessorx/ModelProcessorx.cpp

    r11129 r11347  
    2020
    2121        int   i,analysis_type,dim,verbose;
    22         bool  isthermal,isprognostic,isdiagnostic,isgroundingline;
     22        bool  isthermal,isprognostic,isdiagnostic,isgroundingline,isenthalpy;
    2323       
    2424        /*output: */
    … …  
    3939        iomodel->Constant(&verbose,VerboseEnum);
    4040        iomodel->Constant(&isthermal,TransientIsthermalEnum);
     41        iomodel->Constant(&isenthalpy,ThermalIsenthalpyEnum);
    4142        iomodel->Constant(&isprognostic,TransientIsprognosticEnum);
    4243        iomodel->Constant(&isdiagnostic,TransientIsdiagnosticEnum);
    … …  
    5253                if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && dim==2) continue;
    5354                if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && dim==2) continue;
     55                if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && dim==2) continue;
    5456                if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && isthermal==false) continue;
    5557                if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && isthermal==false) continue;
     58                if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && isthermal==false) continue;
     59                if(solution_type==TransientSolutionEnum && analysis_type==ThermalAnalysisEnum && isenthalpy==true) continue;
     60                if(solution_type==TransientSolutionEnum && analysis_type==MeltingAnalysisEnum && isenthalpy==true) continue;
     61                if(solution_type==TransientSolutionEnum && analysis_type==EnthalpyAnalysisEnum && isenthalpy==false) continue;
    5662                if(solution_type==TransientSolutionEnum && analysis_type==PrognosticAnalysisEnum && isprognostic==false && isgroundingline==false) continue;
    5763                if(solution_type==TransientSolutionEnum && analysis_type==DiagnosticHorizAnalysisEnum && isdiagnostic==false) continue;

  • issm/trunk-jpl/src/c/modules/StringToEnumx/StringToEnumx.cpp

    r11322 r11347  
    163163        else if (strcmp(name,"ThermalSpctemperature")==0) return ThermalSpctemperatureEnum;
    164164        else if (strcmp(name,"ThermalStabilization")==0) return ThermalStabilizationEnum;
     165        else if (strcmp(name,"ThermalIsenthalpy")==0) return ThermalIsenthalpyEnum;
    165166        else if (strcmp(name,"Thickness")==0) return ThicknessEnum;
    166167        else if (strcmp(name,"TimesteppingCflCoefficient")==0) return TimesteppingCflCoefficientEnum;

  • issm/trunk-jpl/src/c/objects/Elements/Penta.cpp

    r11341 r11347  
    40854085        double     rho_ice,heatcapacity,geothermalflux_value;
    40864086        double     basalfriction,alpha2,vx,vy;
    4087         double     scalar;
     4087        double     scalar,enthalpy,enthalpyup;
     4088        double     pressure,pressureup;
    40884089        double     basis[NUMVERTICES];
    40894090        Friction*  friction=NULL;
    40904091        GaussPenta* gauss=NULL;
     4092        GaussPenta* gaussup=NULL;
    40914093
    40924094        /* Geothermal flux on ice sheet base and basal friction */
    … …  
    41064108        Input* vy_input=inputs->GetInput(VyEnum);                         _assert_(vy_input);
    41074109        Input* vz_input=inputs->GetInput(VzEnum);                         _assert_(vz_input);
     4110        Input* pressure_input=inputs->GetInput(PressureEnum);             _assert_(pressure_input);
     4111        Input* enthalpy_input=inputs->GetInput(EnthalpyEnum);             _assert_(enthalpy_input);
    41084112        Input* geothermalflux_input=inputs->GetInput(BasalforcingsGeothermalfluxEnum); _assert_(geothermalflux_input);
    41094113
    … …  
    41134117        /* Start looping on the number of gauss 2d (nodes on the bedrock) */
    41144118        gauss=new GaussPenta(0,1,2,2);
     4119        gaussup=new GaussPenta(3,4,5,2);
    41154120        for(ig=gauss->begin();ig<gauss->end();ig++){
    41164121
    41174122                gauss->GaussPoint(ig);
     4123                gaussup->GaussPoint(ig);
    41184124
    41194125                GetTriaJacobianDeterminant(&Jdet2d, &xyz_list_tria[0][0], gauss);
    41204126                GetNodalFunctionsP1(&basis[0], gauss);
    41214127
    4122                 geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
    4123                 friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
    4124                 vx_input->GetInputValue(&vx,gauss);
    4125                 vy_input->GetInputValue(&vy,gauss);
    4126                 basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
    4127                
    4128                 scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
    4129                 if(dt) scalar=dt*scalar;
    4130 
    4131                 for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     4128                enthalpy_input->GetInputValue(&enthalpy,gauss);
     4129                pressure_input->GetInputValue(&pressure,gauss);
     4130//              if(enthalpy>matpar->PureIceEnthalpy(pressure)){
     4131//                      enthalpy_input->GetInputValue(&enthalpyup,gaussup);
     4132//                      pressure_input->GetInputValue(&pressureup,gaussup);
     4133//                      if(enthalpyup>matpar->PureIceEnthalpy(pressureup)){
     4134//                              //do nothing, don't add heatflux
     4135//                      }
     4136//                      else{
     4137//                              //need to change spcenthalpy according to Aschwanden
     4138//                              //NEED TO UPDATE
     4139//                      }
     4140//              }
     4141//              else{
     4142                        geothermalflux_input->GetInputValue(&geothermalflux_value,gauss);
     4143                        friction->GetAlpha2(&alpha2,gauss,VxEnum,VyEnum,VzEnum);
     4144                        vx_input->GetInputValue(&vx,gauss);
     4145                        vy_input->GetInputValue(&vy,gauss);
     4146                        basalfriction=alpha2*(pow(vx,2.0)+pow(vy,2.0));
     4147
     4148                        scalar=gauss->weight*Jdet2d*(basalfriction+geothermalflux_value)/(heatcapacity*rho_ice);
     4149                        if(dt) scalar=dt*scalar;
     4150
     4151                        for(i=0;i<numdof;i++) pe->values[i]+=scalar*basis[i];
     4152//              }
    41324153        }
    41334154
    41344155        /*Clean up and return*/
    41354156        delete gauss;
     4157        delete gaussup;
    41364158        delete friction;
    41374159        return pe;
    … …  
    43094331                        matpar->EnthalpyToThermal(&temperatures[i],&waterfraction[i],values[i],pressure[i]);
    43104332                        if(waterfraction[i]<0) _error_("Negative water fraction found in solution vector");
    4311                         if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector");
     4333                //      if(waterfraction[i]>1) _error_("Water fraction >1 found in solution vector");
    43124334                }
    43134335                       

  • issm/trunk-jpl/src/c/objects/Materials/Matpar.cpp

    r11301 r11347  
    420420        }
    421421        else{
    422                 return 0.1*thermalconductivity/(rho_ice*heatcapacity);
     422                return 1*thermalconductivity/(rho_ice*heatcapacity);
    423423        }
    424424}

  • issm/trunk-jpl/src/c/solutions/AnalysisConfiguration.cpp

    r10287 r11347  
    9595
    9696                case TransientSolutionEnum:
    97                         numanalyses=8;
     97                        numanalyses=9;
    9898                        analyses=(int*)xmalloc(numanalyses*sizeof(int));
    9999                        analyses[0]=DiagnosticHorizAnalysisEnum;
    … …  
    105105                        analyses[6]=MeltingAnalysisEnum;
    106106                        analyses[7]=PrognosticAnalysisEnum;
     107                        analyses[8]=EnthalpyAnalysisEnum;
    107108                        break;
    108109               

  • issm/trunk-jpl/src/c/solutions/transient_core.cpp

    r11245 r11347  
    2424        /*parameters: */
    2525        double finaltime,dt,yts;
    26         bool   control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline;
     26        bool   control_analysis,isdiagnostic,isprognostic,isthermal,isgroundingline,isenthalpy;
    2727        bool   dakota_analysis=false;
    2828        bool   time_adapt=false;
    … …  
    5151        femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
    5252        femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
     53        femmodel->parameters->FindParam(&isenthalpy,ThermalIsenthalpyEnum);
    5354        if(isgroundingline) femmodel->parameters->FindParam(&groundingline_migration,GroundinglineMigrationEnum);
    5455        femmodel->parameters->FindParam(&numoutputs,TransientNumRequestedOutputsEnum);
    … …  
    9192                        _printf_(VerboseSolution(),"   computing temperatures:\n");
    9293                        #ifdef _HAVE_THERMAL_
    93                         thermal_core_step(femmodel,step,time);
     94                        if(isenthalpy==0){
     95                                thermal_core_step(femmodel,step,time);
     96                        }
     97                        else{
     98                                enthalpy_core_step(femmodel,step,time);
     99                        }
    94100                        #else
    95101                        _error_("ISSM was not compiled with thermal capabilities. Exiting");
    … …  
    135141                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time);
    136142                        if(dim==3 && isthermal) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time);
    137                         InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum,step,time);
     143                        if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WaterfractionEnum,step,time);
     144                        if(isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,EnthalpyEnum,step,time);
     145                        if(!isenthalpy) InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BasalforcingsMeltingRateEnum,step,time);
    138146                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceforcingsMassBalanceEnum,step,time);
    139147                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MaskElementonfloatingiceEnum,step,time);

  • issm/trunk-jpl/src/c/solvers/solver_nonlinear.cpp

    r11322 r11347  
    8585                if(count>=max_nonlinear_iterations){
    8686                        _printf_(true,"   maximum number of iterations (%i) exceeded\n",max_nonlinear_iterations);
     87                        converged=true;
     88                InputUpdateFromConstantx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,converged,ConvergedEnum);
     89                InputUpdateFromSolutionx( femmodel->elements,femmodel->nodes, femmodel->vertices, femmodel->loads, femmodel->materials, femmodel->parameters,ug);
    8790                        break;
    8891                }

  • issm/trunk-jpl/src/m/classes/initialization.m

    r11279 r11347  
    6262                                checkfield(md,'initialization.pressure','NaN',1,'size',[md.mesh.numberofvertices 1]);
    6363                        end
    64                         if ismember(EnthalpyAnalysisEnum,analyses),
     64                        if (ismember(EnthalpyAnalysisEnum,analyses) & md.thermal.isenthalpy),
    6565                                checkfield(md,'initialization.waterfraction','>=',0,'size',[md.mesh.numberofvertices 1]);
    6666                        end

  • issm/trunk-jpl/src/m/classes/thermal.m

    r11265 r11347  
    1212                penalty_lock      = 0;
    1313                penalty_factor    = 0;
     14                isenthalpy        = 0;
    1415        end
    1516        methods
    … …  
    4243                        %factor used to compute the values of the penalties: kappa=max(stiffness matrix)*10^penalty_factor
    4344                        obj.penalty_factor=3;
     45
     46                        %Should we use cold ice (default) or enthalpy formulation
     47                        obj.isenthalpy=0;
    4448                end % }}}
    4549                function checkconsistency(obj,md,solution,analyses) % {{{
    … …  
    5256                        if ismember(EnthalpyAnalysisEnum,analyses),
    5357                        checkfield(md,'thermal.spctemperature','<',md.materials.meltingpoint-md.materials.beta*md.materials.rho_ice*md.constants.g*md.geometry.thickness,'message','spctemperature should be below the adjusted melting point');
     58                        checkfield(md,'thermal.isenthalpy','numel',1,'values',[0 1]);
    5459                        end
    5560                end % }}}
    … …  
    6267                        fielddisplay(obj,'penalty_lock','stabilize unstable thermal constraints that keep zigzagging after n iteration (default is 0, no stabilization)');
    6368                        fielddisplay(obj,'penalty_threshold','threshold to declare convergence of thermal solution (default is 0)');
     69                        fielddisplay(obj,'isenthalpy','use an enthalpy formulation to include temperate ice (default is 0)');
    6470
    6571                end % }}}
    … …  
    7177                        WriteData(fid,'object',obj,'fieldname','penalty_lock','format','Integer');
    7278                        WriteData(fid,'object',obj,'fieldname','penalty_factor','format','Double');
     79                        WriteData(fid,'object',obj,'fieldname','isenthalpy','format','Boolean');
    7380                end % }}}
    7481        end

  • issm/trunk-jpl/src/m/solutions/AnalysisConfiguration.m

    r10286 r11347  
    4242
    4343        case TransientSolutionEnum,
    44                 numanalyses=8;
    45                 analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;PrognosticAnalysisEnum];
     44                numanalyses=9;
     45                analyses=[DiagnosticHorizAnalysisEnum;DiagnosticVertAnalysisEnum;DiagnosticHutterAnalysisEnum;SurfaceSlopeAnalysisEnum;BedSlopeAnalysisEnum;ThermalAnalysisEnum;MeltingAnalysisEnum;PrognosticAnalysisEnum;EnthalpyAnalysisEnum];
    4646
    4747        case FlaimSolutionEnum,

  • issm/trunk-jpl/src/m/solutions/transient_core.m

    r10999 r11347  
    1919        isthermal=femmodel.parameters.TransientIsthermal;
    2020        isgroundingline=femmodel.parameters.TransientIsgroundingline;
     21        isenthalpy=femmodel.parameters.ThermalIsenthalpy;
    2122        groundinglinemigration=femmodel.parameters.GroundinglineMigration;
    2223
    … …  
    5758                if (isthermal & dim==3)
    5859                        issmprintf(VerboseSolution,'\n%s',['   computing temperature']);
    59                         femmodel=thermal_core_step(femmodel);
     60                        if (isenthalpy==0),
     61                                femmodel=thermal_core_step(femmodel);
     62                        else
     63                                femmodel=enthalpy_core_step(femmodel);
     64                        end
    6065                end
    6166
    … …  
    9095                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BedEnum,step,time);
    9196                        if (dim==3 & isthermal), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,TemperatureEnum,step,time);end
     97                        if (dim==3 & isenthalpy), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,WaterfractionEnum,step,time);end
     98                        if (dim==3 & isenthalpy), femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,EnthalpyEnum,step,time);end
    9299                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,BasalforcingsMeltingRateEnum,step,time);
    93100                        femmodel.elements=InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceforcingsMassBalanceEnum,step,time);
Note: See TracChangeset for help on using the changeset viewer.