Changeset 24711


Ignore:
Timestamp:
04/16/20 02:44:38 (5 years ago)
Author:
bdef
Message:

NEW:Finishing implementation of substeps fro 3D

Location:
issm/trunk-jpl
Files:
5 edited

Legend:

Unmodified
Added
Removed

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

    r24709 r24711  
    10061006        /*Update Levelset*/
    10071007        this->AddInput2(distanceenum,&ls[0],P1Enum);
     1008}
     1009/*}}}*/
     1010void       Penta::CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time){/*{{{*/
     1011
     1012        _assert_(end_time>start_time);
     1013
     1014        /*Intermediaries*/
     1015        IssmDouble averaged_values[NUMVERTICES];
     1016        IssmDouble current_values[NUMVERTICES];
     1017        IssmDouble dt;
     1018        int        found,start_offset,end_offset;
     1019        int        averaging_method=0;
     1020
     1021
     1022        /*Get transient input time steps*/
     1023        int         numtimesteps;
     1024        IssmDouble *timesteps    = NULL;
     1025        TransientInput2* transient_input  = this->inputs2->GetTransientInput(transientinput_enum);
     1026        transient_input->GetAllTimes(&timesteps,&numtimesteps);
     1027
     1028        /*go through the timesteps, and grab offset for start and end*/
     1029        found=binary_search(&start_offset,start_time,timesteps,numtimesteps);
     1030        if(!found) _error_("Input not found (is TransientInput sorted ?)");
     1031        found=binary_search(&end_offset,end_time,timesteps,numtimesteps);
     1032        if(!found) _error_("Input not found (is TransientInput sorted ?)");
     1033
     1034        Gauss* gauss=this->NewGauss();
     1035
     1036        /*stack the input for each timestep in the slice*/
     1037        int offset = start_offset;
     1038        while(offset < end_offset ){
     1039
     1040
     1041                if(offset==-1){
     1042                        /*get values for the first time: */
     1043                        _assert_(start_time<timesteps[0]);
     1044                        PentaInput2* input = transient_input->GetPentaInput(0);
     1045
     1046                        int interpolation = input->GetInterpolation();
     1047                        _assert_(interpolation==P1Enum);
     1048                        /*Intermediaries*/
     1049                        int numindices;
     1050                        int indices[NUMVERTICES];
     1051                        numindices = NUMVERTICES;
     1052                        for(int i=0;i<NUMVERTICES;i++) indices[i] = vertices[i]->lid;
     1053                        input->Serve(numindices,&indices[0]);
     1054
     1055                        for(int iv=0;iv<NUMVERTICES;iv++){
     1056                                gauss->GaussVertex(iv);
     1057                                input->GetInputValue(&current_values[iv],gauss);
     1058                        }
     1059                        dt = timesteps[0] - start_time; _assert_(dt>0.);
     1060                }
     1061                else{
     1062                        PentaInput2* input = transient_input->GetPentaInput(offset+1);
     1063                        int interpolation = input->GetInterpolation();
     1064                        _assert_(interpolation==P1Enum);
     1065                        /*Intermediaries*/
     1066                        int numindices;
     1067                        int indices[NUMVERTICES];
     1068                        numindices = NUMVERTICES;
     1069                        for(int i=0;i<NUMVERTICES;i++) indices[i] = vertices[i]->lid;
     1070                        input->Serve(numindices,&indices[0]);
     1071
     1072                        for(int iv=0;iv<NUMVERTICES;iv++){
     1073                                gauss->GaussVertex(iv);
     1074                                input->GetInputValue(&current_values[iv],gauss);
     1075                        }
     1076                        if(offset == numtimesteps-1){
     1077                                dt = end_time - timesteps[offset]; _assert_(dt>0.);
     1078                        }
     1079                        else{
     1080                                dt = timesteps[offset+1] - timesteps[offset]; _assert_(dt>0.);
     1081                        }
     1082                }
     1083
     1084                switch(averaging_method){
     1085                        case 0: /*Arithmetic mean*/
     1086                                if(offset==start_offset){
     1087                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
     1088                                }
     1089                                else{
     1090                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] += dt*current_values[iv];
     1091                                }
     1092                                break;
     1093                        case 1: /*Geometric mean*/
     1094                                if(offset==start_offset){
     1095                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv]  = dt*current_values[iv];
     1096                                }
     1097                                else{
     1098                                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] *= dt*current_values[iv];
     1099                                }
     1100                                break;
     1101                        case 2: /*Harmonic mean*/
     1102                                if(offset==start_offset){
     1103                                        for(int iv=0;iv<NUMVERTICES;iv++){
     1104                                                _assert_(current_values[iv]>1.e-50);
     1105                                                averaged_values[iv]  = dt*1./current_values[iv];
     1106                                        }
     1107                                }
     1108                                else{
     1109                                        for(int iv=0;iv<NUMVERTICES;iv++){
     1110                                                _assert_(current_values[iv]>1.e-50);
     1111                                                averaged_values[iv]  += dt*1./current_values[iv];
     1112                                        }
     1113                                }
     1114                                break;
     1115                        default:
     1116                                _error_("averaging method is not recognised");
     1117                }
     1118
     1119                offset+=1;
     1120        }
     1121
     1122        /*Integration done, now normalize*/
     1123        switch(averaging_method){
     1124                case 0: //Arithmetic mean
     1125                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] =  averaged_values[iv]/(end_time - start_time);
     1126                        break;
     1127                case 1: /*Geometric mean*/
     1128                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = pow(averaged_values[iv], 1./(end_time - start_time));
     1129                        break;
     1130                case 2: /*Harmonic mean*/
     1131                        for(int iv=0;iv<NUMVERTICES;iv++) averaged_values[iv] = 1./(averaged_values[iv]/(end_time - start_time));
     1132                        break;
     1133                default:
     1134                        _error_("averaging method is not recognised");
     1135        }
     1136
     1137        this->AddInput2(averagedinput_enum,&averaged_values[0],P1Enum);
     1138
     1139        /*Cleanup*/
     1140        delete gauss;
     1141        xDelete<IssmDouble>(timesteps);
    10081142}
    10091143/*}}}*/
    … …  
    21732307        delete gauss;
    21742308        return flux;
    2175        
     2309
    21762310}
    21772311/*}}}*/
    … …  
    29293063
    29303064        /*First, serarch the input: */
    2931         Input2* data=this->GetInput2(natureofdataenum); 
     3065        Input2* data=this->GetInput2(natureofdataenum);
    29323066
    29333067        /*figure out if we have the vertex id: */
    … …  
    40414175        /*Get material parameters :*/
    40424176        rho_ice=FindParam(MaterialsRhoIceEnum);
    4043         Input2* floatingmelt_input = this->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input); 
     4177        Input2* floatingmelt_input = this->GetInput2(BasalforcingsFloatingiceMeltingRateEnum); _assert_(floatingmelt_input);
    40444178        Input2* gllevelset_input = this->GetInput2(MaskGroundediceLevelsetEnum); _assert_(gllevelset_input);
    40454179        Input2* scalefactor_input = NULL;
    40464180        if(scaled==true){
    4047                 scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input); 
     4181                scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
    40484182        }
    40494183        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);
    … …  
    40904224        Input2* scalefactor_input  = NULL;
    40914225        if(scaled==true){
    4092                 scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input); 
     4226                scalefactor_input = this->GetInput2(MeshScaleFactorEnum); _assert_(scalefactor_input);
    40934227        }
    40944228        ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES);

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

    r24565 r24711  
    6868                void                            CreateDistanceInputFromSegmentlist(IssmDouble* distances,int distanceenum);
    6969                ElementMatrix* CreateBasalMassMatrix(void);
     70                void           CreateInputTimeAverage(int transientinput_enum,int averagedinput_enum,IssmDouble start_time,IssmDouble end_time);
    7071                void           ElementResponse(IssmDouble* presponse,int response_enum);
    7172                void           ElementSizes(IssmDouble* hx,IssmDouble* hy,IssmDouble* hz);

  • issm/trunk-jpl/src/c/classes/Inputs2/TransientInput2.cpp

    r24659 r24711  
    138138void TransientInput2::AddTriaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/
    139139
    140 
    141140        /*Check whether this is the last time step that we have*/
    142141        if(this->numtimesteps){
    … …  
    181180void TransientInput2::AddPentaTimeInput(IssmDouble time,int numindices,int* indices,IssmDouble* values_in,int interp_in){/*{{{*/
    182181
    183         _error_("not implemented yet, look at TransientInput2::AddTriaTimeInput");
     182        /*Check whether this is the last time step that we have*/
     183        if(this->numtimesteps){
     184                if(fabs(this->timesteps[this->numtimesteps-1]-time)<1.0e-5){
     185                        this->AddPentaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in);
     186                        return;
     187                }
     188        }
     189
     190        /*This is a new time step! we need to add it to the list*/
     191        if(this->numtimesteps>0 && time<this->timesteps[this->numtimesteps-1]) _error_("timestep values must increase sequentially");
     192
     193        IssmDouble *old_timesteps = NULL;
     194        Input2    **old_inputs    = NULL;
     195        if (this->numtimesteps > 0){
     196                old_timesteps=xNew<IssmDouble>(this->numtimesteps);
     197                xMemCpy(old_timesteps,this->timesteps,this->numtimesteps);
     198                xDelete<IssmDouble>(this->timesteps);
     199                old_inputs=xNew<Input2*>(this->numtimesteps);
     200                xMemCpy(old_inputs,this->inputs,this->numtimesteps);
     201                xDelete<Input2*>(this->inputs);
     202        }
     203
     204        this->numtimesteps=this->numtimesteps+1;
     205        this->timesteps=xNew<IssmDouble>(this->numtimesteps);
     206        this->inputs   = xNew<Input2*>(this->numtimesteps);
     207
     208        if (this->numtimesteps > 1){
     209                xMemCpy(this->inputs,old_inputs,this->numtimesteps-1);
     210                xMemCpy(this->timesteps,old_timesteps,this->numtimesteps-1);
     211                xDelete(old_timesteps);
     212                xDelete<Input2*>(old_inputs);
     213        }
     214
     215        /*go ahead and plug: */
     216        this->timesteps[this->numtimesteps-1] = time;
     217        this->inputs[this->numtimesteps-1]    = NULL;
     218        this->AddPentaTimeInput(this->numtimesteps-1,numindices,indices,values_in,interp_in);
    184219
    185220}

  • issm/trunk-jpl/test/NightlyRun/runme.py

    r24565 r24711  
    157157                print(("File {} saved. \n".format(os.path.join('..', 'Archives', archive_name + '.arch'))))
    158158
    159                     #ELSE: CHECK TEST
     159            #ELSE: CHECK TEST
    160160            else:
    161161                #load archive

  • issm/trunk-jpl/test/NightlyRun/test701.py

    r24261 r24711  
    1313md = bamgflowband(model(), x, b + h, b, 'hmax', 80.)
    1414
    15 print(isinstance(md, model))
    16 
    1715#Geometry  #interp1d returns a function to be called on md.mesh.x
    1816md.geometry.surface = np.interp(md.mesh.x, x, b + h)
    … …  
    2119
    2220#mask
    23 md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices, ))
     21md.mask.ice_levelset = -np.ones((md.mesh.numberofvertices))
    2422md.mask.ice_levelset[np.where(md.mesh.vertexflags(2))] = 0.
    25 md.mask.groundedice_levelset = np.zeros((md.mesh.numberofvertices, )) - 0.5
     23md.mask.groundedice_levelset = np.zeros((md.mesh.numberofvertices)) - 0.5
    2624
    2725#materials
    28 md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices, ))
     26md.initialization.temperature = (273. - 20.) * np.ones((md.mesh.numberofvertices))
    2927md.materials.rheology_B = paterson(md.initialization.temperature)
    30 md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements, ))
     28md.materials.rheology_n = 3. * np.ones((md.mesh.numberofelements))
    3129
    3230#friction
    33 md.friction.coefficient = np.zeros((md.mesh.numberofvertices, ))
     31md.friction.coefficient = np.zeros((md.mesh.numberofvertices))
    3432md.friction.coefficient[np.where(md.mesh.vertexflags(1))] = 20.
    35 md.friction.p = np.ones((md.mesh.numberofelements, ))
    36 md.friction.q = np.ones((md.mesh.numberofelements, ))
     33md.friction.p = np.ones((md.mesh.numberofelements))
     34md.friction.q = np.ones((md.mesh.numberofelements))
    3735
    3836#Boundary conditions
    … …  
    4745
    4846#Misc
    49 print(isinstance(md, model))
    50 print(type(md))
    5147md = setflowequation(md, 'FS', 'all')
    5248md.stressbalance.abstol = np.nan
    … …  
    6763for i in ['MINI', 'MINIcondensed', 'TaylorHood', 'LATaylorHood', 'CrouzeixRaviart', 'LACrouzeixRaviart']:
    6864    print(' ')
    69     print(' == == == Testing ' + i + ' Full - Stokes Finite element == == = ')
     65    print('=====Testing ' + i + ' Full-Stokes Finite element=====')
    7066    md.flowequation.fe_FS = i
    7167    md = solve(md, 'Stressbalance')
Note: See TracChangeset for help on using the changeset viewer.