Ignore:
Timestamp:
05/10/18 10:24:27 (7 years ago)
Author:
Mathieu Morlighem
Message:

merged trunk-jpl and trunk for revision 22757

Location:
issm/trunk
Files:
4 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk

  • issm/trunk/src

  • issm/trunk/src/c

    • Property svn:ignore

      •  

        old new  
        2020kriging
        2121issm_slr
         22issm_ocean
         23lnb_param.mod
         24lovenb_sub.mod
         25model.mod
         26util.mod

  • issm/trunk/src/c/analyses/EnthalpyAnalysis.cpp

    r21729 r22758  
    2020        /*Output*/
    2121        IssmDouble *spcvector  = NULL;
     22        IssmDouble *spcvectorstatic  = NULL;
    2223        IssmDouble* times=NULL;
    2324        IssmDouble* values=NULL;
     25        IssmDouble* issurface = NULL;
    2426
    2527        /*Fetch parameters: */
    … …  
    3234
    3335        /*Fetch data: */
     36        iomodel->FetchData(&issurface,&M,&N,"md.mesh.vertexonsurface"); _assert_(N>0); _assert_(M==iomodel->numberofvertices);
    3437        iomodel->FetchData(&spcvector,&M,&N,"md.thermal.spctemperature");
    35 
    36         /*Convert spcs from temperatures to enthalpy*/
    37         _assert_(N>0); _assert_(M>=iomodel->numberofvertices);
    38         for(int i=0;i<iomodel->numberofvertices;i++){
    39                 for(int j=0;j<N;j++){
    40                         spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
    41                 }
    42         }
     38        iomodel->FetchData(&spcvectorstatic,&M,&N,"md.thermal.spctemperature");
    4339
    4440        /*Specific case for PDD, we want the constaints to be updated by the PDD scheme itself*/
    … …  
    5147        }
    5248
     49        /*Convert spcs from temperatures to enthalpy*/
     50        _assert_(N>0); _assert_(M>=iomodel->numberofvertices);
     51        for(int i=0;i<iomodel->numberofvertices;i++){
     52                for(int j=0;j<N;j++){
     53                        if (isdynamic){
     54                                if (issurface[i]==1){
     55                                        spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
     56                                        spcvectorstatic[i*N+j] = NAN;
     57                                }
     58                                else{
     59                                        spcvector[i*N+j] = NAN;
     60                                        spcvectorstatic[i*N+j] = heatcapacity*(spcvectorstatic[i*N+j]-referencetemperature);
     61                                }
     62                        }
     63                        else{
     64                                spcvector[i*N+j] = heatcapacity*(spcvector[i*N+j]-referencetemperature);
     65                        }
     66                }
     67        }
     68
    5369        if(isdynamic){
    54                 IoModelToDynamicConstraintsx(constraints,iomodel,spcvector,M,N,EnthalpyAnalysisEnum,finiteelement);
     70                IoModelToDynamicConstraintsx(constraints,iomodel,spcvector,iomodel->numberofvertices,1,EnthalpyAnalysisEnum,finiteelement);
     71                IoModelToConstraintsx(constraints,iomodel,spcvectorstatic,M,N,EnthalpyAnalysisEnum,finiteelement);
    5572        }
    5673        else{
    … …  
    6077        /*Free ressources:*/
    6178        iomodel->DeleteData(spcvector,"md.thermal.spctemperature");
     79        iomodel->DeleteData(spcvectorstatic,"md.thermal.spctemperature");
     80        iomodel->DeleteData(issurface,"md.mesh.vertexonsurface");
    6281        xDelete<IssmDouble>(times);
    6382        xDelete<IssmDouble>(values);
    … …  
    174193        switch(frictionlaw){
    175194                case 1:
     195                        iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
    176196                        iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
    177197                        iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
    178198                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
     199                        if (FrictionCoupling==1){
     200                                iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
     201                        }
    179202                        break;
    180203                case 2:
    … …  
    187210                        iomodel->FetchDataToInput(elements,"md.friction.As",FrictionAsEnum);
    188211                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
    189                         if (FrictionCoupling==0){
     212                        if (FrictionCoupling==1){
    190213                                iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
    191214                        }
    … …  
    210233                        iomodel->FetchDataToInput(elements,"md.initialization.temperature",TemperatureEnum);
    211234                        break;
     235                case 7:
     236                        iomodel->FindConstant(&FrictionCoupling,"md.friction.coupling");
     237                        iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
     238                        iomodel->FetchDataToInput(elements,"md.friction.coefficientcoulomb",FrictionCoefficientcoulombEnum);
     239                        iomodel->FetchDataToInput(elements,"md.friction.p",FrictionPEnum);
     240                        iomodel->FetchDataToInput(elements,"md.friction.q",FrictionQEnum);
     241                        if (FrictionCoupling==1){
     242                                iomodel->FetchDataToInput(elements,"md.friction.effective_pressure",FrictionEffectivePressureEnum);
     243                        }
     244                        break;
    212245                case 9:
    213246                        iomodel->FetchDataToInput(elements,"md.friction.coefficient",FrictionCoefficientEnum);
    … …  
    217250                        break;
    218251                default:
    219                         _error_("not supported");
     252                        _error_("friction law not supported");
    220253        }
    221254
    … …  
    244277        int frictionlaw;
    245278        iomodel->FindConstant(&frictionlaw,"md.friction.law");
    246         if(frictionlaw==4 || frictionlaw==6) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
    247         if(frictionlaw==3) parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
    248         if(frictionlaw==9) parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
     279        if(frictionlaw==4 || frictionlaw==6){
     280                parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
     281        }
     282        if(frictionlaw==3 || frictionlaw==1 || frictionlaw==7){
     283                parameters->AddObject(iomodel->CopyConstantObject("md.friction.coupling",FrictionCouplingEnum));
     284        }
     285        if(frictionlaw==9){
     286                parameters->AddObject(iomodel->CopyConstantObject("md.friction.gamma",FrictionGammaEnum));
     287                parameters->AddObject(new IntParam(FrictionCouplingEnum,0));
     288        }
    249289}/*}}}*/
    250290
    … …  
    302342                ComputeBasalMeltingrate(element);
    303343        }
     344
     345        /*extrude inputs*/
     346        femmodel->parameters->SetParam(BasalforcingsGroundediceMeltingRateEnum,InputToExtrudeEnum);
     347        extrudefrombase_core(femmodel);
    304348}/*}}}*/
    305349void           EnthalpyAnalysis::ComputeBasalMeltingrate(Element* element){/*{{{*/
    … …  
    525569        IssmDouble*    basis    = xNew<IssmDouble>(numnodes);
    526570        IssmDouble*    dbasis   = xNew<IssmDouble>(3*numnodes);
    527         IssmDouble*    B        = xNew<IssmDouble>(3*numnodes);
    528571        IssmDouble*    Bprime   = xNew<IssmDouble>(3*numnodes);
    529         IssmDouble     D[3][3]  = {0.};
    530572        IssmDouble     K[3][3];
    531573
    … …  
    556598
    557599                element->JacobianDeterminant(&Jdet,xyz_list,gauss);
     600                element->NodalFunctions(basis,gauss);
     601                element->NodalFunctionsDerivatives(dbasis,xyz_list,gauss);
     602
    558603                D_scalar=gauss->weight*Jdet;
    559604                if(dt!=0.) D_scalar=D_scalar*dt;
    560605
    561606                /*Conduction: */
    562                 GetBConduct(B,element,xyz_list,gauss);
    563                 D[0][0]=D_scalar*kappa/rho_ice;
    564                 D[1][1]=D_scalar*kappa/rho_ice;
    565                 D[2][2]=D_scalar*kappa/rho_ice;
    566                 TripleMultiply(B,3,numnodes,1,
    567                                         &D[0][0],3,3,0,
    568                                         B,3,numnodes,0,
    569                                         &Ke->values[0],1);
     607                for(int i=0;i<numnodes;i++){
     608                        for(int j=0;j<numnodes;j++){
     609                                Ke->values[i*numnodes+j] += D_scalar*kappa/rho_ice*(
     610                                                        dbasis[0*numnodes+j]*dbasis[0*numnodes+i] + dbasis[1*numnodes+j]*dbasis[1*numnodes+i] + dbasis[2*numnodes+j]*dbasis[2*numnodes+i]
     611                                                        );
     612                        }
     613                }
    570614
    571615                /*Advection: */
    572                 GetBAdvec(B,element,xyz_list,gauss);
    573                 GetBAdvecprime(Bprime,element,xyz_list,gauss);
    574616                vx_input->GetInputValue(&u,gauss); vxm_input->GetInputValue(&um,gauss); vx=u-um;
    575617                vy_input->GetInputValue(&v,gauss); vym_input->GetInputValue(&vm,gauss); vy=v-vm;
    576618                vz_input->GetInputValue(&w,gauss); vzm_input->GetInputValue(&wm,gauss); vz=w-wm;
    577                 D[0][0]=D_scalar*vx;
    578                 D[1][1]=D_scalar*vy;
    579                 D[2][2]=D_scalar*vz;
    580                 TripleMultiply(B,3,numnodes,1,
    581                                         &D[0][0],3,3,0,
    582                                         Bprime,3,numnodes,0,
    583                                         &Ke->values[0],1);
     619                for(int i=0;i<numnodes;i++){
     620                        for(int j=0;j<numnodes;j++){
     621                                Ke->values[i*numnodes+j] += D_scalar*(
     622                                                        vx*dbasis[0*numnodes+j]*basis[i] + vy*dbasis[1*numnodes+j]*basis[i] +vz*dbasis[2*numnodes+j]*basis[i]
     623                                                        );
     624                        }
     625                }
    584626
    585627                /*Transient: */
    586628                if(dt!=0.){
    587629                        D_scalar=gauss->weight*Jdet;
    588                         element->NodalFunctions(basis,gauss);
    589                         TripleMultiply(basis,numnodes,1,0,
    590                                                 &D_scalar,1,1,0,
    591                                                 basis,1,numnodes,0,
    592                                                 &Ke->values[0],1);
     630                        for(int i=0;i<numnodes;i++){
     631                                for(int j=0;j<numnodes;j++){
     632                                        Ke->values[i*numnodes+j] += D_scalar*basis[j]*basis[i];
     633                                }
     634                        }
    593635                        D_scalar=D_scalar*dt;
    594636                }
    … …  
    634676        xDelete<IssmDouble>(basis);
    635677        xDelete<IssmDouble>(dbasis);
    636         xDelete<IssmDouble>(B);
    637678        xDelete<IssmDouble>(Bprime);
    638679        delete gauss;
    … …  
    10631104                xDelete<IssmDouble>(watercolumn);
    10641105                xDelete<IssmDouble>(drainage_int);
    1065         }
    1066         xDelete<int>(basalnodeindices);
     1106                xDelete<int>(basalnodeindices);
     1107        }
    10671108}/*}}}*/
    10681109void                            EnthalpyAnalysis::DrainageUpdateEnthalpy(FemModel* femmodel){/*{{{*/
    … …  
    15101551                /*Update Rheology only if converged (we must make sure that the temperature is below melting point
    15111552                 * otherwise the rheology could be negative*/
    1512                 element->FindParam(&rheology_law,MaterialsRheologyLawEnum);
     1553                rheology_law=element->GetIntegerMaterialParameter(MaterialsRheologyLawEnum);
    15131554                element->GetInputListOnNodes(&surface[0],SurfaceEnum);
    15141555                switch(rheology_law){
Note: See TracChangeset for help on using the changeset viewer.