Changeset 20608


Ignore:
Timestamp:
05/11/16 20:06:40 (9 years ago)
Author:
Mathieu Morlighem
Message:

CHG: moving viscosity to material in preparation for new flow law

Location:
issm/trunk-jpl/src/c
Files:
11 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk-jpl/src/c/analyses/AdjointHorizAnalysis.cpp

    r19161 r20608  
    105105
    106106                element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    107                 element->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]);
     107                element->material->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    108108                eps1[0]=epsilon[0];   eps2[0]=epsilon[2];   eps3[0]=epsilon[3];
    109109                eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
    … …  
    179179
    180180                element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    181                 element->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]);
     181                element->material->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]);
    182182                eps1[0]=2.*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    183183                eps1[1]=epsilon[2];                 eps2[1]=epsilon[0]+2.*epsilon[1];
    … …  
    283283                thickness_input->GetInputValue(&thickness, gauss);
    284284                basalelement->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    285                 basalelement->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]);
     285                basalelement->material->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]);
    286286                eps1[0]=2.*epsilon[0]+epsilon[1]; eps2[0]=epsilon[2];
    287287                eps1[1]=epsilon[2];               eps2[1]=epsilon[0]+2*epsilon[1];

  • issm/trunk-jpl/src/c/analyses/StressbalanceAnalysis.cpp

    r20459 r20608  
    12001200                thickness_input->GetInputValue(&thickness, gauss);
    12011201                basalelement->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    1202                 basalelement->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]);
     1202                basalelement->material->ViscositySSADerivativeEpsSquare(&mu_prime,&epsilon[0]);
    12031203                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    12041204                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
    … …  
    14741474                this->GetBSSAprime(Bprime,element,dim,xyz_list,gauss);
    14751475
    1476                 element->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
    1477                 element->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
     1476                element->material->ViscositySSA(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
     1477                element->material->ViscositySSA(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
    14781478                thickness_input->GetInputValue(&thickness, gauss);
    14791479
    … …  
    19611961                this->GetBSSAprime(Bprime,basalelement,2,xyz_list,gauss_base);
    19621962
    1963                 element->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
     1963                element->material->ViscosityL1L2(&viscosity,xyz_list,gauss,vx_input,vy_input,surface_input);
    19641964
    19651965                for(int i=0;i<3;i++) D[i*3+i]=2*viscosity*gauss->weight*Jdet;
    … …  
    22572257
    22582258                element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    2259                 element->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]);
     2259                element->material->ViscosityHODerivativeEpsSquare(&mu_prime,&epsilon[0]);
    22602260                eps1[0]=2*epsilon[0]+epsilon[1];   eps2[0]=epsilon[2];
    22612261                eps1[1]=epsilon[2];                eps2[1]=epsilon[0]+2*epsilon[1];
    … …  
    24282428                this->GetBHOprime(Bprime,element,dim,xyz_list,gauss);
    24292429
    2430                 element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
    2431                 element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
     2430                element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
     2431                element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
    24322432
    24332433                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    … …  
    29372937                eps1[1]=epsilon[2];   eps2[1]=epsilon[1];   eps3[1]=epsilon[4];
    29382938                eps1[2]=epsilon[3];   eps2[2]=epsilon[4];   eps3[2]= -epsilon[0] -epsilon[1];
    2939                 element->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]);
     2939                element->material->ViscosityFSDerivativeEpsSquare(&mu_prime,&epsilon[0]);
    29402940
    29412941                for(i=0;i<vnumnodes;i++){
    … …  
    31163116                this->GetBFSprime(Bprime,element,dim,xyz_list,gauss);
    31173117
    3118                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     3118                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    31193119                for(i=0;i<epssize;i++)     D[i*bsize+i] = + 2.*viscosity*gauss->weight*Jdet;
    31203120                for(i=epssize;i<bsize;i++) D[i*bsize+i] = - FSreconditioning*gauss->weight*Jdet;
    … …  
    31913191                this->GetBFSprimevel(Bprime,element,dim,xyz_list,gauss);
    31923192
    3193                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     3193                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    31943194                for(i=0;i<epssize;i++)   D[i*epssize+i] = 2*viscosity*gauss->weight*Jdet;
    31953195
    … …  
    48814881                        }
    48824882
    4883                         element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     4883                        element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    48844884                        epsxx[i] = dvx[0];                sigmapxx[i] = 2.*viscosity*epsxx[i];
    48854885                        epsyy[i] = dvy[1];                sigmapyy[i] = 2.*viscosity*epsyy[i];
    … …  
    55085508                this->GetLprimeFSSSA(&LprimeFSSSA[0][0], element,xyz_list, gauss);
    55095509
    5510                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     5510                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    55115511
    55125512                element->NormalBase(&bed_normal[0],xyz_list_tria);
    … …  
    56255625                this->GetBprimeSSAFS(&Bprime2[0][0], element,xyz_list, gauss);
    56265626
    5627                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     5627                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    56285628
    56295629                D_scalar=2*viscosity*gauss->weight*Jdet;
    … …  
    58055805                this->GetBSSAHO(B, element,xyz_list, gauss);
    58065806                this->GetBSSAprime(Bprime,basaltria,2,xyz_list, gauss_tria);
    5807                 element->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input);
    5808                 element->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input);
     5807                element->material->ViscosityHO(&viscosity,3,xyz_list,gauss,vx_input,vy_input);
     5808                element->material->ViscosityHO(&oldviscosity,3,xyz_list,gauss,vxold_input,vyold_input);
    58095809
    58105810                newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    … …  
    59645964
    59655965                if(approximation==SSAHOApproximationEnum){
    5966                         element->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
    5967                         element->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
     5966                        element->material->ViscosityHO(&viscosity,dim,xyz_list,gauss,vx_input,vy_input);
     5967                        element->material->ViscosityHO(&oldviscosity,dim,xyz_list,gauss,vxold_input,vyold_input);
    59685968                        newviscosity=viscosity+viscosity_overshoot*(viscosity-oldviscosity);
    59695969                }
    59705970                else if (approximation==SSAFSApproximationEnum){
    5971                         element->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     5971                        element->material->ViscosityFS(&newviscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    59725972                }
    59735973                else _error_("approximation " << approximation << " (" << EnumToStringx(approximation) << ") not supported yet");
    … …  
    60696069
    60706070                element->NormalBase(&bed_normal[0],xyz_list_tria);
    6071                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     6071                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    60726072                friction->GetAlpha2(&alpha2_gauss,gauss);
    60736073
    … …  
    61396139                element->NodalFunctionsP1Derivatives(&dbasis[0][0],xyz_list,gauss);
    61406140               
    6141                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     6141                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    61426142                vzHO_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
    61436143
    … …  
    62306230
    62316231                element->NormalBase(&bed_normal[0],xyz_list_tria);
    6232                 element->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
     6232                element->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,vz_input);
    62336233                friction->GetAlpha2(&alpha2_gauss,gauss);
    62346234
    … …  
    62996299
    63006300                vzSSA_input->GetInputDerivativeValue(&dw[0],xyz_list,gauss);
    6301                 element->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
     6301                element->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
    63026302
    63036303                for(i=0;i<6;i++){

  • issm/trunk-jpl/src/c/classes/Elements/Element.cpp

    r20603 r20608  
    28792879        xDelete<IssmDouble>(values);
    28802880}/*}}}*/
    2881 void       Element::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
    2882         /*The effective strain rate is defined in Paterson 3d Ed p 91 eq 9,
    2883          * and Cuffey p 303 eq 8.18:
    2884          *
    2885          *  2 eps_eff^2 = eps_xx^2 + eps_yy^2 + eps_zz^2 + 2(eps_xy^2 + eps_xz^2 + eps_yz^2)
    2886          *
    2887          *  or
    2888          *
    2889          *  eps_eff = 1/sqrt(2) sqrt( \sum_ij eps_ij^2 )
    2890          *
    2891          *          = 1/sqrt(2) ||eps||_F
    2892          *
    2893          *  where ||.||_F is the Frobenius norm */
    2894 
    2895         /*Intermediaries*/
    2896         IssmDouble viscosity;
    2897         IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
    2898         IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];            */
    2899         IssmDouble eps_eff;
    2900         IssmDouble eps0=1.e-27;
    2901 
    2902         if(dim==3){
    2903                 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
    2904                 this->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    2905                 eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
    2906         }
    2907         else{
    2908                 /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
    2909                 this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
    2910                 eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
    2911         }
    2912 
    2913         /*Get viscosity*/
    2914         material->GetViscosity(&viscosity,eps_eff);
    2915 
    2916         /*Assign output pointer*/
    2917         *pviscosity=viscosity;
    2918 }
    2919 /*}}}*/
    2920 void       Element::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
    2921         this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
    2922 }/*}}}*/
    2923 void       Element::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
    2924 
    2925         /*Intermediaries*/
    2926         IssmDouble viscosity;
    2927         IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
    2928         IssmDouble epsilon2d[2];/* epsilon=[exx,exy];            */
    2929         IssmDouble eps_eff;
    2930 
    2931         if(dim==3){
    2932                 /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
    2933                 this->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
    2934                 eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] +  epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]);
    2935         }
    2936         else{
    2937                 /* eps_eff^2 = 1/2 (exx^2 + 2*exy^2 )*/
    2938                 this->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
    2939                 eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1]);
    2940         }
    2941 
    2942         /*Get viscosity*/
    2943         material->GetViscosity(&viscosity,eps_eff);
    2944 
    2945         /*Assign output pointer*/
    2946         *pviscosity=viscosity;
    2947 }/*}}}*/
    2948 void       Element::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
    2949         this->material->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
    2950 }/*}}}*/
    2951 void       Element::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
    2952         /*Compute the L1L2 viscosity
    2953          *
    2954          *      1
    2955          * mu = - A^-1 (sigma'_e)^(1-n)
    2956          *      2
    2957          *
    2958          * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18)
    2959          *
    2960          * L1L2 assumptions:
    2961          *
    2962          * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//|
    2963          * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2
    2964          *
    2965          * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|)
    2966          *
    2967          * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0     */
    2968 
    2969         IssmDouble z,s,viscosity,p,q,delta;
    2970         IssmDouble tau_perp,tau_par,eps_b,A;
    2971         IssmDouble epsilon[5];   /*exx eyy exy exz eyz*/
    2972         IssmDouble slope[3];
    2973 
    2974         /*Check that both inputs have been found*/
    2975         if (!vx_input || !vy_input || !surface_input) _error_("Input missing");
    2976 
    2977         /*Get tau_perp*/
    2978         surface_input->GetInputValue(&s,gauss);
    2979         surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
    2980         z=this->GetZcoord(xyz_list,gauss);
    2981         tau_perp = matpar->GetMaterialParameter(MaterialsRhoIceEnum) * matpar->GetMaterialParameter(ConstantsGEnum) * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
    2982 
    2983         /* Get eps_b*/
    2984         this->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    2985         eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
    2986         if(eps_b==0.){
    2987                 *pviscosity = 2.5e+17;
    2988                 return;
    2989         }
    2990 
    2991         /*Get A*/
    2992         _assert_(material->GetN()==3.0);
    2993         A=material->GetA();
    2994 
    2995         /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/
    2996         p     = tau_perp *tau_perp;
    2997         q     = - eps_b/A;
    2998         delta = q *q + p*p*p*4./27.;
    2999         _assert_(delta>0);
    3000         tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.);
    3001 
    3002         /*Viscosity*/
    3003         viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp));
    3004         _assert_(!xIsNan(viscosity));
    3005         _assert_(viscosity > 0.);
    3006 
    3007         /*Assign output pointer*/
    3008         *pviscosity = viscosity;
    3009 }/*}}}*/
    3010 void       Element::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
    3011 
    3012         /*Intermediaries*/
    3013         IssmDouble viscosity;
    3014         IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
    3015         IssmDouble epsilon1d;   /* epsilon=[exx];    */
    3016         IssmDouble eps_eff;
    3017 
    3018          if(dim==2){
    3019                  /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
    3020                  this->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
    3021                  eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
    3022          }
    3023          else{
    3024                  /* eps_eff^2 = 1/2 exx^2*/
    3025                  this->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
    3026                  eps_eff = sqrt(epsilon1d*epsilon1d/2.);
    3027          }
    3028 
    3029         /*Get viscosity*/
    3030         material->GetViscosityBar(&viscosity,eps_eff);
    3031 
    3032         /*Assign output pointer*/
    3033         *pviscosity=viscosity;
    3034 }/*}}}*/
    3035 void       Element::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
    3036         this->material->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon);
    3037 }/*}}}*/
    30382881void       Element::ViscousHeatingCreateInput(void){/*{{{*/
    30392882
    … …  
    30642907
    30652908                this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    3066                 this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
     2909                this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
    30672910                this->GetPhi(&phi,&epsilon[0],viscosity);
    30682911

  • issm/trunk-jpl/src/c/classes/Elements/Element.h

    r20461 r20608  
    167167                void               TransformStiffnessMatrixCoord(ElementMatrix* Ke,Node** nodes,int numnodes,int* cs_array);
    168168                void               TransformStiffnessMatrixCoord(ElementMatrix* Ke,int numnodes,int* transformenum_list){_error_("not implemented yet");};/*Tiling only*/
    169                 void               ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
    170                 void               ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
    171                 void               ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    172                 void               ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
    173                 void               ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
    174                 void               ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
    175                 void               ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
    176169                void               ViscousHeatingCreateInput(void);
    177170

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

    r20461 r20608  
    289289                /*Compute strain rate viscosity and pressure: */
    290290                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    291                 this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     291                this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    292292                pressure_input->GetInputValue(&pressure,gauss);
    293293
    … …  
    346346                /*Compute strain rate viscosity and pressure: */
    347347                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    348                 this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     348                this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    349349
    350350                /*Compute Stress*/
    … …  
    398398                /*Compute strain rate viscosity and pressure: */
    399399                this->StrainRateFS(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    400                 this->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
     400                this->material->ViscosityFS(&viscosity,3,&xyz_list[0][0],gauss,vx_input,vy_input,vz_input);
    401401                pressure_input->GetInputValue(&pressure,gauss);
    402402
    … …  
    640640                        /*Check for basal force only if grounded and touching GL*/
    641641                        this->StrainRateFS(&epsilon[0],xyz_list,gauss,vx_input,vy_input,vz_input);
    642                         this->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
     642                        this->material->ViscosityFS(&viscosity,3,xyz_list,gauss,vx_input,vy_input,vz_input);
    643643                        pressure_input->GetInputValue(&pressure, gauss);
    644644                        base_input->GetInputValue(&base, gauss); _assert_(base<0.);
    … …  
    31743174        _assert_(gauss->Enum()==GaussPentaEnum);
    31753175        this->StrainRateFS(&epsilon[0],xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
    3176         this->ViscosityFS(&viscosity,3,xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
     3176        this->material->ViscosityFS(&viscosity,3,xyz_list,(GaussPenta*)gauss,vx_input,vy_input,vz_input);
    31773177        GetPhi(&phi,&epsilon[0],viscosity);
    31783178

  • issm/trunk-jpl/src/c/classes/Elements/Tetra.cpp

    r19395 r20608  
    997997        _assert_(gauss->Enum()==GaussTetraEnum);
    998998        this->StrainRateFS(&epsilon[0],xyz_list,(GaussTetra*)gauss,vx_input,vy_input,vz_input);
    999         this->ViscosityFS(&viscosity,3,xyz_list,(GaussTetra*)gauss,vx_input,vy_input,vz_input);
     999        this->material->ViscosityFS(&viscosity,3,xyz_list,(GaussTetra*)gauss,vx_input,vy_input,vz_input);
    10001000        GetPhi(&phi,&epsilon[0],viscosity);
    10011001

  • issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

    r20461 r20608  
    498498                /*Compute strain rate and viscosity: */
    499499                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    500                 this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
     500                this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
    501501
    502502                /*Compute Stress*/
    … …  
    555555                        /*Compute strain rate viscosity and pressure: */
    556556                        this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    557                         this->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,NULL);
     557                        this->material->ViscosityFS(&viscosity,dim,xyz_list,gauss,vx_input,vy_input,NULL);
    558558                        pressure_input->GetInputValue(&pressure,gauss);
    559559
    … …  
    611611                /*Compute strain rate viscosity and pressure: */
    612612                this->StrainRateSSA(&epsilon[0],&xyz_list[0][0],gauss,vx_input,vy_input);
    613                 this->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
     613                this->material->ViscositySSA(&viscosity,dim,&xyz_list[0][0],gauss,vx_input,vy_input);
    614614                pressure_input->GetInputValue(&pressure,gauss);
    615615
    … …  
    906906                        //              if(this->inputs->Min(MaskGroundediceLevelsetEnum)==0.){
    907907                        this->StrainRateSSA(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
    908                         this->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
     908                        this->material->ViscosityFS(&viscosity,2,xyz_list,gauss,vx_input,vy_input,NULL);
    909909                        pressure_input->GetInputValue(&pressure, gauss);
    910910                        base_input->GetInputValue(&base, gauss);

  • issm/trunk-jpl/src/c/classes/Materials/Material.h

    r18946 r20608  
    1414class Element;
    1515class Elements;
     16class Gauss;
     17class Input;
    1618/*}}}*/
    1719
    … …  
    4345                virtual void       ResetHooks()=0;
    4446
     47                virtual void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input)=0;
     48                virtual void       ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;
     49                virtual void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
     50                virtual void       ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;
     51                virtual void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf)=0;
     52                virtual void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input)=0;
     53                virtual void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon)=0;
     54
    4555};
    4656#endif

  • issm/trunk-jpl/src/c/classes/Materials/Matice.cpp

    r19282 r20608  
    337337                                                  2 eps_eff ^[(n-1)/n]
    338338
    339           where viscosity is the viscotiy, B the flow law parameter , eps_eff is the effective strain rate
    340           and n the flow law exponent.
     339          where B the flow law parameter, eps_eff is the effective strain rate and n the flow law exponent.
    341340
    342341          If eps_eff = 0 , it means this is the first time SystemMatrices is being run, and we
    … …  
    580579}
    581580/*}}}*/
     581void  Matice::ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){/*{{{*/
     582        /*The effective strain rate is defined in Paterson 3d Ed p 91 eq 9,
     583         * and Cuffey p 303 eq 8.18:
     584         *
     585         *  2 eps_eff^2 = eps_xx^2 + eps_yy^2 + eps_zz^2 + 2(eps_xy^2 + eps_xz^2 + eps_yz^2)
     586         *
     587         *  or
     588         *
     589         *  eps_eff = 1/sqrt(2) sqrt( \sum_ij eps_ij^2 )
     590         *
     591         *          = 1/sqrt(2) ||eps||_F
     592         *
     593         *  where ||.||_F is the Frobenius norm */
     594
     595        /*Intermediaries*/
     596        IssmDouble viscosity;
     597        IssmDouble epsilon3d[6]; /* epsilon=[exx,eyy,ezz,exy,exz,eyz];*/
     598        IssmDouble epsilon2d[3]; /* epsilon=[exx,eyy,exy];            */
     599        IssmDouble eps_eff;
     600        IssmDouble eps0=1.e-27;
     601
     602        if(dim==3){
     603                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
     604                element->StrainRateFS(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input,vz_input);
     605                eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[3]*epsilon3d[3] +  epsilon3d[4]*epsilon3d[4] + epsilon3d[5]*epsilon3d[5] + epsilon3d[0]*epsilon3d[1]+eps0*eps0);
     606        }
     607        else{
     608                /* eps_eff^2 = 1/2 ( exx^2 + eyy^2 + 2*exy^2 )*/
     609                element->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
     610                eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + 2.*epsilon2d[2]*epsilon2d[2]);
     611        }
     612
     613        /*Get viscosity*/
     614        this->GetViscosity(&viscosity,eps_eff);
     615
     616        /*Assign output pointer*/
     617        *pviscosity=viscosity;
     618}
     619/*}}}*/
     620void  Matice::ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
     621        this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
     622}/*}}}*/
     623void  Matice::ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     624
     625        /*Intermediaries*/
     626        IssmDouble viscosity;
     627        IssmDouble epsilon3d[5];/* epsilon=[exx,eyy,exy,exz,eyz];*/
     628        IssmDouble epsilon2d[2];/* epsilon=[exx,exy];            */
     629        IssmDouble eps_eff;
     630
     631        if(dim==3){
     632                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exz^2 + eyz^2 + exx*eyy */
     633                element->StrainRateHO(&epsilon3d[0],xyz_list,gauss,vx_input,vy_input);
     634                eps_eff = sqrt(epsilon3d[0]*epsilon3d[0] + epsilon3d[1]*epsilon3d[1] + epsilon3d[2]*epsilon3d[2] +  epsilon3d[3]*epsilon3d[3] + epsilon3d[4]*epsilon3d[4] + epsilon3d[0]*epsilon3d[1]);
     635        }
     636        else{
     637                /* eps_eff^2 = 1/2 (exx^2 + 2*exy^2 )*/
     638                element->StrainRateHO2dvertical(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
     639                eps_eff = 1./sqrt(2.)*sqrt(epsilon2d[0]*epsilon2d[0] + 2.*epsilon2d[1]*epsilon2d[1]);
     640        }
     641
     642        /*Get viscosity*/
     643        this->GetViscosity(&viscosity,eps_eff);
     644
     645        /*Assign output pointer*/
     646        *pviscosity=viscosity;
     647}/*}}}*/
     648void  Matice::ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
     649        this->GetViscosityDerivativeEpsSquare(pmu_prime,epsilon);
     650}/*}}}*/
     651void  Matice::ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surface_input){/*{{{*/
     652        /*Compute the L1L2 viscosity
     653         *
     654         *      1
     655         * mu = - A^-1 (sigma'_e)^(1-n)
     656         *      2
     657         *
     658         * sigma'_e^2 = |sigma'_//|^2 + |sigma'_perp|^2 (see Perego 2012 eq. 17,18)
     659         *
     660         * L1L2 assumptions:
     661         *
     662         * (1) |eps_b|_// = A (|sigma'_//|^2 + |sigma'_perp|^2)^((n-1)/2) |sigma'_//|
     663         * (2) |sigma'_perp|^2 = |rho g (s-z) grad(s)|^2
     664         *
     665         * Assuming that n = 3, we have a polynom of degree 3 to solve (the only unkown is X=|sigma'_//|)
     666         *
     667         * A X^3 + A |rho g (s-z) grad(s)|^2 X - |eps_b|_// = 0     */
     668
     669        IssmDouble z,s,viscosity,p,q,delta;
     670        IssmDouble tau_perp,tau_par,eps_b,A;
     671        IssmDouble epsilon[5];   /*exx eyy exy exz eyz*/
     672        IssmDouble slope[3];
     673
     674        /*Check that both inputs have been found*/
     675        if (!vx_input || !vy_input || !surface_input) _error_("Input missing");
     676
     677        /*Get tau_perp*/
     678        surface_input->GetInputValue(&s,gauss);
     679        surface_input->GetInputDerivativeValue(&slope[0],xyz_list,gauss);
     680        z=this->element->GetZcoord(xyz_list,gauss);
     681        tau_perp = element->matpar->GetMaterialParameter(MaterialsRhoIceEnum) * element->matpar->GetMaterialParameter(ConstantsGEnum) * fabs(s-z)*sqrt(slope[0]*slope[0]+slope[1]*slope[1]);
     682
     683        /* Get eps_b*/
     684        element->StrainRateHO(&epsilon[0],xyz_list,gauss,vx_input,vy_input);
     685        eps_b = sqrt(epsilon[0]*epsilon[0] + epsilon[1]*epsilon[1] + epsilon[0]*epsilon[1] + epsilon[2]*epsilon[2]);
     686        if(eps_b==0.){
     687                *pviscosity = 2.5e+17;
     688                return;
     689        }
     690
     691        /*Get A*/
     692        _assert_(this->GetN()==3.0);
     693        A=this->GetA();
     694
     695        /*Solve for tau_perp (http://fr.wikipedia.org/wiki/Méthode_de_Cardan)*/
     696        p     = tau_perp *tau_perp;
     697        q     = - eps_b/A;
     698        delta = q *q + p*p*p*4./27.;
     699        _assert_(delta>0);
     700        tau_par = pow(0.5*(-q+sqrt(delta)),1./3.) - pow(0.5*(q+sqrt(delta)),1./3.);
     701
     702        /*Viscosity*/
     703        viscosity = 1./(2.*A*(tau_par*tau_par + tau_perp*tau_perp));
     704        _assert_(!xIsNan(viscosity));
     705        _assert_(viscosity > 0.);
     706
     707        /*Assign output pointer*/
     708        *pviscosity = viscosity;
     709}/*}}}*/
     710void  Matice::ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){/*{{{*/
     711
     712        /*Intermediaries*/
     713        IssmDouble viscosity;
     714        IssmDouble epsilon2d[3];/* epsilon=[exx,eyy,exy];    */
     715        IssmDouble epsilon1d;   /* epsilon=[exx];    */
     716        IssmDouble eps_eff;
     717
     718        if(dim==2){
     719                /* eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy*/
     720                element->StrainRateSSA(&epsilon2d[0],xyz_list,gauss,vx_input,vy_input);
     721                eps_eff = sqrt(epsilon2d[0]*epsilon2d[0] + epsilon2d[1]*epsilon2d[1] + epsilon2d[2]*epsilon2d[2] + epsilon2d[0]*epsilon2d[1]);
     722        }
     723        else{
     724                /* eps_eff^2 = 1/2 exx^2*/
     725                element->StrainRateSSA1d(&epsilon1d,xyz_list,gauss,vx_input);
     726                eps_eff = sqrt(epsilon1d*epsilon1d/2.);
     727        }
     728
     729        /*Get viscosity*/
     730        this->GetViscosityBar(&viscosity,eps_eff);
     731
     732        /*Assign output pointer*/
     733        *pviscosity=viscosity;
     734}/*}}}*/
     735void  Matice::ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){/*{{{*/
     736        this->GetViscosity2dDerivativeEpsSquare(pmu_prime,epsilon);
     737}/*}}}*/
    582738void  Matice::ResetHooks(){/*{{{*/
    583739

  • issm/trunk-jpl/src/c/classes/Materials/Matice.h

    r19216 r20608  
    1818class Materials;
    1919class Parameters;
     20class Gauss;
     21class Input;
    2022/*}}}*/
    2123
    … …  
    7274                bool       IsDamage();
    7375                void       ResetHooks();
     76
     77                void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input);
     78                void       ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
     79                void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     80                void       ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
     81                void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf);
     82                void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input);
     83                void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon);
    7484                /*}}}*/
    7585};

  • issm/trunk-jpl/src/c/classes/Materials/Matpar.h

    r19984 r20608  
    115115                bool       IsDamage(){_error_("not supported");};
    116116                void       ResetHooks();
     117
     118                void       ViscosityFS(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* vz_input){_error_("not supported");};
     119                void       ViscosityFSDerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");};
     120                void       ViscosityHO(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
     121                void       ViscosityHODerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");};
     122                void       ViscosityL1L2(IssmDouble* pviscosity,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input,Input* surf){_error_("not supported");};
     123                void       ViscositySSA(IssmDouble* pviscosity,int dim,IssmDouble* xyz_list,Gauss* gauss,Input* vx_input,Input* vy_input){_error_("not supported");};
     124                void       ViscositySSADerivativeEpsSquare(IssmDouble* pmu_prime,IssmDouble* epsilon){_error_("not supported");};
    117125                /*}}}*/
    118126                /*Numerics: {{{*/
Note: See TracChangeset for help on using the changeset viewer.