Changeset 16773


Ignore:
Timestamp:
11/14/13 20:04:08 (11 years ago)
Author:
seroussi
Message:

BUG: fixed InputUpdateFromSolution for SSA HO tiling

Location:
issm/trunk-jpl/src/c/analyses
Files:
2 edited

Legend:

Unmodified
Added
Removed

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

    r16769 r16773  
    943943                        InputUpdateFromSolutionHO(solution,element);
    944944                        return;
    945                 case SIAApproximationEnum:
    946                         InputUpdateFromSolutionHO(solution,element);
    947                         return;
    948945                case L1L2ApproximationEnum:
    949946                        InputUpdateFromSolutionSSA(solution,element);
    … …  
    961958        }
    962959}/*}}}*/
    963 void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
     960void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
     961
     962        int          i,dim,meshtype;
     963        int*         vdoflist=NULL;
     964        int*         pdoflist=NULL;
     965        IssmDouble   FSreconditioning;
     966
     967        element->FindParam(&meshtype,MeshTypeEnum);
     968        element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
     969        switch(meshtype){
     970                case Mesh2DverticalEnum: dim = 2; break;
     971                case Mesh3DEnum:         dim = 3; break;
     972                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     973        }
     974
     975        /*Fetch number of nodes and dof for this finite element*/
     976        int vnumnodes = element->GetNumberOfNodesVelocity();
     977        int pnumnodes = element->GetNumberOfNodesPressure();
     978        int vnumdof   = vnumnodes*dim;
     979        int pnumdof   = pnumnodes*1;
     980
     981        /*Initialize values*/
     982        IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
     983        IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
     984        IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
     985        IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
     986        IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
     987        IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
     988
     989        /*Prepare coordinate system list*/
     990        int* cs_list = xNew<int>(vnumnodes+pnumnodes);
     991        if(dim==2){
     992                for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
     993        }
     994        else{
     995                for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
     996        }
     997        for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
     998
     999        /*Get dof list: */
     1000        element->GetDofListVelocity(&vdoflist,GsetEnum);
     1001        element->GetDofListPressure(&pdoflist,GsetEnum);
     1002
     1003        /*Use the dof list to index into the solution vector: */
     1004        for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
     1005        for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
     1006
     1007        /*Transform solution in Cartesian Space*/
     1008        element->TransformSolutionCoord(&values[0],&cs_list[0]);
     1009
     1010        /*Ok, we have vx and vy in values, fill in all arrays: */
     1011        for(i=0;i<vnumnodes;i++){
     1012                vx[i] = values[i*dim+0];
     1013                vy[i] = values[i*dim+1];
     1014                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
     1015                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     1016
     1017                if(dim==3){
     1018                        vz[i] = values[i*dim+2];
     1019                        if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
     1020                }
     1021        }
     1022        for(i=0;i<pnumnodes;i++){
     1023                pressure[i] = values[vnumdof+i];
     1024                if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
     1025        }
     1026
     1027        /*Recondition pressure and compute vel: */
     1028        for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
     1029        if(dim==3) for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     1030        else       for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
     1031
     1032        /*Now, we have to move the previous inputs  to old
     1033         * status, otherwise, we'll wipe them off: */
     1034        element->InputChangeName(VxEnum,VxPicardEnum);
     1035        element->InputChangeName(VyEnum,VyPicardEnum);
     1036        element->InputChangeName(PressureEnum,PressurePicardEnum);
     1037        if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);
     1038
     1039        /*Add vx and vy as inputs to the tria element: */
     1040        element->AddInput(VxEnum,vx,P1Enum);
     1041        element->AddInput(VyEnum,vy,P1Enum);
     1042        element->AddInput(VelEnum,vel,P1Enum);
     1043        element->AddInput(PressureEnum,pressure,P1Enum);
     1044        if(dim==3) element->AddInput(VzEnum,vz,P1Enum);
     1045
     1046        /*Free ressources:*/
     1047        xDelete<IssmDouble>(pressure);
     1048        xDelete<IssmDouble>(vel);
     1049        xDelete<IssmDouble>(vz);
     1050        xDelete<IssmDouble>(vy);
     1051        xDelete<IssmDouble>(vx);
     1052        xDelete<IssmDouble>(values);
     1053        xDelete<int>(vdoflist);
     1054        xDelete<int>(pdoflist);
     1055        xDelete<int>(cs_list);
     1056}/*}}}*/
     1057void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
     1058
     1059        int         i,meshtype;
     1060        IssmDouble  rho_ice,g;
     1061        int*        doflist=NULL;
     1062        IssmDouble* xyz_list=NULL;
     1063
     1064        /*Fetch number of nodes and dof for this finite element*/
     1065        int numnodes = element->GetNumberOfNodes();
     1066        int numdof   = numnodes*2;
     1067
     1068        /*Fetch dof list and allocate solution vectors*/
     1069        element->GetDofList(&doflist,HOApproximationEnum,GsetEnum);
     1070        IssmDouble* values    = xNew<IssmDouble>(numdof);
     1071        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
     1072        IssmDouble* vy        = xNew<IssmDouble>(numnodes);
     1073        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
     1074        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
     1075        IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
     1076        IssmDouble* surface   = xNew<IssmDouble>(numnodes);
     1077
     1078        /*Use the dof list to index into the solution vector: */
     1079        for(i=0;i<numdof;i++) values[i]=solution[doflist[i]];
     1080
     1081        /*Transform solution in Cartesian Space*/
     1082        element->TransformSolutionCoord(&values[0],XYEnum);
     1083        element->FindParam(&meshtype,MeshTypeEnum);
     1084
     1085        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
     1086        for(i=0;i<numnodes;i++){
     1087                vx[i]=values[i*2+0];
     1088                vy[i]=values[i*2+1];
     1089
     1090                /*Check solution*/
     1091                if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
     1092                if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
     1093        }
     1094
     1095        /*Get Vz and compute vel*/
     1096        element->GetInputListOnNodes(&vz[0],VzEnum,0.);
     1097        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
     1098
     1099        /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
     1100         *so the pressure is just the pressure at the bedrock: */
     1101        rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
     1102        g       = element->GetMaterialParameter(ConstantsGEnum);
     1103        element->GetVerticesCoordinates(&xyz_list);
     1104        element->GetInputListOnNodes(&surface[0],SurfaceEnum);
     1105        for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
     1106
     1107        /*Now, we have to move the previous Vx and Vy inputs  to old
     1108         * status, otherwise, we'll wipe them off: */
     1109        element->InputChangeName(VxEnum,VxPicardEnum);
     1110        element->InputChangeName(VyEnum,VyPicardEnum);
     1111        element->InputChangeName(PressureEnum,PressurePicardEnum);
     1112
     1113        /*Add vx and vy as inputs to the tria element: */
     1114        element->AddInput(VxEnum,vx,P1Enum);
     1115        element->AddInput(VyEnum,vy,P1Enum);
     1116        element->AddBasalInput(VelEnum,vel,P1Enum);
     1117        element->AddBasalInput(PressureEnum,pressure,P1Enum);
     1118
     1119        /*Free ressources:*/
     1120        xDelete<IssmDouble>(surface);
     1121        xDelete<IssmDouble>(pressure);
     1122        xDelete<IssmDouble>(vel);
     1123        xDelete<IssmDouble>(vz);
     1124        xDelete<IssmDouble>(vy);
     1125        xDelete<IssmDouble>(vx);
     1126        xDelete<IssmDouble>(values);
     1127        xDelete<IssmDouble>(xyz_list);
     1128        xDelete<int>(doflist);
     1129}/*}}}*/
     1130void StressbalanceAnalysis::InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element){/*{{{*/
    9641131
    9651132        int         i,meshtype;
    … …  
    10601227        if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
    10611228}/*}}}*/
    1062 void StressbalanceAnalysis::InputUpdateFromSolutionHO(IssmDouble* solution,Element* element){/*{{{*/
     1229void StressbalanceAnalysis::InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element){/*{{{*/
    10631230
    10641231        int         i,meshtype;
    … …  
    10661233        int*        doflist=NULL;
    10671234        IssmDouble* xyz_list=NULL;
     1235        Element*    basalelement=NULL;
     1236
     1237        /*Deal with pressure first*/
     1238        int numvertices = element->GetNumberOfVertices();
     1239        IssmDouble* pressure  = xNew<IssmDouble>(numvertices);
     1240        IssmDouble* thickness = xNew<IssmDouble>(numvertices);
     1241        IssmDouble* surface   = xNew<IssmDouble>(numvertices);
     1242
     1243        element->FindParam(&meshtype,MeshTypeEnum);
     1244        rho_ice =element->GetMaterialParameter(MaterialsRhoIceEnum);
     1245        g       =element->GetMaterialParameter(ConstantsGEnum);
     1246        switch(meshtype){
     1247                case Mesh2DhorizontalEnum:
     1248                        element->GetInputListOnVertices(thickness,ThicknessEnum);
     1249                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*thickness[i];
     1250                        break;
     1251                case Mesh3DEnum:
     1252                        element->GetVerticesCoordinates(&xyz_list);
     1253                        element->GetInputListOnVertices(surface,SurfaceEnum);
     1254                        for(i=0;i<numvertices;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
     1255                        break;
     1256                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1257        }
     1258        element->AddInput(PressureEnum,pressure,P1Enum);
     1259        xDelete<IssmDouble>(pressure);
     1260        xDelete<IssmDouble>(thickness);
     1261        xDelete<IssmDouble>(surface);
     1262
     1263        /*Get basal element*/
     1264        switch(meshtype){
     1265                case Mesh2DhorizontalEnum:
     1266                        basalelement = element;
     1267                        break;
     1268                case Mesh3DEnum:
     1269                        if(!element->IsOnBed()){xDelete<IssmDouble>(xyz_list); return;}
     1270                        basalelement=element->SpawnBasalElement();
     1271                        break;
     1272                default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
     1273        }
    10681274
    10691275        /*Fetch number of nodes and dof for this finite element*/
    1070         int numnodes = element->GetNumberOfNodes();
     1276        int numnodes = basalelement->GetNumberOfNodes();
    10711277        int numdof   = numnodes*2;
    10721278
    10731279        /*Fetch dof list and allocate solution vectors*/
    1074         element->GetDofList(&doflist,NoneApproximationEnum,GsetEnum);
     1280        basalelement->GetDofList(&doflist,SSAApproximationEnum,GsetEnum);
    10751281        IssmDouble* values    = xNew<IssmDouble>(numdof);
    10761282        IssmDouble* vx        = xNew<IssmDouble>(numnodes);
    … …  
    10781284        IssmDouble* vz        = xNew<IssmDouble>(numnodes);
    10791285        IssmDouble* vel       = xNew<IssmDouble>(numnodes);
    1080         IssmDouble* pressure  = xNew<IssmDouble>(numnodes);
    1081         IssmDouble* surface   = xNew<IssmDouble>(numnodes);
    10821286
    10831287        /*Use the dof list to index into the solution vector: */
    … …  
    10851289
    10861290        /*Transform solution in Cartesian Space*/
    1087         element->TransformSolutionCoord(&values[0],XYEnum);
    1088         element->FindParam(&meshtype,MeshTypeEnum);
     1291        basalelement->TransformSolutionCoord(&values[0],XYEnum);
     1292        basalelement->FindParam(&meshtype,MeshTypeEnum);
    10891293
    10901294        /*Ok, we have vx and vy in values, fill in vx and vy arrays: */
    … …  
    10991303
    11001304        /*Get Vz and compute vel*/
    1101         element->GetInputListOnNodes(&vz[0],VzEnum,0.);
     1305        basalelement->GetInputListOnNodes(&vz[0],VzEnum,0.);
    11021306        for(i=0;i<numnodes;i++) vel[i]=sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    1103 
    1104         /*For pressure: we have not computed pressure in this analysis, for this element. We are in 2D,
    1105          *so the pressure is just the pressure at the bedrock: */
    1106         rho_ice = element->GetMaterialParameter(MaterialsRhoIceEnum);
    1107         g       = element->GetMaterialParameter(ConstantsGEnum);
    1108         element->GetVerticesCoordinates(&xyz_list);
    1109         element->GetInputListOnNodes(&surface[0],SurfaceEnum);
    1110         for(i=0;i<numnodes;i++) pressure[i]=rho_ice*g*(surface[i]-xyz_list[i*3+2]);
    11111307
    11121308        /*Now, we have to move the previous Vx and Vy inputs  to old
    … …  
    11141310        element->InputChangeName(VxEnum,VxPicardEnum);
    11151311        element->InputChangeName(VyEnum,VyPicardEnum);
    1116         element->InputChangeName(PressureEnum,PressurePicardEnum);
    11171312
    11181313        /*Add vx and vy as inputs to the tria element: */
    1119         element->AddInput(VxEnum,vx,P1Enum);
    1120         element->AddInput(VyEnum,vy,P1Enum);
     1314        element->AddBasalInput(VxEnum,vx,P1Enum);
     1315        element->AddBasalInput(VyEnum,vy,P1Enum);
    11211316        element->AddBasalInput(VelEnum,vel,P1Enum);
    1122         element->AddBasalInput(PressureEnum,pressure,P1Enum);
    11231317
    11241318        /*Free ressources:*/
    1125         xDelete<IssmDouble>(surface);
    1126         xDelete<IssmDouble>(pressure);
    11271319        xDelete<IssmDouble>(vel);
    11281320        xDelete<IssmDouble>(vz);
    … …  
    11321324        xDelete<IssmDouble>(xyz_list);
    11331325        xDelete<int>(doflist);
    1134 }/*}}}*/
    1135 void StressbalanceAnalysis::InputUpdateFromSolutionFS(IssmDouble* solution,Element* element){/*{{{*/
    1136 
    1137         int          i,dim,meshtype;
    1138         int*         vdoflist=NULL;
    1139         int*         pdoflist=NULL;
    1140         IssmDouble   FSreconditioning;
    1141 
    1142         element->FindParam(&meshtype,MeshTypeEnum);
    1143         element->FindParam(&FSreconditioning,StressbalanceFSreconditioningEnum);
    1144         switch(meshtype){
    1145                 case Mesh2DverticalEnum: dim = 2; break;
    1146                 case Mesh3DEnum:         dim = 3; break;
    1147                 default: _error_("mesh "<<EnumToStringx(meshtype)<<" not supported yet");
    1148         }
    1149 
    1150         /*Fetch number of nodes and dof for this finite element*/
    1151         int vnumnodes = element->GetNumberOfNodesVelocity();
    1152         int pnumnodes = element->GetNumberOfNodesPressure();
    1153         int vnumdof   = vnumnodes*dim;
    1154         int pnumdof   = pnumnodes*1;
    1155 
    1156         /*Initialize values*/
    1157         IssmDouble* values   = xNew<IssmDouble>(vnumdof+pnumdof);
    1158         IssmDouble* vx       = xNew<IssmDouble>(vnumnodes);
    1159         IssmDouble* vy       = xNew<IssmDouble>(vnumnodes);
    1160         IssmDouble* vz       = xNew<IssmDouble>(vnumnodes);
    1161         IssmDouble* vel      = xNew<IssmDouble>(vnumnodes);
    1162         IssmDouble* pressure = xNew<IssmDouble>(pnumnodes);
    1163 
    1164         /*Prepare coordinate system list*/
    1165         int* cs_list = xNew<int>(vnumnodes+pnumnodes);
    1166         if(dim==2){
    1167                 for(i=0;i<vnumnodes;i++) cs_list[i] = XYEnum;
    1168         }
    1169         else{
    1170                 for(i=0;i<vnumnodes;i++) cs_list[i] = XYZEnum;
    1171         }
    1172         for(i=0;i<pnumnodes;i++) cs_list[vnumnodes+i] = PressureEnum;
    1173 
    1174         /*Get dof list: */
    1175         element->GetDofListVelocity(&vdoflist,GsetEnum);
    1176         element->GetDofListPressure(&pdoflist,GsetEnum);
    1177 
    1178         /*Use the dof list to index into the solution vector: */
    1179         for(i=0;i<vnumdof;i++) values[i]        =solution[vdoflist[i]];
    1180         for(i=0;i<pnumdof;i++) values[vnumdof+i]=solution[pdoflist[i]];
    1181 
    1182         /*Transform solution in Cartesian Space*/
    1183         element->TransformSolutionCoord(&values[0],&cs_list[0]);
    1184 
    1185         /*Ok, we have vx and vy in values, fill in all arrays: */
    1186         for(i=0;i<vnumnodes;i++){
    1187                 vx[i] = values[i*dim+0];
    1188                 vy[i] = values[i*dim+1];
    1189                 if(xIsNan<IssmDouble>(vx[i])) _error_("NaN found in solution vector");
    1190                 if(xIsNan<IssmDouble>(vy[i])) _error_("NaN found in solution vector");
    1191 
    1192                 if(dim==3){
    1193                         vz[i] = values[i*dim+2];
    1194                         if(xIsNan<IssmDouble>(vz[i])) _error_("NaN found in solution vector");
    1195                 }
    1196         }
    1197         for(i=0;i<pnumnodes;i++){
    1198                 pressure[i] = values[vnumdof+i];
    1199                 if(xIsNan<IssmDouble>(pressure[i])) _error_("NaN found in solution vector");
    1200         }
    1201 
    1202         /*Recondition pressure and compute vel: */
    1203         for(i=0;i<pnumnodes;i++) pressure[i] = pressure[i]*FSreconditioning;
    1204         if(dim==3) for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i] + vz[i]*vz[i]);
    1205         else       for(i=0;i<vnumnodes;i++) vel[i] = sqrt(vx[i]*vx[i] + vy[i]*vy[i]);
    1206 
    1207         /*Now, we have to move the previous inputs  to old
    1208          * status, otherwise, we'll wipe them off: */
    1209         element->InputChangeName(VxEnum,VxPicardEnum);
    1210         element->InputChangeName(VyEnum,VyPicardEnum);
    1211         element->InputChangeName(PressureEnum,PressurePicardEnum);
    1212         if(dim==3) element->InputChangeName(VzEnum,VzPicardEnum);
    1213 
    1214         /*Add vx and vy as inputs to the tria element: */
    1215         element->AddInput(VxEnum,vx,P1Enum);
    1216         element->AddInput(VyEnum,vy,P1Enum);
    1217         element->AddInput(VelEnum,vel,P1Enum);
    1218         element->AddInput(PressureEnum,pressure,P1Enum);
    1219         if(dim==3) element->AddInput(VzEnum,vz,P1Enum);
    1220 
    1221         /*Free ressources:*/
    1222         xDelete<IssmDouble>(pressure);
    1223         xDelete<IssmDouble>(vel);
    1224         xDelete<IssmDouble>(vz);
    1225         xDelete<IssmDouble>(vy);
    1226         xDelete<IssmDouble>(vx);
    1227         xDelete<IssmDouble>(values);
    1228         xDelete<int>(vdoflist);
    1229         xDelete<int>(pdoflist);
    1230         xDelete<int>(cs_list);
     1326        if(meshtype!=Mesh2DhorizontalEnum) delete basalelement;
    12311327}/*}}}*/
    12321328void StressbalanceAnalysis::InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element){/*{{{*/
    … …  
    12521348        /*Fetch dof list and allocate solution vectors*/
    12531349        basalelement->GetDofList(&SSAdoflist,SSAApproximationEnum,GsetEnum);
    1254         element->GetDofList(     &HOdoflist, HOApproximationEnum, GsetEnum);
     1350        element     ->GetDofList(&HOdoflist, HOApproximationEnum, GsetEnum);
    12551351        IssmDouble* HOvalues  = xNew<IssmDouble>(numdof);
    12561352        IssmDouble* SSAvalues = xNew<IssmDouble>(numdof);

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

    r16762 r16773  
    2424                void InputUpdateFromSolutionFS(IssmDouble* solution,Element* element);
    2525                void InputUpdateFromSolutionHO(IssmDouble* solution,Element* element);
     26                void InputUpdateFromSolutionL1L2(IssmDouble* solution,Element* element);
    2627                void InputUpdateFromSolutionSSA(IssmDouble* solution,Element* element);
    2728                void InputUpdateFromSolutionSSAHO(IssmDouble* solution,Element* element);
Note: See TracChangeset for help on using the changeset viewer.