Changeset 7640


Ignore:
Timestamp:
03/15/11 09:01:44 (14 years ago)
Author:
Eric.Larour
Message:

New hydrology solution, with new modules, new element routines, new data for materials, new model processing, and new solution configurations

Location:
issm/trunk/src/c
Files:
9 added
15 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/c/Makefile.am

    r7447 r7640  
    395395                                        ./modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\
    396396                                        ./modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp\
     397                                        ./modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp\
     398                                        ./modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp\
     399                                        ./modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp\
     400                                        ./modules/ModelProcessorx/Hydrology/CreateLoadsHydrology.cpp\
    397401                                        ./modules/ModelProcessorx/Melting/UpdateElementsMelting.cpp\
    398402                                        ./modules/ModelProcessorx/Melting/CreateNodesMelting.cpp\
    … …  
    420424                                        ./modules/OutputResultsx/OutputResultsx.h\
    421425                                        ./modules/OutputResultsx/OutputResultsx.cpp\
     426                                        ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.h\
     427                                        ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsx.cpp\
     428                                        ./modules/PointCloudFindNeighborsx/PointCloudFindNeighborsxt.cpp\
    422429                                        ./modules/MinVelx/MinVelx.h\
    423430                                        ./modules/MinVelx/MinVelx.cpp\
    … …  
    970977                                        ./modules/ModelProcessorx/Thermal/CreateConstraintsThermal.cpp\
    971978                                        ./modules/ModelProcessorx/Thermal/CreateLoadsThermal.cpp\
     979                                        ./modules/ModelProcessorx/Hydrology/UpdateElementsHydrology.cpp\
     980                                        ./modules/ModelProcessorx/Hydrology/CreateNodesHydrology.cpp\
     981                                        ./modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp\
     982                                        ./modules/ModelProcessorx/Hydrology/CreateLoadsHydrology.cpp\
    972983                                        ./modules/ModelProcessorx/Melting/UpdateElementsMelting.cpp\
    973984                                        ./modules/ModelProcessorx/Melting/CreateNodesMelting.cpp\
    … …  
    11711182                                        ./solutions/surfaceslope_core.cpp\
    11721183                                        ./solutions/bedslope_core.cpp\
     1184                                        ./solutions/hydrology_core.cpp\
     1185                                        ./solutions/hydrology_core_step.cpp\
    11731186                                        ./solutions/transient2d_core.cpp\
    11741187                                        ./solutions/transient3d_core.cpp\
    … …  
    11821195                                        ./solvers/solver_linear.cpp\
    11831196                                        ./solvers/solver_adjoint_linear.cpp\
    1184                                         ./solvers/solver_diagnostic_nonlinear.cpp\
     1197                                        ./solvers/solver_nonlinear.cpp\
    11851198                                        ./solvers/solver_stokescoupling_nonlinear.cpp\
    11861199                                        ./solvers/solver_thermal_nonlinear.cpp\

  • issm/trunk/src/c/modules/ModelProcessorx/CreateDataSets.cpp

    r6412 r7640  
    7272                        UpdateElementsThermal(elements,iomodel,iomodel_handle,analysis_counter,analysis_type);
    7373                        break;
     74               
     75                case HydrologyAnalysisEnum:
     76                        CreateNodesHydrology(pnodes, iomodel,iomodel_handle);
     77                        CreateConstraintsHydrology(pconstraints,iomodel,iomodel_handle);
     78                        CreateLoadsHydrology(ploads,iomodel,iomodel_handle);
     79                        UpdateElementsHydrology(elements,iomodel,iomodel_handle,analysis_counter,analysis_type);
     80                        break;
    7481
    7582                case MeltingAnalysisEnum:

  • issm/trunk/src/c/modules/ModelProcessorx/DistributeNumDofs.cpp

    r6412 r7640  
    7575                numdofs=1;
    7676        }
     77        else if (analysis_type==HydrologyAnalysisEnum){
     78                numdofs=1;
     79        }
    7780        else if (analysis_type==MeltingAnalysisEnum){
    7881                numdofs=1;

  • issm/trunk/src/c/modules/ModelProcessorx/ModelProcessorx.h

    r7159 r7640  
    6161void    UpdateElementsThermal(Elements* elements,IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_counter,int analysis_type);
    6262
     63/*hydrology:*/
     64void    CreateNodesHydrology(Nodes** pnodes,IoModel* iomodel,ConstDataHandle iomodel_handle);
     65void    CreateConstraintsHydrology(Constraints** pconstraints,IoModel* iomodel,ConstDataHandle iomodel_handle);
     66void  CreateLoadsHydrology(Loads** ploads, IoModel* iomodel, ConstDataHandle iomodel_handle);
     67void    UpdateElementsHydrology(Elements* elements,IoModel* iomodel,ConstDataHandle iomodel_handle,int analysis_counter,int analysis_type);
     68
     69
    6370/*melting:*/
    6471void    CreateNodesMelting(Nodes** pnodes,IoModel* iomodel,ConstDataHandle iomodel_handle);

  • issm/trunk/src/c/modules/modules.h

    r7079 r7640  
    7474#include "./OutputRiftsx/OutputRiftsx.h"
    7575#include "./PenaltyConstraintsx/PenaltyConstraintsx.h"
     76#include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h"
    7677#include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h"
    7778#include "./Qmux/Qmux.h"

  • issm/trunk/src/c/objects/Bamg/Geometry.cpp

    r6412 r7640  
    527527                        /*if there is a vertex found that is to close to vertices[i] -> error*/
    528528                        if( v && Norme1(v->r - vertices[i].r) < eps ){
    529                                 _error_("two points of the geometry are very closed to each other");
     529                                //_error_("two points of the geometry are very closed to each other");
     530                                printf("%g %g\n",v->ReferenceNumber,vertices[i].ReferenceNumber);
    530531                        }
    531532

  • issm/trunk/src/c/objects/Elements/Tria.cpp

    r7391 r7640  
    405405                        Ke=CreateKMatrixPrognostic();
    406406                        break;
     407                case HydrologyAnalysisEnum:
     408                        Ke=CreateKMatrixHydrology();
     409                        break;
    407410                case BalancedthicknessAnalysisEnum:
    408411                        Ke=CreateKMatrixBalancedthickness();
    … …  
    11821185}
    11831186/*}}}*/
     1187/*FUNCTION Tria::CreateKMatrixHydrology {{{1*/
     1188ElementMatrix* Tria::CreateKMatrixHydrology(void){
     1189
     1190        /*Constants*/
     1191        const int    numdof=NDOF1*NUMVERTICES;
     1192
     1193        /*Intermediaries */
     1194        int        i,j,ig;
     1195        double     Jdettria,dt;
     1196        double     xyz_list[NUMVERTICES][3];
     1197        double     L[NUMVERTICES];
     1198        double     dL[2][NUMVERTICES];
     1199        double     dLx[NUMVERTICES];
     1200        double     dLy[NUMVERTICES];
     1201        double     K[2];
     1202        double     Ke_gg_gaussian1[numdof][numdof]  ={0.0};
     1203        double     Ke_gg_gaussian2[numdof][numdof]  ={0.0};
     1204        double     Ke_gg_gaussian3[numdof][numdof]  ={0.0};
     1205        double     gamma;
     1206        double     DL_scalar1;
     1207        double     DL_scalar2;
     1208        double     DL_scalar3;
     1209        GaussTria *gauss=NULL;
     1210
     1211        /*Initialize Element matrix*/
     1212        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     1213
     1214        /*Retrieve all inputs and parameters*/
     1215        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     1216        this->parameters->FindParam(&dt,DtEnum);
     1217       
     1218        /*retrieve material parameters: */
     1219        gamma=matpar->GetGamma();
     1220
     1221        /* Start  looping on the number of gaussian points: */
     1222        gauss=new GaussTria(2);
     1223        for (ig=gauss->begin();ig<gauss->end();ig++){
     1224
     1225                gauss->GaussPoint(ig);
     1226
     1227                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     1228                GetNodalFunctions(&L[0],gauss);
     1229                GetNodalFunctionsDerivatives(&dL[0][0],&xyz_list[0][0],gauss);
     1230
     1231                for(i=0;i<NUMVERTICES;i++)dLx[i]=dL[0][i];
     1232                for(i=0;i<NUMVERTICES;i++)dLy[i]=dL[1][i];
     1233
     1234                /*Get K parameter: */
     1235                GetHydrologyK(&K[0],&xyz_list[0][0],gauss);
     1236       
     1237                if(dt){
     1238                        DL_scalar1=gauss->weight*Jdettria;
     1239                        TripleMultiply( &L[0],1,numdof,1,
     1240                                        &DL_scalar1,1,1,0,
     1241                                        &L[0],1,numdof,0,
     1242                                        &Ke_gg_gaussian1[0][0],0);
     1243                }
     1244
     1245                if(dt){
     1246                        DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt; //don't forget the -
     1247                        DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt; //don't forget the -
     1248
     1249                }
     1250                else{
     1251                        DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]; //don't forget the -
     1252                        DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]; //don't forget the -
     1253                }
     1254               
     1255                TripleMultiply( &dLx[0],1,numdof,1, &DL_scalar2,1,1,0, &L[0],1,numdof,0, &Ke_gg_gaussian2[0][0],0);
     1256                TripleMultiply( &dLy[0],1,numdof,1, &DL_scalar3,1,1,0, &L[0],1,numdof,0, &Ke_gg_gaussian3[0][0],0);
     1257
     1258                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian1[i][j];
     1259                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian2[i][j];
     1260                for(i=0;i<numdof;i++) for(j=0;j<numdof;j++) Ke->values[i*numdof+j]+=Ke_gg_gaussian3[i][j];
     1261        }
     1262
     1263        /*Clean up and return*/
     1264        delete gauss;
     1265        return Ke;
     1266}
     1267/*}}}*/
    11841268/*FUNCTION Tria::CreateKMatrixPrognostic {{{1*/
    11851269ElementMatrix* Tria::CreateKMatrixPrognostic(void){
    … …  
    15301614                case PrognosticAnalysisEnum:
    15311615                        pe=CreatePVectorPrognostic();
     1616                        break;
     1617                case HydrologyAnalysisEnum:
     1618                        pe=CreatePVectorHydrology();
    15321619                        break;
    15331620                case BalancedthicknessAnalysisEnum:
    … …  
    22932380}
    22942381/*}}}*/
     2382/*FUNCTION Tria::CreatePVectorHydrology {{{1*/
     2383ElementVector* Tria::CreatePVectorHydrology(void){
     2384
     2385        /*Constants*/
     2386        const int    numdof=NDOF1*NUMVERTICES;
     2387
     2388        /*Intermediaries */
     2389        int        i,j,ig;
     2390        double     Jdettria,dt;
     2391        double     melting_g;
     2392        double     old_watercolumn_g;
     2393        double     xyz_list[NUMVERTICES][3];
     2394        double     L[NUMVERTICES];
     2395        GaussTria* gauss=NULL;
     2396
     2397        /*Initialize Element vector*/
     2398        ElementVector* pe=new ElementVector(nodes,NUMVERTICES,this->parameters);
     2399
     2400        /*Retrieve all inputs and parameters*/
     2401        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     2402        this->parameters->FindParam(&dt,DtEnum);
     2403        Input* melting_input=inputs->GetInput(MeltingRateEnum);           _assert_(melting_input);
     2404        Input* old_watercolumn_input=inputs->GetInput(OldWaterColumnEnum);           _assert_(old_watercolumn_input);
     2405
     2406        /*Initialize melting_correction_g to 0, do not forget!:*/
     2407        /* Start  looping on the number of gaussian points: */
     2408        gauss=new GaussTria(2);
     2409        for(ig=gauss->begin();ig<gauss->end();ig++){
     2410
     2411                gauss->GaussPoint(ig);
     2412
     2413                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     2414                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
     2415
     2416                melting_input->GetParameterValue(&melting_g,gauss);
     2417                old_watercolumn_input->GetParameterValue(&old_watercolumn_g,gauss);
     2418       
     2419
     2420               
     2421                if(dt)for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*melting_g)*L[i];
     2422                else  for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*melting_g*L[i];
     2423
     2424        }
     2425               
     2426        /*Clean up and return*/
     2427        delete gauss;
     2428        return pe;
     2429}
     2430/*}}}*/
    22952431/*FUNCTION Tria::CreatePVectorPrognostic_CG {{{1*/
    22962432ElementVector* Tria::CreatePVectorPrognostic_CG(void){
    … …  
    28953031        else if (analysis_type==DiagnosticHutterAnalysisEnum)
    28963032         GetSolutionFromInputsDiagnosticHutter(solution);
     3033        else if (analysis_type==HydrologyAnalysisEnum)
     3034         GetSolutionFromInputsHydrology(solution);
    28973035        else
    28983036         _error_("analysis: %s not supported yet",EnumToString(analysis_type));
    … …  
    29693107                values[i*NDOF2+0]=vx;
    29703108                values[i*NDOF2+1]=vy;
     3109        }
     3110
     3111        VecSetValues(solution,numdof,doflist,(const double*)values,INSERT_VALUES);
     3112
     3113        /*Free ressources:*/
     3114        delete gauss;
     3115        xfree((void**)&doflist);
     3116}
     3117/*}}}*/
     3118/*FUNCTION Tria::GetSolutionFromInputsHydrology{{{1*/
     3119void  Tria::GetSolutionFromInputsHydrology(Vec solution){
     3120
     3121        const int    numdof=NDOF1*NUMVERTICES;
     3122
     3123        int i,dummy;
     3124        int*         doflist=NULL;
     3125        double       watercolumn;
     3126        double       values[numdof];
     3127        GaussTria*   gauss=NULL;
     3128
     3129        /*Get dof list: */
     3130        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     3131
     3132        /*Get inputs*/
     3133        Input* watercolumn_input=inputs->GetInput(WaterColumnEnum); _assert_(watercolumn_input);
     3134
     3135        /*Ok, we have watercolumn values, fill in watercolumn array: */
     3136        /*P1 element only for now*/
     3137        gauss=new GaussTria();
     3138        for(i=0;i<NUMVERTICES;i++){
     3139
     3140                gauss->GaussVertex(i);
     3141
     3142                /*Recover watercolumn*/
     3143                watercolumn_input->GetParameterValue(&watercolumn,gauss);
     3144                values[i]=watercolumn;
     3145                if((values[i])<0)values[i]=0;
    29713146        }
    29723147
    … …  
    30153190                }
    30163191        }
     3192}
     3193/*}}}*/
     3194/*FUNCTION Tria::GetHydrologyK {{{1*/
     3195void Tria::GetHydrologyK(double* K,double* xyz_list,GaussTria* gauss){
     3196
     3197
     3198        /*material parameters: */
     3199        double rho_ice;
     3200        double rho_water;
     3201        double g;
     3202        double dsdx, dsdy;
     3203        double dbdx, dbdy;
     3204        double surface_slope;
     3205        double kn;
     3206        double w;
     3207
     3208        Input* surfaceslopex_input=NULL;
     3209        Input* surfaceslopey_input=NULL;
     3210        Input* bedslopex_input=NULL;
     3211        Input* bedslopey_input=NULL;
     3212        Input* watercolumn_input=NULL;
     3213
     3214       
     3215        /*recover parameters and inputs: */
     3216        rho_ice=matpar->GetRhoIce();
     3217        rho_water=matpar->GetRhoWater();
     3218        g=matpar->GetG();
     3219        kn=matpar->GetKn();
     3220
     3221        surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
     3222        surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
     3223        bedslopex_input=inputs->GetInput(BedSlopeXEnum); _assert_(bedslopex_input);
     3224        bedslopey_input=inputs->GetInput(BedSlopeYEnum); _assert_(bedslopey_input);
     3225        watercolumn_input=inputs->GetInput(WaterColumnEnum); _assert_(watercolumn_input);
     3226
     3227        /*recover slopes: */
     3228        surfaceslopex_input->GetParameterValue(&dsdx,gauss);
     3229        surfaceslopey_input->GetParameterValue(&dsdy,gauss);
     3230        bedslopex_input->GetParameterValue(&dbdx,gauss);
     3231        bedslopey_input->GetParameterValue(&dbdy,gauss);
     3232
     3233        /*magnitude of surface slope: */
     3234        surface_slope=sqrt( pow(dsdx,2) + pow(dsdy,2));
     3235
     3236        /*recover water column: */
     3237        watercolumn_input->GetParameterValue(&w,gauss);
     3238
     3239        K[0]=pow(w,2)*(rho_ice*g*dsdx+(rho_water/rho_ice-1)*rho_ice*g*dbdx) - rho_ice * g * kn* w * (dsdx - dbdx ) * surface_slope;
     3240        K[1]=pow(w,2)*(rho_ice*g*dsdy+(rho_water/rho_ice-1)*rho_ice*g*dbdy) - rho_ice * g * kn *w * (dsdy - dbdy ) * surface_slope;
     3241
    30173242}
    30183243/*}}}*/
    … …  
    36623887                this->inputs->AddInput(new TriaVertexInput(MeltingRateCorrectionEnum,nodeinputs));
    36633888        }
     3889        if (iomodel->watercolumn){
     3890                for(i=0;i<3;i++)nodeinputs[i]=iomodel->watercolumn[tria_vertex_ids[i]-1];
     3891                this->inputs->AddInput(new TriaVertexInput(WaterColumnEnum,nodeinputs));
     3892                this->inputs->AddInput(new TriaVertexInput(OldWaterColumnEnum,nodeinputs));
     3893        }
     3894
    36643895        if (iomodel->accumulation_rate) {
    36653896                for(i=0;i<3;i++)nodeinputs[i]=iomodel->accumulation_rate[tria_vertex_ids[i]-1]/iomodel->yts;
    … …  
    37924023                        InputUpdateFromSolutionPrognostic(solution);
    37934024                        break;
     4025                case HydrologyAnalysisEnum:
     4026                        InputUpdateFromSolutionHydrology(solution);
     4027                        break;
    37944028                case BalancedthicknessAnalysisEnum:
    37954029                        InputUpdateFromSolutionOneDof(solution,ThicknessEnum);
    … …  
    41004334}
    41014335/*}}}*/
     4336/*FUNCTION Tria::InputUpdateFromSolutionHydrology{{{1*/
     4337void  Tria::InputUpdateFromSolutionHydrology(double* solution){
     4338
     4339        /*Intermediaries*/
     4340        const int numdof = NDOF1*NUMVERTICES;
     4341
     4342        int       i;
     4343        int*      doflist=NULL;
     4344        double    values[numdof];
     4345
     4346        /*Get dof list: */
     4347        GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     4348
     4349        /*Use the dof list to index into the solution vector: */
     4350        for(i=0;i<numdof;i++){
     4351                values[i]=solution[doflist[i]];
     4352                if(isnan(values[i])) _error_("NaN found in solution vector");
     4353                //if (values[i]<pow(10,-10))values[i]=pow(10,-10);
     4354        }
     4355
     4356        /*Now, we have to move the previous WaterColumn input  to Picard
     4357         * status, otherwise, we'll wipe them off: */
     4358        this->inputs->ChangeEnum(WaterColumnEnum,PicardWaterColumnEnum);
     4359
     4360        /*Add input to the element: */
     4361        this->inputs->AddInput(new TriaVertexInput(WaterColumnEnum,values));
     4362
     4363        /*Free ressources:*/
     4364        xfree((void**)&doflist);
     4365}
     4366/*}}}*/
    41024367/*FUNCTION Tria::InputUpdateFromVector(double* vector, int name, int type);{{{1*/
    41034368void  Tria::InputUpdateFromVector(double* vector, int name, int type){
    … …  
    42334498                                name==SurfaceSlopeYEnum ||
    42344499                                name==MeltingRateEnum ||
     4500                                name==WaterColumnEnum ||
    42354501                                name==AccumulationRateEnum ||
    42364502                                name==SurfaceAreaEnum||

  • issm/trunk/src/c/objects/Elements/Tria.h

    r7391 r7640  
    154154                ElementMatrix* CreateKMatrixDiagnosticVertSurface(void);
    155155                ElementMatrix* CreateKMatrixMelting(void);
     156                ElementMatrix* CreateKMatrixHydrology(void);
    156157                ElementMatrix* CreateKMatrixPrognostic(void);
    157158                ElementMatrix* CreateKMatrixPrognostic_CG(void);
    … …  
    169170                ElementVector* CreatePVectorAdjointBalancedthickness(void);
    170171                ElementVector* CreatePVectorDiagnosticHutter(void);
     172                ElementVector* CreatePVectorHydrology(void);
    171173                ElementVector* CreatePVectorPrognostic(void);
    172174                ElementVector* CreatePVectorPrognostic_CG(void);
    … …  
    185187                void      GetSolutionFromInputsDiagnosticHoriz(Vec solution);
    186188                void      GetSolutionFromInputsDiagnosticHutter(Vec solution);
     189                void      GetSolutionFromInputsHydrology(Vec solution);
    187190                void    GetStrainRate2d(double* epsilon,double* xyz_list, GaussTria* gauss, Input* vx_input, Input* vy_input);
    188191                void      GradjDragStokes(Vec gradient);
    … …  
    193196                void      InputUpdateFromSolutionOneDof(double* solution,int enum_type);
    194197                void      InputUpdateFromSolutionPrognostic(double* solution);
     198                void      InputUpdateFromSolutionHydrology(double* solution);
    195199                bool      IsInput(int name);
    196200                void      SetClone(int* minranks);
    197201                void      SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
     202                void      GetHydrologyK(double* K,double* xyz_list,GaussTria* gauss);
    198203                /*}}}*/
    199204

  • issm/trunk/src/c/objects/IoModel.cpp

    r7444 r7640  
    8080        xfree((void**)&this->spcvelocity);
    8181        xfree((void**)&this->spcthickness);
     82        xfree((void**)&this->spcwatercolumn);
    8283        xfree((void**)&this->spctemperature);
    8384        xfree((void**)&this->edges);
    8485        xfree((void**)&this->geothermalflux);
    8586        xfree((void**)&this->melting_rate);
     87        xfree((void**)&this->watercolumn);
    8688        xfree((void**)&this->melting_rate_correction);
    8789        xfree((void**)&this->accumulation_rate);
    … …  
    199201        /*!Get thermal parameters: */
    200202        IoModelFetchData(&this->beta,iomodel_handle,"beta");
     203        IoModelFetchData(&this->gamma,iomodel_handle,"gamma");
     204        IoModelFetchData(&this->kn,iomodel_handle,"kn");
    201205        IoModelFetchData(&this->meltingpoint,iomodel_handle,"meltingpoint");
    202206        IoModelFetchData(&this->latentheat,iomodel_handle,"latentheat");
    … …  
    274278        this->gl_melting_rate=0;
    275279        this->melting_rate=NULL;
     280        this->watercolumn=NULL;
    276281        this->melting_rate_correction=NULL;
    277282        this->melting_rate_correction_apply=0;
    … …  
    306311        this-> spctemperature=NULL;
    307312        this-> spcthickness=NULL;
     313        this-> spcwatercolumn=NULL;
    308314        this->numberofedges=0;
    309315        this->edges=NULL;

  • issm/trunk/src/c/objects/IoModel.h

    r7444 r7640  
    103103                double* spctemperature;
    104104                double* spcthickness;
     105                double* spcwatercolumn;
    105106                double* geothermalflux;
    106107                int     numberofedges;
    … …  
    163164                /*thermal parameters: */
    164165                double beta;
     166                double gamma;
     167                double kn;
    165168                double meltingpoint;
    166169                double latentheat;
    … …  
    182185                /*basal: */
    183186                double*  melting_rate;
     187                double*  watercolumn;
    184188                double   gl_melting_rate;
    185189                double*  melting_rate_correction;

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

    r5320 r7640  
    3535        double  matpar_thermal_exchange_velocity;
    3636        double  matpar_g;
     37        double  matpar_gamma;
     38        double  matpar_kn;
    3739
    3840        matpar_g=iomodel->g;
    … …  
    4648        matpar_mixed_layer_capacity=iomodel->mixed_layer_capacity;
    4749        matpar_thermal_exchange_velocity=iomodel->thermal_exchange_velocity;
     50        matpar_gamma=iomodel->gamma;
     51        matpar_kn=iomodel->kn;
    4852
    4953        this->mid=matpar_mid;
    … …  
    5862        this->thermal_exchange_velocity=matpar_thermal_exchange_velocity;
    5963        this->g=matpar_g;
     64        this->gamma=matpar_gamma;
     65        this->kn=matpar_kn;
    6066
    6167}
    … …  
    8389        printf("   thermal_exchange_velocity: %g\n",thermal_exchange_velocity);
    8490        printf("   g: %g\n",g);
     91        printf("   gamma: %g\n",gamma);
     92        printf("   kn: %g\n",kn);
    8593        return;
    8694}
    … …  
    101109        printf("   thermal_exchange_velocity: %g\n",thermal_exchange_velocity);
    102110        printf("   g: %g\n",g);
     111        printf("   gamma: %g\n",gamma);
     112        printf("   kn: %g\n",kn);
    103113        return;
    104114}               
    … …  
    140150        memcpy(marshalled_dataset,&thermal_exchange_velocity,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity);
    141151        memcpy(marshalled_dataset,&g,sizeof(g));marshalled_dataset+=sizeof(g);
     152        memcpy(marshalled_dataset,&gamma,sizeof(gamma));marshalled_dataset+=sizeof(gamma);
     153        memcpy(marshalled_dataset,&kn,sizeof(kn));marshalled_dataset+=sizeof(kn);
    142154
    143155        *pmarshalled_dataset=marshalled_dataset;
    … …  
    159171                sizeof(thermal_exchange_velocity)+
    160172                sizeof(g)+
     173                sizeof(gamma)+
     174                sizeof(kn)+
    161175                sizeof(int); //sizeof(int) for enum type
    162176}
    … …  
    184198        memcpy(&thermal_exchange_velocity,marshalled_dataset,sizeof(thermal_exchange_velocity));marshalled_dataset+=sizeof(thermal_exchange_velocity);
    185199        memcpy(&g,marshalled_dataset,sizeof(g));marshalled_dataset+=sizeof(g);
     200        memcpy(&gamma,marshalled_dataset,sizeof(gamma));marshalled_dataset+=sizeof(gamma);
     201        memcpy(&kn,marshalled_dataset,sizeof(kn));marshalled_dataset+=sizeof(kn);
    186202
    187203        /*return: */
    … …  
    359375}
    360376/*}}}1*/
    361 
    362 
     377/*FUNCTION Matpar::GetGamma {{{1*/
     378double Matpar::GetGamma(){
     379        return gamma;
     380}
     381/*}}}1*/
     382/*FUNCTION Matpar::GetKn {{{1*/
     383double Matpar::GetKn(){
     384        return kn;
     385}
     386/*}}}1*/

  • issm/trunk/src/c/objects/Materials/Matpar.h

    r6412 r7640  
    2626                double  thermal_exchange_velocity;
    2727                double  g;
     28                /*hydrology: */
     29                double  kn;
     30                double  gamma;
    2831
    2932        public:
    … …  
    7376                double GetBeta();
    7477                double GetMeltingPoint();
     78                double GetGamma();
     79                double GetKn();
    7580                /*}}}*/
    7681

  • issm/trunk/src/c/solutions/CorePointerFromSolutionEnum.cpp

    r7307 r7640  
    5858                        solutioncore=&groundinglinemigration2d_core;
    5959                        break;
    60 
     60                case HydrologySolutionEnum:
     61                        solutioncore=&hydrology_core;
     62                        break;
    6163                default:
    6264                        _error_("%s%s%s"," solution type: ",EnumToString(solutiontype)," not supported yet!");

  • issm/trunk/src/c/solutions/SolutionConfiguration.cpp

    r7307 r7640  
    5959                        analyses[0]=ThermalAnalysisEnum;
    6060                        analyses[1]=MeltingAnalysisEnum;
     61                        break;
     62               
     63                case HydrologySolutionEnum:
     64                        numanalyses=3;
     65                        analyses=(int*)xmalloc(numanalyses*sizeof(int));
     66                        analyses[0]=HydrologyAnalysisEnum;
     67                        analyses[1]=SurfaceSlopeAnalysisEnum;
     68                        analyses[2]=BedSlopeAnalysisEnum;
    6169                        break;
    6270

  • issm/trunk/src/c/solutions/solutions.h

    r7307 r7640  
    1717void gradient_core(FemModel* femmodel,int step=0, double search_scalar=0);
    1818void diagnostic_core(FemModel* femmodel);
     19void hydrology_core(FemModel* femmodel);
     20void hydrology_core_step(FemModel* femmodel,int step, double time);
    1921void thermal_core(FemModel* femmodel);
    2022void thermal_core_step(FemModel* femmodel,int step, double time);
Note: See TracChangeset for help on using the changeset viewer.