Changeset 9271


Ignore:
Timestamp:
08/11/11 08:29:42 (14 years ago)
Author:
Mathieu Morlighem
Message:

Implemented hydrological model with Anne LeBrocq's model

Location:
issm/trunk/src
Files:
3 added
15 edited

Legend:

Unmodified
Added
Removed

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

    r9218 r9271  
    484484        InputfilenameEnum,
    485485        SpcvzEnum,
    486         list,
     486        Fake33Enum,
    487487        NumberOfNodes2DEnum,
    488488        NodeOnStokesEnum,
    … …  
    548548        ViscousHeatingEnum,
    549549        QmuTemperatureEnum,
    550         QmuRheologyBEnum
     550        QmuRheologyBEnum,
     551        HydrologyWaterVxEnum,
     552        HydrologyWaterVyEnum
    551553};
    552554

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

    r9218 r9271  
    426426                case InputfilenameEnum : return "Inputfilename";
    427427                case SpcvzEnum : return "Spcvz";
     428                case Fake33Enum : return "Fake33";
    428429                case NumberOfNodes2DEnum : return "NumberOfNodes2D";
    429430                case NodeOnStokesEnum : return "NodeOnStokes";
    … …  
    490491                case QmuTemperatureEnum : return "QmuTemperature";
    491492                case QmuRheologyBEnum : return "QmuRheologyB";
     493                case HydrologyWaterVxEnum : return "HydrologyWaterVx";
     494                case HydrologyWaterVyEnum : return "HydrologyWaterVy";
    492495                default : return "unknown";
    493496

  • issm/trunk/src/c/modules/ModelProcessorx/Hydrology/CreateConstraintsHydrology.cpp

    r8926 r9271  
    44
    55#include "../../../Container/Container.h"
     6#include "../../../modules/modules.h"
    67#include "../../../toolkits/toolkits.h"
    78#include "../../../EnumDefinitions/EnumDefinitions.h"
    … …  
    1314void    CreateConstraintsHydrology(Constraints** pconstraints, IoModel* iomodel,FILE* iomodel_handle){
    1415
    15         /*Intermediary*/
    16         int i;
    17         int count;
    18        
    1916        /*Output*/
    2017        Constraints *constraints = NULL;
    21         Spc         *spc         = NULL;
    2218
    2319        /*Recover pointer: */
    … …  
    2622        /*Create constraints if they do not exist yet*/
    2723        if(!constraints) constraints = new Constraints(ConstraintsEnum);
    28 
    29         /*Fetch data: */
    30         IoModelFetchData(&iomodel->spcwatercolumn,NULL,NULL,iomodel_handle,SpcwatercolumnEnum);
    31 
    32         /*Initialize counter*/
    33         count=0;
    34 
    35         /*Create spcs from x,y,z, as well as the spc values on those spcs: */
    36         for (i=0;i<iomodel->numberofvertices;i++){
    37        
    38                 /*keep only this partition's nodes:*/
    39                 if((iomodel->my_vertices[i])){
    40 
    41                         if ((int)iomodel->spcwatercolumn[2*i]){
    42 
    43                                 constraints->AddObject(new Spc(iomodel->constraintcounter+count+1,iomodel->nodecounter+i+1,1,*(iomodel->spcwatercolumn+2*i+1),HydrologyAnalysisEnum));
    44                                 count++;
    45                         }
    46                 } //if((my_vertices[i]))
    47         }
    48 
    49         /*Free data: */
    50         xfree((void**)&iomodel->spcwatercolumn);
     24        IoModelToConstraintsx(constraints,iomodel,iomodel_handle,SpcwatercolumnEnum,HydrologyAnalysisEnum);
    5125       
    5226        /*Assign output pointer: */

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

    r9218 r9271  
    424424        else if (strcmp(name,"Inputfilename")==0) return InputfilenameEnum;
    425425        else if (strcmp(name,"Spcvz")==0) return SpcvzEnum;
     426        else if (strcmp(name,"Fake33")==0) return Fake33Enum;
    426427        else if (strcmp(name,"NumberOfNodes2D")==0) return NumberOfNodes2DEnum;
    427428        else if (strcmp(name,"NodeOnStokes")==0) return NodeOnStokesEnum;
    … …  
    488489        else if (strcmp(name,"QmuTemperature")==0) return QmuTemperatureEnum;
    489490        else if (strcmp(name,"QmuRheologyB")==0) return QmuRheologyBEnum;
     491        else if (strcmp(name,"HydrologyWaterVx")==0) return HydrologyWaterVxEnum;
     492        else if (strcmp(name,"HydrologyWaterVy")==0) return HydrologyWaterVyEnum;
    490493        else _error_("Enum %s not found",name);
    491494

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

    r9256 r9271  
    376376}
    377377/*}}}*/
     378/*FUNCTION Tria::CreateHydrologyWaterVelocityInput {{{1*/
     379void Tria::CreateHydrologyWaterVelocityInput(void){
     380
     381        /*material parameters: */
     382        double mu_water=0.001787;  //unit= [N s/m2]
     383        double w;
     384        double rho_ice, rho_water, g;
     385        double dsdx,dsdy,dbdx,dbdy;
     386        double vx[NUMVERTICES];
     387        double vy[NUMVERTICES];
     388        GaussTria *gauss = NULL;
     389
     390        /*Retrieve all inputs and parameters*/
     391        rho_ice=matpar->GetRhoIce();
     392        rho_water=matpar->GetRhoWater();
     393        g=matpar->GetG();
     394        Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
     395        Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
     396        Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum);         _assert_(bedslopex_input);
     397        Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum);         _assert_(bedslopey_input);
     398        Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);     _assert_(watercolumn_input);
     399
     400        gauss=new GaussTria();
     401        for (int iv=0;iv<NUMVERTICES;iv++){
     402                gauss->GaussVertex(iv);
     403                surfaceslopex_input->GetParameterValue(&dsdx,gauss);
     404                surfaceslopey_input->GetParameterValue(&dsdy,gauss);
     405                bedslopex_input->GetParameterValue(&dbdx,gauss);
     406                bedslopey_input->GetParameterValue(&dbdy,gauss);
     407                watercolumn_input->GetParameterValue(&w,gauss);
     408
     409                /* Water velocity x and y components */
     410                vx[iv]= pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
     411                vy[iv]= pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
     412                //vx[iv]=0.001;
     413                //vy[iv]=0;
     414        }
     415
     416        /*clean-up*/
     417        delete gauss;
     418
     419        /*Add to inputs*/
     420        this->inputs->AddInput(new TriaVertexInput(HydrologyWaterVxEnum,vx));
     421        this->inputs->AddInput(new TriaVertexInput(HydrologyWaterVyEnum,vy));
     422}
     423/*}}}*/
    378424/*FUNCTION Tria::CreateKMatrix {{{1*/
    379425void  Tria::CreateKMatrix(Mat Kff, Mat Kfs,Vec df){
    … …  
    783829}
    784830/*}}}*/
     831/*FUNCTION Tria::CreateKMatrixHydrology{{{1*/
     832ElementMatrix* Tria::CreateKMatrixHydrology(void){
     833
     834        /*Constants*/
     835        const int    numdof=NDOF1*NUMVERTICES;
     836
     837        /*Intermediaries */
     838        int        artdiff;
     839        int        i,j,ig;
     840        double     Jdettria,DL_scalar,dt,h;
     841        double     vx,vy,vel,dvxdx,dvydy;
     842        double     dvx[2],dvy[2];
     843        double     v_gauss[2]={0.0};
     844        double     xyz_list[NUMVERTICES][3];
     845        double     L[NUMVERTICES];
     846        double     B[2][NUMVERTICES];
     847        double     Bprime[2][NUMVERTICES];
     848        double     K[2][2]                        ={0.0};
     849        double     KDL[2][2]                      ={0.0};
     850        double     DL[2][2]                        ={0.0};
     851        double     DLprime[2][2]                   ={0.0};
     852        GaussTria *gauss=NULL;
     853
     854        /*Initialize Element matrix*/
     855        ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
     856
     857        /*Create water velocity vx and vy from current inputs*/
     858        CreateHydrologyWaterVelocityInput();
     859
     860        /*Retrieve all inputs and parameters*/
     861        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
     862        this->parameters->FindParam(&dt,DtEnum);
     863        this->parameters->FindParam(&artdiff,ArtDiffEnum);
     864        Input* vx_input=inputs->GetInput(HydrologyWaterVxEnum); _assert_(vx_input);
     865        Input* vy_input=inputs->GetInput(HydrologyWaterVyEnum); _assert_(vy_input);
     866        h=this->GetArea();
     867
     868        /* Start  looping on the number of gaussian points: */
     869        gauss=new GaussTria(2);
     870        for (ig=gauss->begin();ig<gauss->end();ig++){
     871
     872                gauss->GaussPoint(ig);
     873
     874                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
     875                GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
     876
     877                vx_input->GetParameterValue(&vx,gauss);
     878                vy_input->GetParameterValue(&vy,gauss);
     879                vx_input->GetParameterDerivativeValue(&dvx[0],&xyz_list[0][0],gauss);
     880                vy_input->GetParameterDerivativeValue(&dvy[0],&xyz_list[0][0],gauss);
     881
     882                DL_scalar=gauss->weight*Jdettria;
     883
     884                TripleMultiply( &L[0],1,numdof,1,
     885                                        &DL_scalar,1,1,0,
     886                                        &L[0],1,numdof,0,
     887                                        &Ke->values[0],1);
     888
     889                GetBPrognostic(&B[0][0], &xyz_list[0][0], gauss);
     890                GetBprimePrognostic(&Bprime[0][0], &xyz_list[0][0], gauss);
     891
     892                dvxdx=dvx[0];
     893                dvydy=dvy[1];
     894                DL_scalar=dt*gauss->weight*Jdettria;
     895
     896                DL[0][0]=DL_scalar*dvxdx;
     897                DL[1][1]=DL_scalar*dvydy;
     898                DLprime[0][0]=DL_scalar*vx;
     899                DLprime[1][1]=DL_scalar*vy;
     900
     901                TripleMultiply( &B[0][0],2,numdof,1,
     902                                        &DL[0][0],2,2,0,
     903                                        &B[0][0],2,numdof,0,
     904                                        &Ke->values[0],1);
     905
     906                TripleMultiply( &B[0][0],2,numdof,1,
     907                                        &DLprime[0][0],2,2,0,
     908                                        &Bprime[0][0],2,numdof,0,
     909                                        &Ke->values[0],1);
     910
     911                if(artdiff){
     912                        vel=sqrt(pow(vx,2.)+pow(vy,2.));
     913                        K[0][0]=h/(2*vel)*vx*vx;
     914                        K[1][0]=h/(2*vel)*vy*vx;
     915                        K[0][1]=h/(2*vel)*vx*vy;
     916                        K[1][1]=h/(2*vel)*vy*vy;
     917                        KDL[0][0]=DL_scalar*K[0][0];
     918                        KDL[1][0]=DL_scalar*K[1][0];
     919                        KDL[0][1]=DL_scalar*K[0][1];
     920                        KDL[1][1]=DL_scalar*K[1][1];
     921
     922                        TripleMultiply( &Bprime[0][0],2,numdof,1,
     923                                                &KDL[0][0],2,2,0,
     924                                                &Bprime[0][0],2,numdof,0,
     925                                                &Ke->values[0],1);
     926                }
     927        }
     928
     929        /*Clean up and return*/
     930        delete gauss;
     931        return Ke;
     932}
     933/*}}}*/
    785934/*FUNCTION Tria::CreateKMatrixMelting {{{1*/
    786935ElementMatrix* Tria::CreateKMatrixMelting(void){
    … …  
    820969                                        &L[0],1,numdof,0,
    821970                                        &Ke->values[0],1);
    822         }
    823 
    824         /*Clean up and return*/
    825         delete gauss;
    826         return Ke;
    827 }
    828 /*}}}*/
    829 /*FUNCTION Tria::CreateKMatrixHydrology {{{1*/
    830 ElementMatrix* Tria::CreateKMatrixHydrology(void){
    831 
    832         /*Constants*/
    833         const int    numdof=NDOF1*NUMVERTICES;
    834 
    835         /*Intermediaries */
    836         int        i,j,ig;
    837         double     Jdettria,dt;
    838         double     xyz_list[NUMVERTICES][3];
    839         double     L[NUMVERTICES];
    840         double     dL[2][NUMVERTICES];
    841         double     dLx[NUMVERTICES];
    842         double     dLy[NUMVERTICES];
    843         double     K[2];
    844         double     hydro_p,hydro_q,hydro_CR,hydro_n,rho_water,g,gamma;
    845         double     DL_scalar1;
    846         double     DL_scalar2;
    847         double     DL_scalar3;
    848         GaussTria *gauss=NULL;
    849 
    850         /*Initialize Element matrix*/
    851         ElementMatrix* Ke=new ElementMatrix(nodes,NUMVERTICES,this->parameters,NoneApproximationEnum);
    852 
    853         /*Retrieve all inputs and parameters*/
    854         GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    855         this->parameters->FindParam(&dt,DtEnum);
    856 
    857         /*retrieve material parameters: */
    858         g        =matpar->GetG();
    859         rho_water=matpar->GetRhoWater();
    860         hydro_p  =matpar->GetHydroP();
    861         hydro_q  =matpar->GetHydroQ();
    862         hydro_CR =matpar->GetHydroCR();
    863         hydro_n  =matpar->GetHydroN();
    864         gamma=1/(hydro_n*pow(hydro_CR,hydro_p)*pow((rho_water*g),hydro_q));
    865        
    866         /* Start  looping on the number of gaussian points: */
    867         gauss=new GaussTria(2);
    868         for (ig=gauss->begin();ig<gauss->end();ig++){
    869 
    870                 gauss->GaussPoint(ig);
    871 
    872                 GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    873                 GetNodalFunctions(&L[0],gauss);
    874                 GetNodalFunctionsDerivatives(&dL[0][0],&xyz_list[0][0],gauss);
    875                 for(i=0;i<NUMVERTICES;i++)dLx[i]=dL[0][i];
    876                 for(i=0;i<NUMVERTICES;i++)dLy[i]=dL[1][i];
    877 
    878                 /*Get K parameter: */
    879                 GetHydrologyK(&K[0],&xyz_list[0][0],gauss);
    880 
    881                 if(dt){
    882                         DL_scalar1=gauss->weight*Jdettria;
    883                         TripleMultiply( &L[0],1,numdof,1,
    884                                         &DL_scalar1,1,1,0,
    885                                         &L[0],1,numdof,0,
    886                                         &Ke->values[0],1);
    887                 }
    888 
    889                 if(dt){
    890                         DL_scalar2=-gauss->weight*Jdettria*gamma*K[0]*dt;
    891                         DL_scalar3=-gauss->weight*Jdettria*gamma*K[1]*dt;
    892 
    893                 }
    894                 else{
    895                         DL_scalar2=-gauss->weight*Jdettria*gamma*K[0];
    896                         DL_scalar3=-gauss->weight*Jdettria*gamma*K[1];
    897                 }
    898                
    899                 TripleMultiply(&dLx[0],1,numdof,1, &DL_scalar2,1,1,0, &L[0],1,numdof,0,&Ke->values[0],1);
    900                 TripleMultiply(&dLy[0],1,numdof,1, &DL_scalar3,1,1,0, &L[0],1,numdof,0,&Ke->values[0],1);
    901971        }
    902972
    … …  
    18651935        double     old_watercolumn_g;
    18661936        double     xyz_list[NUMVERTICES][3];
    1867         double     L[NUMVERTICES];
     1937        double     basis[numdof];
    18681938        GaussTria* gauss=NULL;
    18691939
    … …  
    18741944        GetVerticesCoordinates(&xyz_list[0][0], nodes, NUMVERTICES);
    18751945        this->parameters->FindParam(&dt,DtEnum);
    1876         Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum);           _assert_(basal_melting_input);
     1946        Input* basal_melting_input=inputs->GetInput(BasalMeltingRateEnum);_assert_(basal_melting_input);
    18771947        Input* old_watercolumn_input=inputs->GetInput(WaterColumnOldEnum);_assert_(old_watercolumn_input);
    18781948
    … …  
    18851955
    18861956                GetJacobianDeterminant2d(&Jdettria, &xyz_list[0][0],gauss);
    1887                 GetL(&L[0], &xyz_list[0][0], gauss,NDOF1);
     1957                GetNodalFunctions(basis, gauss);
    18881958
    18891959                basal_melting_input->GetParameterValue(&basal_melting_g,gauss);
    18901960                old_watercolumn_input->GetParameterValue(&old_watercolumn_g,gauss);
    1891                 //printf("ig=%d, old_watercolumn_g=%f,Jdettria=%f\n",ig,old_watercolumn_g,Jdettria);//bk
    1892                 //if(old_watercolumn_g<0) printf("old_watercolumn_g=%g,Jdettria=%f\n",old_watercolumn_g,Jdettria);//bk
    1893 
    1894                 if(dt)for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*L[i];
    1895                 else  for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*L[i];
     1961
     1962                if(dt)for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*(old_watercolumn_g+dt*basal_melting_g)*basis[i];
     1963                else  for(i=0;i<numdof;i++) pe->values[i]+=Jdettria*gauss->weight*basal_melting_g*basis[i];
    18961964        }
    18971965               
    … …  
    23342402        return &this->horizontalneighborsids[0];
    23352403
    2336 }
    2337 /*}}}*/
    2338 /*FUNCTION Tria::GetHydrologyK {{{1*/
    2339 void Tria::GetHydrologyK(double* K,double* xyz_list,GaussTria* gauss){
    2340 
    2341         /*material parameters: */
    2342         double rho_ice, rho_water, g,hydro_p,hydro_q;
    2343         double dsdx, dsdy;
    2344         double dbdx, dbdy;
    2345         double surface_slope;
    2346         double kn, w;
    2347        
    2348         /*Retrieve all inputs and parameters*/
    2349         rho_ice=matpar->GetRhoIce();
    2350         rho_water=matpar->GetRhoWater();
    2351         g=matpar->GetG();
    2352         kn=matpar->GetKn();
    2353         hydro_p=matpar->GetHydroP();
    2354         hydro_q=matpar->GetHydroQ();
    2355         Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
    2356         Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
    2357         Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum);         _assert_(bedslopex_input);
    2358         Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum);         _assert_(bedslopey_input);
    2359         Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);     _assert_(watercolumn_input);
    2360 
    2361         /*recover slopes: */
    2362         surfaceslopex_input->GetParameterValue(&dsdx,gauss);
    2363         surfaceslopey_input->GetParameterValue(&dsdy,gauss);
    2364         bedslopex_input->GetParameterValue(&dbdx,gauss);
    2365         bedslopey_input->GetParameterValue(&dbdy,gauss);
    2366 
    2367         /*magnitude of surface slope: */
    2368         surface_slope=sqrt( pow(dsdx,2) + pow(dsdy,2));
    2369 
    2370         /*recover water column: */
    2371         watercolumn_input->GetParameterValue(&w,gauss);
    2372 
    2373         K[0]=pow(w,hydro_p)*pow((rho_ice*g*dsdx+(rho_water/rho_ice-1)*rho_ice*g*dbdx) - rho_ice * g * kn/w * (dsdx - dbdx ) * surface_slope,hydro_q);
    2374         K[1]=pow(w,hydro_p)*pow((rho_ice*g*dsdy+(rho_water/rho_ice-1)*rho_ice*g*dbdy) - rho_ice * g * kn/w * (dsdy - dbdy ) * surface_slope,hydro_q);
    2375        
    2376 }
    2377 /*}}}*/
    2378 /*FUNCTION Tria::GetHydrologyWaterVelocity {{{1*/
    2379 void Tria::GetHydrologyWaterVelocity(double* v,GaussTria* gauss){
    2380         /*material parameters: */
    2381         double mu_water=0.001787;  //unit= [N s/m2]
    2382         double w;
    2383         double rho_ice, rho_water, g;
    2384         double dsdx,dsdy,dbdx,dbdy;
    2385 
    2386         /*Retrieve all inputs and parameters*/
    2387         rho_ice=matpar->GetRhoIce();
    2388         rho_water=matpar->GetRhoWater();
    2389         g=matpar->GetG();
    2390         Input* surfaceslopex_input=inputs->GetInput(SurfaceSlopeXEnum); _assert_(surfaceslopex_input);
    2391         Input* surfaceslopey_input=inputs->GetInput(SurfaceSlopeYEnum); _assert_(surfaceslopey_input);
    2392         Input* bedslopex_input=inputs->GetInput(BedSlopeXEnum);         _assert_(bedslopex_input);
    2393         Input* bedslopey_input=inputs->GetInput(BedSlopeYEnum);         _assert_(bedslopey_input);
    2394         Input* watercolumn_input=inputs->GetInput(WatercolumnEnum);     _assert_(watercolumn_input);
    2395 
    2396         /*recover slopes: */
    2397         surfaceslopex_input->GetParameterValue(&dsdx,gauss);
    2398         surfaceslopey_input->GetParameterValue(&dsdy,gauss);
    2399         bedslopex_input->GetParameterValue(&dbdx,gauss);
    2400         bedslopey_input->GetParameterValue(&dbdy,gauss);
    2401 
    2402         /*recover water column: */
    2403         watercolumn_input->GetParameterValue(&w,gauss);
    2404 
    2405         /* Water velocity x and y components */
    2406         v[0]=pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdx+(rho_water-rho_ice)*g*dbdx);
    2407         v[1]=pow(w,2)/(12 * mu_water)*(rho_ice*g*dsdy+(rho_water-rho_ice)*g*dbdy);
    24082404}
    24092405/*}}}*/

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

    r9256 r9271  
    181181                int     GetElementType(void);
    182182                void      GetDofList(int** pdoflist,int approximation_enum,int setenum);
    183                 void    GetHydrologyWaterVelocity(double* v, GaussTria* gauss);
     183                void    CreateHydrologyWaterVelocityInput(void);
    184184                void      GetDofList1(int* doflist);
    185185                void    GetSidList(int* sidlist);
    … …  
    202202                void      SetClone(int* minranks);
    203203                void      SurfaceNormal(double* surface_normal, double xyz_list[3][3]);
    204                 void      GetHydrologyK(double* K,double* xyz_list,GaussTria* gauss);
    205204                /*}}}*/
    206205

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

    r9013 r9271  
    223223        IoModelFetchData(&this->gl_melting_rate,iomodel_handle,GlMeltingRateEnum);
    224224        IoModelFetchData(&this->rheology_law,iomodel_handle,RheologyLawEnum);
    225        
    226         /*!Get hydrology parameters: */
    227         IoModelFetchData(&this->hydro_kn,iomodel_handle,HydroKnEnum);
    228         IoModelFetchData(&this->hydro_p,iomodel_handle,HydroPEnum);
    229         IoModelFetchData(&this->hydro_q,iomodel_handle,HydroQEnum);
    230         IoModelFetchData(&this->hydro_CR,iomodel_handle,HydroCREnum);
    231         IoModelFetchData(&this->hydro_n,iomodel_handle,HydroNEnum);
    232225
    233226        /*qmu: */

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

    r9013 r9271  
    215215                int     loadcounter; //keep track of how many loads are being created in each analysis type
    216216                int     constraintcounter; //keep track of how many constraints are being created in each analysis type
    217 
    218                 /*hydrology parameters: */
    219                 double hydro_kn;
    220                 double hydro_p;
    221                 double hydro_q;
    222                 double hydro_CR;
    223                 double hydro_n;
    224217                 /*}}}*/
    225218                /*Methods: {{{1*/

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

    r8926 r9271  
    3737        this->thermal_exchange_velocity = iomodel->thermal_exchange_velocity;
    3838        this->g                         = iomodel->g;
    39         this->kn                        = iomodel->hydro_kn;
    40         this->hydro_p                   = iomodel->hydro_p;
    41         this->hydro_q                   = iomodel->hydro_q;
    42         this->hydro_CR                  = iomodel->hydro_CR;
    43         this->hydro_n                   = iomodel->hydro_n;
    4439}
    4540/*}}}1*/
    … …  
    361356}
    362357/*}}}1*/
    363 /*FUNCTION Matpar::GetKn {{{1*/
    364 double Matpar::GetKn(){
    365         return kn;
    366 }
    367 /*}}}1*/
    368 /*FUNCTION Matpar::GetHydroP {{{1*/
    369 double Matpar::GetHydroP(){
    370         return hydro_p;
    371 }
    372 /*}}}1*/
    373 /*FUNCTION Matqar::GetHydroQ {{{1*/
    374 double Matpar::GetHydroQ(){
    375         return hydro_q;
    376 }
    377 /*}}}1*/
    378 /*FUNCTION Matpar::GetHydroCR {{{1*/
    379 double Matpar::GetHydroCR(){
    380            return hydro_CR;
    381 }
    382 /*}}}1*/
    383 /*FUNCTION Matpar::GetHydroN {{{1*/
    384 double Matpar::GetHydroN(){
    385            return hydro_n;
    386 }
    387 /*}}}1*/
    388358/*FUNCTION Matpar::TMeltingPoint {{{1*/
    389359double Matpar::TMeltingPoint(double pressure){

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

    r8561 r9271  
    2727                double  thermal_exchange_velocity;
    2828                double  g;
    29 
    30                 /*hydrology: */
    31                 double  kn;
    32                 double  hydro_p;
    33                 double  hydro_q;
    34                 double  hydro_CR;
    35                 double  hydro_n;
    3629
    3730        public:

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

    r8926 r9271  
    4545                if(nsteps)_printf_(VerboseSolution(),"time step:%i/%i\n",i+1,nsteps);
    4646                time=(i+1)*dt;
     47                femmodel->parameters->SetParam(time,TimeEnum);
    4748
    4849                /*call hydrology_core step: */
    … …  
    5253                        _printf_(VerboseSolution(),"   saving results\n");
    5354                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,WatercolumnEnum,i+1,time);
     55                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,HydrologyWaterVxEnum,i+1,time);
     56                        InputToResultx(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,HydrologyWaterVyEnum,i+1,time);
    5457                }
    5558

  • issm/trunk/src/c/solvers/solver_adjoint_linear.cpp

    r8821 r9271  
    2222        /*Recover parameters: */
    2323        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     24        UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
    2425
    2526        SystemMatricesx(&Kff, &Kfs, &pf, &df, NULL,femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters);

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

    r8821 r9271  
    3131        int  configuration_type;
    3232
    33 
    3433        /*Recover parameters: */
    3534        femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum);
    3635        femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum);
    3736        femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
     37        UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
    3838
    3939        /*Were loads requested as output? : */

  • issm/trunk/src/c/solvers/solver_stokescoupling_nonlinear.cpp

    r8834 r9271  
    3535        femmodel->parameters->FindParam(&min_mechanical_constraints,MinMechanicalConstraintsEnum);
    3636        femmodel->parameters->FindParam(&max_nonlinear_iterations,MaxNonlinearIterationsEnum);
     37        UpdateConstraintsx(femmodel->nodes,femmodel->constraints,femmodel->parameters);
    3738       
    3839        count=1;

  • issm/trunk/src/m/classes/model.m

    r9267 r9271  
    268268                 basal_melting_rate_correction_apply = {0,true,'Integer'};
    269269                 pressure                            = {NaN,true,'DoubleMat',1};
    270                  %Hydrology
    271270                 watercolumn                         = {NaN,true,'DoubleMat',1};
    272                  hydro_n                             = {0,true,'Double'};
    273                  hydro_CR                            = {0,true,'Double'};
    274                  hydro_p                             = {0,true,'Double'};
    275                  hydro_q                             = {0,true,'Double'};
    276                  hydro_kn                            = {0,true,'Double'};
    277271                 
    278272                 %Parallelisation
    … …  
    849843                         md.petscoptions=addoptions(md.petscoptions,DiagnosticVertAnalysisEnum,mumpsoptions);
    850844
    851 
    852                          %hydrology:  from Johnson's 2002 thesis, section 3.5.4
    853                          md.hydro_n=.02;
    854                          md.hydro_CR=2;
    855                          md.hydro_p=2;
    856                          md.hydro_q=1;
    857                          md.hydro_kn=0;
    858 
    859845                         %Rheology law: what is the temperature dependence of B with T
    860846                         %available: None, Paterson and Arrhenius
Note: See TracChangeset for help on using the changeset viewer.