Index: /issm/trunk-jpl/m4/issm_options.m4 =================================================================== --- /issm/trunk-jpl/m4/issm_options.m4 (revision 25626) +++ /issm/trunk-jpl/m4/issm_options.m4 (revision 25627) @@ -2475,9 +2475,19 @@ dnl optimization{{{ dnl -- bypass standard optimization -g -O2 -fPIC -std=c++11? + dnl NOTE: Added -Wno-deprecated-register flag to clear, + dnl + dnl externalpackages/dakota/install/include/MPIPackBuffer.hpp:301:8: warning: + dnl 'register' storage class specifier is deprecated and incompatible with + dnl C++17 [-Wdeprecated-register] + dnl for (register typename ContainerT::size_type i=0; iAddObject(iomodel->CopyConstantObject("md.solidearth.settings.computesealevelchange",SolidearthSettingsComputesealevelchangeEnum)); parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.planetradius",SolidearthPlanetRadiusEnum)); + parameters->AddObject(new DoubleParam(CumBslrEnum,0.0)); + parameters->AddObject(new DoubleParam(CumBslrIceEnum,0.0)); + parameters->AddObject(new DoubleParam(CumBslrHydroEnum,0.0)); + parameters->AddObject(new DoubleParam(CumGmtslrEnum,0.0)); /*compute planet area and plug into parameters:*/ @@ -226,145 +231,149 @@ } /*}}}*/ /*Deal with elasticity {{{*/ - /*love numbers: */ - iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); - iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); - iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l"); - iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th"); - iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); - iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); - parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); - - parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); - parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1)); - parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1)); - parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); - - /*compute elastic green function for a range of angles*/ - iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); - M=reCast(180./degacc+1.); - - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_rigid=xNew(M,"t"); - G_elastic=xNew(M,"t"); - U_elastic=xNew(M,"t"); - H_elastic=xNew(M,"t"); - #else - G_rigid=xNew(M); - G_elastic=xNew(M); - U_elastic=xNew(M); - H_elastic=xNew(M); - #endif - - /*compute combined legendre + love number (elastic green function:*/ - m=DetermineLocalSize(M,IssmComm::GetComm()); - GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); - // AD performance is sensitive to calls to ensurecontiguous. - // // Providing "t" will cause ensurecontiguous to be called. - #ifdef _HAVE_AD_ - G_elastic_local=xNew(m,"t"); - G_rigid_local=xNew(m,"t"); - U_elastic_local=xNew(m,"t"); - H_elastic_local=xNew(m,"t"); - #else - G_elastic_local=xNew(m); - G_rigid_local=xNew(m); - U_elastic_local=xNew(m); - H_elastic_local=xNew(m); - #endif - - for(int i=lower_row;i(i)*degacc * PI / 180.0; - - G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); - G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; - U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; - H_elastic_local[i-lower_row]= 0; - IssmDouble Pn = 0.; - IssmDouble Pn1 = 0.; - IssmDouble Pn2 = 0.; - IssmDouble Pn_p = 0.; - IssmDouble Pn_p1 = 0.; - IssmDouble Pn_p2 = 0.; - - for (int n=0;nFetchData(&elastic,"md.solidearth.settings.elastic"); +// if(elastic){ + /*love numbers: */ + iomodel->FetchData(&love_h,&nl,NULL,"md.solidearth.lovenumbers.h"); + iomodel->FetchData(&love_k,&nl,NULL,"md.solidearth.lovenumbers.k"); + iomodel->FetchData(&love_l,&nl,NULL,"md.solidearth.lovenumbers.l"); + iomodel->FetchData(&love_th,&nl,NULL,"md.solidearth.lovenumbers.th"); + iomodel->FetchData(&love_tk,&nl,NULL,"md.solidearth.lovenumbers.tk"); + iomodel->FetchData(&love_tl,&nl,NULL,"md.solidearth.lovenumbers.tl"); + parameters->AddObject(iomodel->CopyConstantObject("md.solidearth.lovenumbers.tk2secular",TidalLoveK2SecularEnum)); + + parameters->AddObject(new DoubleMatParam(LoadLoveHEnum,love_h,nl,1)); + parameters->AddObject(new DoubleMatParam(LoadLoveKEnum,love_k,nl,1)); + parameters->AddObject(new DoubleMatParam(LoadLoveLEnum,love_l,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveHEnum,love_th,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveKEnum,love_tk,nl,1)); + parameters->AddObject(new DoubleMatParam(TidalLoveLEnum,love_tl,nl,1)); + + /*compute elastic green function for a range of angles*/ + iomodel->FetchData(°acc,"md.solidearth.settings.degacc"); + M=reCast(180./degacc+1.); + + // AD performance is sensitive to calls to ensurecontiguous. + // // Providing "t" will cause ensurecontiguous to be called. + #ifdef _HAVE_AD_ + G_rigid=xNew(M,"t"); + G_elastic=xNew(M,"t"); + U_elastic=xNew(M,"t"); + H_elastic=xNew(M,"t"); + #else + G_rigid=xNew(M); + G_elastic=xNew(M); + U_elastic=xNew(M); + H_elastic=xNew(M); + #endif + + /*compute combined legendre + love number (elastic green function:*/ + m=DetermineLocalSize(M,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,m,IssmComm::GetComm()); + // AD performance is sensitive to calls to ensurecontiguous. + // // Providing "t" will cause ensurecontiguous to be called. + #ifdef _HAVE_AD_ + G_elastic_local=xNew(m,"t"); + G_rigid_local=xNew(m,"t"); + U_elastic_local=xNew(m,"t"); + H_elastic_local=xNew(m,"t"); + #else + G_elastic_local=xNew(m); + G_rigid_local=xNew(m); + U_elastic_local=xNew(m); + H_elastic_local=xNew(m); + #endif + + for(int i=lower_row;i(i)*degacc * PI / 180.0; + + G_rigid_local[i-lower_row]= .5/sin(alpha/2.0); + G_elastic_local[i-lower_row]= (love_k[nl-1]-love_h[nl-1])*G_rigid_local[i-lower_row]; + U_elastic_local[i-lower_row]= (love_h[nl-1])*G_rigid_local[i-lower_row]; + H_elastic_local[i-lower_row]= 0; + IssmDouble Pn = 0.; + IssmDouble Pn1 = 0.; + IssmDouble Pn2 = 0.; + IssmDouble Pn_p = 0.; + IssmDouble Pn_p1 = 0.; + IssmDouble Pn_p2 = 0.; + + for (int n=0;n(IssmComm::GetSize()); - int* displs=xNew(IssmComm::GetSize()); - - //recvcounts: - ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*displs: */ - ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); - - /*All gather:*/ - ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); - /*free ressources: */ - xDelete(recvcounts); - xDelete(displs); - - /*}}}*/ - - /*Avoid singularity at 0: */ - G_rigid[0]=G_rigid[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); - G_elastic[0]=G_elastic[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M)); - U_elastic[0]=U_elastic[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M)); - H_elastic[0]=H_elastic[1]; - parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M)); - - /*free ressources: */ - xDelete(love_h); - xDelete(love_k); - xDelete(love_l); - xDelete(love_th); - xDelete(love_tk); - xDelete(love_tl); - xDelete(G_rigid); - xDelete(G_rigid_local); - xDelete(G_elastic); - xDelete(G_elastic_local); - xDelete(U_elastic); - xDelete(U_elastic_local); - xDelete(H_elastic); - xDelete(H_elastic_local); + } + + /*merge G_elastic_local into G_elastic; U_elastic_local into U_elastic; H_elastic_local to H_elastic:{{{*/ + int* recvcounts=xNew(IssmComm::GetSize()); + int* displs=xNew(IssmComm::GetSize()); + + //recvcounts: + ISSM_MPI_Allgather(&m,1,ISSM_MPI_INT,recvcounts,1,ISSM_MPI_INT,IssmComm::GetComm()); + + /*displs: */ + ISSM_MPI_Allgather(&lower_row,1,ISSM_MPI_INT,displs,1,ISSM_MPI_INT,IssmComm::GetComm()); + + /*All gather:*/ + ISSM_MPI_Allgatherv(G_rigid_local, m, ISSM_MPI_DOUBLE, G_rigid, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + ISSM_MPI_Allgatherv(G_elastic_local, m, ISSM_MPI_DOUBLE, G_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + ISSM_MPI_Allgatherv(U_elastic_local, m, ISSM_MPI_DOUBLE, U_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + ISSM_MPI_Allgatherv(H_elastic_local, m, ISSM_MPI_DOUBLE, H_elastic, recvcounts, displs, ISSM_MPI_DOUBLE,IssmComm::GetComm()); + + /*free resources: */ + xDelete(recvcounts); + xDelete(displs); + + /*}}}*/ + + /*Avoid singularity at 0: */ + G_rigid[0]=G_rigid[1]; + parameters->AddObject(new DoubleVecParam(SealevelriseGRigidEnum,G_rigid,M)); + G_elastic[0]=G_elastic[1]; + parameters->AddObject(new DoubleVecParam(SealevelriseGElasticEnum,G_elastic,M)); + U_elastic[0]=U_elastic[1]; + parameters->AddObject(new DoubleVecParam(SealevelriseUElasticEnum,U_elastic,M)); + H_elastic[0]=H_elastic[1]; + parameters->AddObject(new DoubleVecParam(SealevelriseHElasticEnum,H_elastic,M)); + + /*free resources: */ + xDelete(love_h); + xDelete(love_k); + xDelete(love_l); + xDelete(love_th); + xDelete(love_tk); + xDelete(love_tl); + xDelete(G_rigid); + xDelete(G_rigid_local); + xDelete(G_elastic); + xDelete(G_elastic_local); + xDelete(U_elastic); + xDelete(U_elastic_local); + xDelete(H_elastic); + xDelete(H_elastic_local); + // } /*Indicate we have not yet ran the Geometry Core module: */ Index: /issm/trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 25626) +++ /issm/trunk-jpl/src/c/classes/Elements/Element.h (revision 25627) @@ -377,5 +377,7 @@ virtual IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks)=0; virtual void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks)=0; - virtual void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks,IssmDouble oceanarea)=0; + virtual IssmDouble SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks,IssmDouble oceanarea)=0; + virtual IssmDouble SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks,IssmDouble oceanarea)=0; + virtual void SealevelriseEustaticBottomPressure(IssmDouble* Sgi, SealevelMasks* masks)=0; virtual void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz)=0; virtual void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* mask)=0; Index: /issm/trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 25626) +++ /issm/trunk-jpl/src/c/classes/Elements/Penta.h (revision 25627) @@ -217,5 +217,7 @@ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks, IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks, IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelriseEustaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 25626) +++ /issm/trunk-jpl/src/c/classes/Elements/Seg.h (revision 25627) @@ -173,5 +173,7 @@ void SetSealevelMasks(SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks, IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks, IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelriseEustaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 25626) +++ /issm/trunk-jpl/src/c/classes/Elements/Tetra.h (revision 25627) @@ -179,5 +179,7 @@ void SealevelriseMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx, IssmDouble* yy, IssmDouble* zz){_error_("not implemented yet!");}; - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks,IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks, IssmDouble oceanarea){_error_("not implemented yet!");}; + IssmDouble SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks, IssmDouble oceanarea){_error_("not implemented yet!");}; + void SealevelriseEustaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseNonEustatic(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void SealevelriseGeodetic(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 25627) @@ -5468,4 +5468,10 @@ this->parameters->FindParam(&horiz,SolidearthSettingsHorizEnum); + /*compute area and add to inputs:*/ + area=GetAreaSpherical(); + this->inputs->SetDoubleInput(AreaEnum,this->lid,area); + this->AddInput(SealevelAreaEnum,&area,P0Enum); + if(!computerigid & !computeelastic)return; + /*recover precomputed green function kernels:*/ DoubleVecParam* parameter = static_cast(this->parameters->FindParamObject(SealevelriseGRigidEnum)); _assert_(parameter); @@ -5475,9 +5481,9 @@ parameter->GetParameterValueByPointer((IssmDouble**)&G_elastic_precomputed,&M); + parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); + parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); + parameter = static_cast(this->parameters->FindParamObject(SealevelriseHElasticEnum)); _assert_(parameter); parameter->GetParameterValueByPointer((IssmDouble**)&H_elastic_precomputed,&M); - - parameter = static_cast(this->parameters->FindParamObject(SealevelriseUElasticEnum)); _assert_(parameter); - parameter->GetParameterValueByPointer((IssmDouble**)&U_elastic_precomputed,&M); /*allocate indices:*/ @@ -5490,8 +5496,5 @@ } - /*compute area:*/ - area=GetAreaSpherical(); - - + /* Where is the centroid of this element:*/ @@ -5558,4 +5561,6 @@ this->inputs->SetArrayInput(SealevelriseGEEnum,this->lid,GE,gsize); } + + /* NOTE: May need to remove these two lines (from merge) */ this->inputs->SetDoubleInput(AreaEnum,this->lid,area); this->AddInput(SealevelAreaEnum,&area,P0Enum); @@ -5583,22 +5588,5 @@ } /*}}}*/ -void Tria::SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble oceanarea){ /*{{{*/ - - int bp_compute_fingerprints= 0; - - /*Compute bottom pressure contribution from ocean if requested:*/ - this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum); - if(bp_compute_fingerprints)this->SealevelriseBottomPressure(Sgi,masks); - - /*Compute eustatic ice contribution to sea level rise: */ - this->SealevelriseEustaticIce(Sgi,peustatic,masks, oceanarea); - - /*Compute hydrological contribution to sea level rise: */ - this->SealevelriseEustaticHydro(Sgi,peustatic,masks, oceanarea); - - -} -/*}}}*/ -void Tria::SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble oceanarea){ /*{{{*/ +IssmDouble Tria::SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks, IssmDouble oceanarea){ /*{{{*/ /*diverse:*/ @@ -5609,4 +5597,8 @@ bool notfullygrounded=false; bool scaleoceanarea= false; + bool computeelastic= true; + + /*output: */ + IssmDouble bslrice=0; /*elastic green function:*/ @@ -5618,7 +5610,4 @@ /*constants:*/ IssmDouble constant=0; - - /*Initialize eustatic component: do not skip this step :):*/ - IssmDouble eustatic = 0.; /*early return if we are not on an ice cap:*/ @@ -5627,6 +5616,6 @@ constant=0; this->AddInput(SealevelEustaticMaskEnum,&constant,P0Enum); #endif - *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! - return; + bslrice=0; + return bslrice; } @@ -5637,6 +5626,6 @@ this->AddInput(SealevelEustaticMaskEnum,&constant,P0Enum); #endif - *peustatic=0; //do not forget to assign this pointer, otherwise, global eustatic will be garbage! - return; + bslrice=0; + return bslrice; } @@ -5650,13 +5639,14 @@ #endif - /*recover material parameters: */ + /*recover some parameters:*/ rho_ice=FindParam(MaterialsRhoIceEnum); rho_water=FindParam(MaterialsRhoSeawaterEnum); - - /*recover ocean area scaling: */ + this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); /*retrieve precomputed G:*/ - this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); + + /* NOTE: May need to remove condition here (from merge) */ + if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ @@ -5702,30 +5692,32 @@ /*}}}*/ - /*Compute eustatic component:*/ + /*Compute barystatic contribution:*/ _assert_(oceanarea>0.); if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 - eustatic += rho_ice*area*phi*I/(oceanarea*rho_water); - - /*convert from m to kg/m^2:*/ - I=I*rho_ice*phi; - - /*convolve:*/ - for(int i=0;i(eustatic)); - *peustatic=eustatic; - return; -} -/*}}}*/ -void Tria::SealevelriseEustaticHydro(IssmDouble* Sgi, IssmDouble* peustatic, SealevelMasks* masks, IssmDouble oceanarea){ /*{{{*/ + bslrice = rho_ice*area*phi*I/(oceanarea*rho_water); + _assert_(!xIsNan(bslrice)); + + if(computeelastic){ + /*convert from m to kg/m^2:*/ + I=I*rho_ice*phi; + + /*convolve:*/ + for(int i=0;iisiceonly[this->lid]){ return; } + if(masks->isiceonly[this->lid]){ + bslrhydro=0; + return bslrhydro; + } /*early return if we are fully floating:*/ - if (masks->isfullyfloating[this->lid]){ return; } + if (masks->isfullyfloating[this->lid]){ + bslrhydro=0; + return bslrhydro; + } /*If we are here, we are on earth, not on ice: */ - /*recover material parameters: */ + /*recover parameters: */ rho_water=FindParam(MaterialsRhoSeawaterEnum); rho_freshwater=FindParam(MaterialsRhoFreshwaterEnum); - - /*recover ocean area scaling: */ + this->parameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); /*retrieve precomputed G:*/ - this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); + + /* NOTE: May need to remove condition here (from merge) */ + if(computeelastic)this->inputs->GetArrayPtr(SealevelriseGEnum,this->lid,&G,&gsize); /*Get area of element: precomputed in the sealevelrise_core_geometry:*/ @@ -5766,22 +5764,23 @@ deltathickness_input->GetInputAverage(&W); - /*Compute eustatic component:*/ + /*Compute barystatic component:*/ _assert_(oceanarea>0.); if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 - eustatic += rho_freshwater*area*phi*W/(oceanarea*rho_water); - - /*convert from m to kg/m^2:*/ - W=W*rho_freshwater*phi; - - /*convolve:*/ - for(int i=0;i(eustatic)); - *peustatic=eustatic; //do not forget the +, otherwise, you'll wipe out previous eustatic contributions! - return; -} -/*}}}*/ -void Tria::SealevelriseBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/ + bslrhydro = rho_freshwater*area*phi*W/(oceanarea*rho_water); + _assert_(!xIsNan(bslrhydro)); + + if(computeelastic){ + /*convert from m to kg/m^2:*/ + W=W*rho_freshwater*phi; + + /*convolve:*/ + for(int i=0;iparameters->FindParam(&computeelastic,SolidearthSettingsElasticEnum); + if(!computeelastic)return; /*we are here to compute fingerprints originating fromn bottom pressure changes:*/ Index: /issm/trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 25626) +++ /issm/trunk-jpl/src/c/classes/Elements/Tria.h (revision 25627) @@ -165,8 +165,7 @@ void SetSealevelMasks(SealevelMasks* masks); void SealevelriseGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz); - void SealevelriseEustatic(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble oceanarea); - void SealevelriseEustaticIce(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble oceanarea); - void SealevelriseEustaticHydro(IssmDouble* Sgi, IssmDouble* peustatic,SealevelMasks* masks, IssmDouble oceanarea); - void SealevelriseBottomPressure(IssmDouble* Sgi,SealevelMasks* masks); + IssmDouble SealevelriseEustaticIce(IssmDouble* Sgi, SealevelMasks* masks, IssmDouble oceanarea); + IssmDouble SealevelriseEustaticHydro(IssmDouble* Sgi, SealevelMasks* masks, IssmDouble oceanarea); + void SealevelriseEustaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks); void SealevelriseNonEustatic(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks); void SealevelriseGeodetic(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks); Index: /issm/trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/classes/FemModel.cpp (revision 25627) @@ -53,4 +53,7 @@ /*Object constructors and destructor*/ +FemModel::FemModel(void){ /*{{{*/ + /*do nothing:*/ +} /*}}}*/ FemModel::FemModel(int argc,char** argv,ISSM_MPI_Comm incomm,bool trace){/*{{{*/ @@ -66,4 +69,5 @@ char *restartfilename = NULL; char *rootpath = NULL; + char *modelname = NULL; /*First things first, store the communicator, and set it as a global variable: */ @@ -81,9 +85,9 @@ /*From command line arguments, retrieve different filenames needed to create the FemModel: */ - ProcessArguments(&solution_type,&binfilename,&outbinfilename,&petscfilename,&lockfilename,&restartfilename,&rootpath,argc,argv); + ProcessArguments(&solution_type,&binfilename,&outbinfilename,&petscfilename,&lockfilename,&restartfilename,&rootpath,&modelname,argc,argv); /*Create femmodel from input files: */ profiler->Start(MPROCESSOR); - this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,restartfilename, solution_type,trace,NULL); + this->InitFromFiles(rootpath,binfilename,outbinfilename,petscfilename,lockfilename,restartfilename, modelname, solution_type,trace,NULL); profiler->Stop(MPROCESSOR); @@ -149,5 +153,5 @@ } /*}}}*/ -FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X){ /*{{{*/ +FemModel::FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, char* modelname, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X){ /*{{{*/ bool traceon=true; @@ -160,5 +164,5 @@ /*Create femmodel from input files, with trace activated: */ profiler->Start(MPROCESSOR); - this->InitFromFiles(rootpath,inputfilename,outputfilename,toolkitsfilename,lockfilename,restartfilename, solution_type,traceon,X); + this->InitFromFiles(rootpath,inputfilename,outputfilename,toolkitsfilename,lockfilename,restartfilename, modelname,solution_type,traceon,X); profiler->Stop(MPROCESSOR); @@ -418,5 +422,5 @@ } /*}}}*/ -void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, const int in_solution_type,bool trace,IssmPDouble* X){/*{{{*/ +void FemModel::InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, char* modelname, const int in_solution_type,bool trace,IssmPDouble* X){/*{{{*/ /*intermediary*/ @@ -442,4 +446,5 @@ /*Now save all of these file names into parameters, you never know when you might need them: */ this->parameters->AddObject(new StringParam(ToolkitsFileNameEnum,toolkitsfilename)); + this->parameters->AddObject(new StringParam(ModelnameEnum,modelname)); this->parameters->AddObject(new StringParam(RootPathEnum,rootpath)); this->parameters->AddObject(new StringParam(InputFileNameEnum,inputfilename)); @@ -4742,13 +4747,15 @@ #endif #ifdef _HAVE_SEALEVELRISE_ -void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks) { /*{{{*/ +void FemModel::SealevelriseEustatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslr,IssmDouble* pbslrice, IssmDouble* pbslrhydro, SealevelMasks* masks) { /*{{{*/ /*serialized vectors:*/ - IssmDouble eustatic = 0.; - IssmDouble eustatic_cpu = 0.; - IssmDouble eustatic_cpu_e = 0.; + IssmDouble bslrice = 0.; + IssmDouble bslrice_cpu = 0.; + IssmDouble bslrhydro = 0.; + IssmDouble bslrhydro_cpu = 0.; IssmDouble area = 0.; IssmDouble oceanarea = 0.; IssmDouble oceanarea_cpu = 0.; + int bp_compute_fingerprints= 0; /*Initialize temporary vector that will be used to sum eustatic components @@ -4771,9 +4778,25 @@ _assert_(oceanarea>0.); - /*Call the sea level rise core: */ + /*Call the sea level rise core for ice : */ + bslrice_cpu=0; for(Object* & object : this->elements->objects){ Element* element = xDynamicCast(object); - element->SealevelriseEustatic(RSLgi,&eustatic_cpu_e,masks, oceanarea); - eustatic_cpu+=eustatic_cpu_e; + bslrice_cpu+=element->SealevelriseEustaticIce(RSLgi,masks, oceanarea); + } + + /*Call the sea level rise core for hydro: */ + bslrhydro_cpu=0; + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + bslrhydro_cpu+=element->SealevelriseEustaticHydro(RSLgi,masks, oceanarea); + } + + /*Call the sea level rise core for bottom pressures: */ + this->parameters->FindParam(&bp_compute_fingerprints,DslComputeFingerprintsEnum); + if(bp_compute_fingerprints){ + for(int i=0;iSize();i++){ + Element* element=xDynamicCast(elements->GetObjectByOffset(i)); + element->SealevelriseEustaticBottomPressure(RSLgi,masks); + } } @@ -4783,7 +4806,11 @@ /*Sum all eustatic components from all cpus:*/ - ISSM_MPI_Reduce (&eustatic_cpu,&eustatic,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&eustatic,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - _assert_(!xIsNan(eustatic)); + ISSM_MPI_Reduce (&bslrice_cpu,&bslrice,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&bslrice,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + _assert_(!xIsNan(bslrice)); + + ISSM_MPI_Reduce (&bslrhydro_cpu,&bslrhydro,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); + ISSM_MPI_Bcast(&bslrhydro,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); + _assert_(!xIsNan(bslrhydro)); /*Free ressources:*/ @@ -4793,5 +4820,7 @@ /*Assign output pointers:*/ *poceanarea = oceanarea; - *peustatic = eustatic; + *pbslrice = bslrice; + *pbslrhydro = bslrhydro; + *pbslr=bslrhydro+bslrice; } Index: /issm/trunk-jpl/src/c/classes/FemModel.h =================================================================== --- /issm/trunk-jpl/src/c/classes/FemModel.h (revision 25626) +++ /issm/trunk-jpl/src/c/classes/FemModel.h (revision 25627) @@ -67,6 +67,7 @@ /*constructors, destructors: */ + FemModel(void); FemModel(int argc,char** argv,ISSM_MPI_Comm comm_init,bool trace=false); - FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X); + FemModel(char* rootpath, char* inputfilename, char* outputfilename, char* toolkitsfilename, char* lockfilename, char* restartfilename, char* modelname, ISSM_MPI_Comm incomm, int solution_type,IssmPDouble* X); ~FemModel(); @@ -79,5 +80,5 @@ void Echo(); int GetElementsWidth(){return 3;};//just tria elements in this first version - void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, char* restartfilename, const int solution_type,bool trace,IssmPDouble* X=NULL); + void InitFromFiles(char* rootpath, char* inputfilename, char* outputfilename, char* petscfilename, char* lockfilename, char* restartfilename, char* modelname, const int solution_type,bool trace,IssmPDouble* X=NULL); void InitFromFids(char* rootpath, FILE* IOMODEL, FILE* toolkitsoptionsfid, int in_solution_type, bool trace, IssmPDouble* X=NULL); void Marshall(MarshallHandle* marshallhandle); @@ -165,5 +166,5 @@ #endif #ifdef _HAVE_SEALEVELRISE_ - void SealevelriseEustatic(Vector* pSgi, IssmDouble* poceanarea, IssmDouble* peustatic, SealevelMasks* masks); + void SealevelriseEustatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslr,IssmDouble* pbslrice, IssmDouble* pbslrhydro, SealevelMasks* masks); void SealevelriseNonEustatic(Vector* pSgo, Vector* pSg_old, SealevelMasks* masks,bool verboseconvolution); void SealevelriseRotationalFeedback(Vector* pRSLgo_rot, Vector* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks); Index: /issm/trunk-jpl/src/c/cores/ProcessArguments.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/ProcessArguments.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/cores/ProcessArguments.cpp (revision 25627) @@ -8,5 +8,5 @@ #include "../shared/shared.h" -void ProcessArguments(int* solution_type,char** pbinfilename,char** poutbinfilename,char** ptoolkitsfilename,char** plockfilename,char** prestartfilename, char** prootpath, int argc,char **argv){ +void ProcessArguments(int* solution_type,char** pbinfilename,char** poutbinfilename,char** ptoolkitsfilename,char** plockfilename,char** prestartfilename, char** prootpath, char** pmodelname, int argc,char **argv){ /*Check input arguments*/ @@ -18,5 +18,6 @@ *solution_type = StringToEnumx(argv[1]); char* rootpatharg = argv[2]; - char* modelname = argv[3]; + char* modelname = xNew(strlen(argv[3])+1); + xMemCpy(modelname,argv[3],strlen(argv[3])+1); /*Recover myrank and length of string "my_rank" */ @@ -42,4 +43,5 @@ *prestartfilename=restartfilename; *prootpath=rootpath; + *pmodelname=modelname; } Index: /issm/trunk-jpl/src/c/cores/cores.h =================================================================== --- /issm/trunk-jpl/src/c/cores/cores.h (revision 25626) +++ /issm/trunk-jpl/src/c/cores/cores.h (revision 25627) @@ -73,5 +73,5 @@ //diverse -void ProcessArguments(int* solution,char** pbinname,char** poutbinname,char** ptoolkitsname,char** plockname,char** prestartname, char** prootpath,int argc,char **argv); +void ProcessArguments(int* solution,char** pbinname,char** poutbinname,char** ptoolkitsname,char** plockname,char** prestartname, char** prootpath,char** pmodelname, int argc,char **argv); void WriteLockFile(char* filename); void ResetBoundaryConditions(FemModel* femmodel, int analysis_type); Index: /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 25627) @@ -268,4 +268,10 @@ IssmDouble dt; int geodetic=0; + int step; + IssmDouble time; + IssmDouble cumgmtslr=0; + IssmDouble cumbslr=0; + IssmDouble cumgmslr=0; + IssmDouble gmtslr=0; /*Retrieve parameters:*/ @@ -274,4 +280,8 @@ femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum); femmodel->parameters->FindParam(&geodetic,SolidearthSettingsComputesealevelchangeEnum); + femmodel->parameters->FindParam(&cumgmtslr,CumGmtslrEnum); + femmodel->parameters->FindParam(&cumbslr,CumBslrEnum); + femmodel->parameters->FindParam(&step,StepEnum); + femmodel->parameters->FindParam(&time,TimeEnum); /*in case we are running SealevelriseSolutionEnum, then bypass transient settings:*/ @@ -296,4 +306,13 @@ } + /*Save steric rate, and cumulate it through time, cumulate with barystatic sea level rise, and cumulate through time:*/ + gmtslr=steric_rate_g->Max()*dt; + cumgmtslr+=gmtslr; + cumgmslr=cumbslr+cumgmtslr; + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,GmtslrEnum,gmtslr,step,time)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumGmtslrEnum,cumgmtslr,step,time)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumGmslrEnum,cumgmslr,step,time)); + femmodel->parameters->SetParam(cumgmtslr,CumGmtslrEnum); + /*compute: sea level change = initial sea level + (N_gia_rate+N_esa_rate) * dt + steric_rate + dynamic_rate dt*/ if(geodetic){ @@ -376,5 +395,5 @@ - /*Run sealevelrie geometry routine in elements:*/ + /*Run sealevelrise geometry routine in elements:*/ for(Object* & object : femmodel->elements->objects){ Element* element=xDynamicCast(object); @@ -410,10 +429,14 @@ IssmDouble time; - /*outputs:*/ - IssmDouble eustatic; + /*barystatic contribution:*/ + IssmDouble bslr; + IssmDouble bslrice; + IssmDouble bslrhydro; + IssmDouble cumbslr; + IssmDouble cumbslrice; + IssmDouble cumbslrhydro; - if(VerboseSolution()) _printf0_(" computing eustatic components on ice\n"); - - + if(VerboseSolution()) _printf0_(" computing bslr components on ice\n"); + /*Figure out size of g-set deflection vector and allocate solution vector: */ gsize = femmodel->nodes->NumberOfDofs(GsetEnum); @@ -422,21 +445,38 @@ femmodel->parameters->FindParam(&step,StepEnum); femmodel->parameters->FindParam(&time,TimeEnum); + femmodel->parameters->FindParam(&cumbslr,CumBslrEnum); + femmodel->parameters->FindParam(&cumbslrice,CumBslrIceEnum); + femmodel->parameters->FindParam(&cumbslrhydro,CumBslrHydroEnum); /*Initialize:*/ RSLgi = new Vector(gsize); - /*call the eustatic main module: */ - femmodel->SealevelriseEustatic(RSLgi,&oceanarea,&eustatic, masks); //this computes + /*call the bslr main module: */ + femmodel->SealevelriseEustatic(RSLgi,&oceanarea,&bslr, &bslrice, &bslrhydro, masks); //this computes /*we need to average RSLgi over the ocean: RHS term 4 in Eq.4 of Farrel and clarke. Only the elements can do that: */ RSLgi_oceanaverage=femmodel->SealevelriseOceanAverage(RSLgi,masks, oceanarea); - /*RSLg is the sum of the pure eustatic component (term 3) and the contribution from the perturbation to the graviation potential due to the + /*RSLg is the sum of the pure bslr component (term 3) and the contribution from the perturbation to the graviation potential due to the * presence of ice (terms 1 and 4 in Eq.4 of Farrel and Clarke):*/ - RSLgi->Shift(-eustatic-RSLgi_oceanaverage); - - /*save eustatic value for results: */ - femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,SealevelRSLEustaticEnum,-eustatic,step,time)); - + RSLgi->Shift(-bslr-RSLgi_oceanaverage); + + /*save bslr and cumulated bslr value for results:{{{ */ + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,BslrEnum,-bslr,step,time)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,BslrIceEnum,-bslrice,step,time)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,BslrHydroEnum,-bslrhydro,step,time)); + + cumbslr=cumbslr-bslr; + cumbslrice=cumbslrice-bslrice; + cumbslrhydro=cumbslrhydro-bslrhydro; + femmodel->parameters->SetParam(cumbslr,CumBslrEnum); + femmodel->parameters->SetParam(cumbslrice,CumBslrIceEnum); + femmodel->parameters->SetParam(cumbslrhydro,CumBslrHydroEnum); + + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumBslrEnum,cumbslr,step,time)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumBslrIceEnum,cumbslrice,step,time)); + femmodel->results->AddResult(new GenericExternalResult(femmodel->results->Size()+1,CumBslrHydroEnum,cumbslrhydro,step,time)); + /*}}}*/ + /*Assign output pointers and return: */ *poceanarea=oceanarea; @@ -465,5 +505,4 @@ IssmDouble eps_rel; IssmDouble eps_abs; - IssmDouble eustatic; IssmDouble Ixz, Iyz, Izz; Index: /issm/trunk-jpl/src/c/issm_post =================================================================== --- /issm/trunk-jpl/src/c/issm_post (revision 25627) +++ /issm/trunk-jpl/src/c/issm_post (revision 25627) @@ -0,0 +1,210 @@ +#! /bin/sh + +# issm_post - temporary wrapper script for .libs/issm_post.exe +# Generated by libtool (GNU libtool) 2.4.2 +# +# The issm_post program cannot be directly executed until all the libtool +# libraries that it depends on are installed. +# +# This wrapper script should never be moved out of the build directory. +# If it is, it will not operate correctly. + +# Sed substitution that helps us do robust quoting. It backslashifies +# metacharacters that are still active within double-quoted strings. +sed_quote_subst='s/\([`"$\\]\)/\\\1/g' + +# Be Bourne compatible +if test -n "${ZSH_VERSION+set}" && (emulate sh) >/dev/null 2>&1; then + emulate sh + NULLCMD=: + # Zsh 3.x and 4.x performs word splitting on ${1+"$@"}, which + # is contrary to our usage. Disable this feature. + alias -g '${1+"$@"}'='"$@"' + setopt NO_GLOB_SUBST +else + case `(set -o) 2>/dev/null` in *posix*) set -o posix;; esac +fi +BIN_SH=xpg4; export BIN_SH # for Tru64 +DUALCASE=1; export DUALCASE # for MKS sh + +# The HP-UX ksh and POSIX shell print the target directory to stdout +# if CDPATH is set. +(unset CDPATH) >/dev/null 2>&1 && unset CDPATH + +relink_command="" + +# This environment variable determines our operation mode. +if test "$libtool_install_magic" = "%%%MAGIC variable%%%"; then + # install mode needs the following variables: + generated_by_libtool_version='2.4.2' + notinst_deplibs=' ./libISSMCore.la ./libISSMOverload.la ./libISSMModules.la /Users/jdq/Repos/issm/issm/trunk-jpl-new/src/c/libISSMCore.la' +else + # When we are sourced in execute mode, $file and $ECHO are already set. + if test "$libtool_execute_magic" != "%%%MAGIC variable%%%"; then + file="$0" + +# A function that is used when there is no print builtin or printf. +func_fallback_echo () +{ + eval 'cat <<_LTECHO_EOF +$1 +_LTECHO_EOF' +} + ECHO="printf %s\\n" + fi + +# Very basic option parsing. These options are (a) specific to +# the libtool wrapper, (b) are identical between the wrapper +# /script/ and the wrapper /executable/ which is used only on +# windows platforms, and (c) all begin with the string --lt- +# (application programs are unlikely to have options which match +# this pattern). +# +# There are only two supported options: --lt-debug and +# --lt-dump-script. There is, deliberately, no --lt-help. +# +# The first argument to this parsing function should be the +# script's ../../libtool value, followed by no. +lt_option_debug= +func_parse_lt_options () +{ + lt_script_arg0=$0 + shift + for lt_opt + do + case "$lt_opt" in + --lt-debug) lt_option_debug=1 ;; + --lt-dump-script) + lt_dump_D=`$ECHO "X$lt_script_arg0" | /usr/bin/sed -e 's/^X//' -e 's%/[^/]*$%%'` + test "X$lt_dump_D" = "X$lt_script_arg0" && lt_dump_D=. + lt_dump_F=`$ECHO "X$lt_script_arg0" | /usr/bin/sed -e 's/^X//' -e 's%^.*/%%'` + cat "$lt_dump_D/$lt_dump_F" + exit 0 + ;; + --lt-*) + $ECHO "Unrecognized --lt- option: '$lt_opt'" 1>&2 + exit 1 + ;; + esac + done + + # Print the debug banner immediately: + if test -n "$lt_option_debug"; then + echo "issm_post.exe:issm_post:${LINENO}: libtool wrapper (GNU libtool) 2.4.2" 1>&2 + fi +} + +# Used when --lt-debug. Prints its arguments to stdout +# (redirection is the responsibility of the caller) +func_lt_dump_args () +{ + lt_dump_args_N=1; + for lt_arg + do + $ECHO "issm_post.exe:issm_post:${LINENO}: newargv[$lt_dump_args_N]: $lt_arg" + lt_dump_args_N=`expr $lt_dump_args_N + 1` + done +} + +# Core function for launching the target application +func_exec_program_core () +{ + + if test -n "$lt_option_debug"; then + $ECHO "issm_post.exe:issm_post:${LINENO}: newargv[0]: $progdir/$program" 1>&2 + func_lt_dump_args ${1+"$@"} 1>&2 + fi + exec "$progdir/$program" ${1+"$@"} + + $ECHO "$0: cannot exec $program $*" 1>&2 + exit 1 +} + +# A function to encapsulate launching the target application +# Strips options in the --lt-* namespace from $@ and +# launches target application with the remaining arguments. +func_exec_program () +{ + case " $* " in + *\ --lt-*) + for lt_wr_arg + do + case $lt_wr_arg in + --lt-*) ;; + *) set x "$@" "$lt_wr_arg"; shift;; + esac + shift + done ;; + esac + func_exec_program_core ${1+"$@"} +} + + # Parse options + func_parse_lt_options "$0" ${1+"$@"} + + # Find the directory that this script lives in. + thisdir=`$ECHO "$file" | /usr/bin/sed 's%/[^/]*$%%'` + test "x$thisdir" = "x$file" && thisdir=. + + # Follow symbolic links until we get to the real thisdir. + file=`ls -ld "$file" | /usr/bin/sed -n 's/.*-> //p'` + while test -n "$file"; do + destdir=`$ECHO "$file" | /usr/bin/sed 's%/[^/]*$%%'` + + # If there was a directory component, then change thisdir. + if test "x$destdir" != "x$file"; then + case "$destdir" in + [\\/]* | [A-Za-z]:[\\/]*) thisdir="$destdir" ;; + *) thisdir="$thisdir/$destdir" ;; + esac + fi + + file=`$ECHO "$file" | /usr/bin/sed 's%^.*/%%'` + file=`ls -ld "$thisdir/$file" | /usr/bin/sed -n 's/.*-> //p'` + done + + # Usually 'no', except on cygwin/mingw when embedded into + # the cwrapper. + WRAPPER_SCRIPT_BELONGS_IN_OBJDIR=no + if test "$WRAPPER_SCRIPT_BELONGS_IN_OBJDIR" = "yes"; then + # special case for '.' + if test "$thisdir" = "."; then + thisdir=`pwd` + fi + # remove .libs from thisdir + case "$thisdir" in + *[\\/].libs ) thisdir=`$ECHO "$thisdir" | /usr/bin/sed 's%[\\/][^\\/]*$%%'` ;; + .libs ) thisdir=. ;; + esac + fi + + # Try to get the absolute directory name. + absdir=`cd "$thisdir" && pwd` + test -n "$absdir" && thisdir="$absdir" + + program='issm_post.exe' + progdir="$thisdir/.libs" + + + if test -f "$progdir/$program"; then + # Add our own library path to DYLD_LIBRARY_PATH + DYLD_LIBRARY_PATH="/Users/jdq/Repos/issm/issm/trunk-jpl-new/src/c/.libs:$DYLD_LIBRARY_PATH" + + # Some systems cannot cope with colon-terminated DYLD_LIBRARY_PATH + # The second colon is a workaround for a bug in BeOS R4 sed + DYLD_LIBRARY_PATH=`$ECHO "$DYLD_LIBRARY_PATH" | /usr/bin/sed 's/::*$//'` + + export DYLD_LIBRARY_PATH + + if test "$libtool_execute_magic" != "%%%MAGIC variable%%%"; then + # Run the actual program with our arguments. + func_exec_program ${1+"$@"} + fi + else + # The program doesn't exist. + $ECHO "$0: error: \`$progdir/$program' does not exist" 1>&2 + $ECHO "This script is just a wrapper for $program." 1>&2 + $ECHO "See the libtool documentation for more information." 1>&2 + exit 1 + fi +fi Index: /issm/trunk-jpl/src/c/main/issm_dakota.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/issm_dakota.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/main/issm_dakota.cpp (revision 25627) @@ -14,5 +14,5 @@ #endif -int main(int argc,char **argv){ +int main(int argc,char **argv){ /*{{{*/ #if defined(_HAVE_DAKOTA_) && _DAKOTA_MAJOR_ >= 6 @@ -22,4 +22,5 @@ char* dakota_output_file = NULL; char* dakota_error_file = NULL; + bool statistics=false; /*Define MPI_DEBUG in dakota_global_defs.cpp to cause a hold here*/ @@ -38,4 +39,7 @@ dakota_error_file=xNew((strlen(argv[2])+strlen(argv[3])+strlen(".qmu.err")+2)); sprintf(dakota_error_file,"%s/%s%s",argv[2],argv[3],".qmu.err"); + + /*Create directory structure for model outputs:*/ + statistics=DakotaDirStructure(argc,argv); /* Parse input and construct Dakota LibraryEnvironment, performing input data checks*/ @@ -83,4 +87,8 @@ env.execute(); + /* Run statistics if requested:*/ + if(statistics)DakotaStatistics(argc,argv); + + /*free allocations:*/ xDelete(dakota_input_file); xDelete(dakota_output_file); @@ -94,3 +102,3 @@ #endif -} +} /*}}}*/ Index: /issm/trunk-jpl/src/c/main/issm_post.cpp =================================================================== --- /issm/trunk-jpl/src/c/main/issm_post.cpp (revision 25627) +++ /issm/trunk-jpl/src/c/main/issm_post.cpp (revision 25627) @@ -0,0 +1,24 @@ +/*!\file: issm_post.cpp + * \brief: ISSM DAKOTA post-processing of statistics + */ + +#include "./issm.h" +#include + +int main(int argc,char **argv){ /*{{{*/ + + char* dakota_input_file=NULL; + char* dakota_output_file = NULL; + char* dakota_error_file = NULL; + bool statistics=false; + + /*Initialize MPI: */ + ISSM_MPI_Init(&argc, &argv); // initialize MPI + + /*Run statistics:*/ + DakotaStatistics(argc,argv); + + /*Return unix success: */ + return 0; + +} /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (revision 25627) @@ -8,5 +8,5 @@ #include "./ModelProcessorx.h" -#define MAXCONNECTIVITY 5 +#define MAXCONNECTIVITY 15 bool IsVertexInRank(int* vertices_ranks,int* vertices_proc_count,int vid,int rank){/*{{{*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateNodes.cpp (revision 25627) @@ -9,5 +9,5 @@ #include "./ModelProcessorx.h" -#define MAXCONNECTIVITY 5 +#define MAXCONNECTIVITY 15 bool IsNodeInRank(int* nodes_ranks,int* nodes_proc_count,int nid,int rank){/*{{{*/ Index: /issm/trunk-jpl/src/c/modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/modules/ModelProcessorx/Dakota/CreateParametersDakota.cpp (revision 25627) @@ -41,4 +41,9 @@ int M,N; + //qmu statistics + bool statistics = false; + int numdirectories = 0; + int nfilesperdirectory = 0; + /*recover parameters: */ iomodel->FindConstant(&dakota_analysis,"md.qmu.isdakota"); @@ -72,4 +77,15 @@ /*Ok, we have all the response descriptors. Build a parameter with it: */ parameters->AddObject(new StringArrayParam(QmuResponsedescriptorsEnum,responsedescriptors,numresponsedescriptors)); + + /*Deal with statistics: */ + iomodel->FindConstant(&statistics,"md.qmu.statistics"); + parameters->AddObject(new BoolParam(QmuStatisticsEnum,statistics)); + if(statistics){ + iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories"); + parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories)); + + iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory"); + parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory)); + } /*Load partitioning vectors specific to variables:*/ @@ -109,5 +125,4 @@ /*}}}*/ - /*Deal with data needed because of qmu variables*/ DataSet* dataset_variable_descriptors = new DataSet(QmuVariableDescriptorsEnum); @@ -143,5 +158,5 @@ delete dataset_variable_descriptors; - /*clean-up*/ + /*clean-up {{{*/ for(i=0;i(qmuerrname); xDelete(qmuoutname); - - - + xDelete(name); + /*}}}*/ } - /*Free data*/ - xDelete(name); } Index: /issm/trunk-jpl/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/modules/NodeConnectivityx/NodeConnectivityx.cpp (revision 25627) @@ -19,5 +19,5 @@ int i,j,n; - const int maxels=50; + const int maxels=100; const int width=maxels+1; Index: /issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/modules/OutputResultsx/OutputResultsx.cpp (revision 25627) @@ -63,7 +63,25 @@ else{ if(my_rank==0){ + /*a little bit complicated. Either statistic computations are requested, which means we + * put our outbin files in subidirectories with numbers, or we don't, and we dump our + * outbins directly in the current directory:*/ int currEvalId ; + int nfilesperdirectory; + bool statistics=false; + char* root=NULL; + char* modelname=NULL; + femmodel->parameters->FindParam(&currEvalId,QmuCurrEvalIdEnum); - sprintf(outputfilename2,"%s.%i",outputfilename,currEvalId); + femmodel->parameters->FindParam(&statistics,QmuStatisticsEnum); + + if(statistics){ + femmodel->parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum); + femmodel->parameters->FindParam(&root,RootPathEnum); + femmodel->parameters->FindParam(&modelname,ModelnameEnum); + sprintf(outputfilename2,"%s/%i/%s.outbin.%i",root,(int)(floor((currEvalId-1)/nfilesperdirectory)+1),modelname,currEvalId); + } + else{ + sprintf(outputfilename2,"%s.%i",outputfilename,currEvalId); + } fid=pfopen0(outputfilename2,"ab+"); } Index: /issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp =================================================================== --- /issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp (revision 25627) +++ /issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.cpp (revision 25627) @@ -0,0 +1,1407 @@ +/*!\file: QmuStatisticsx routines + */ +/*includes and prototypes:*/ +#include +#include "./QmuStatisticsx.h" +#include "../OutputResultsx/OutputResultsx.h" + +int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step){ /*{{{*/ + + int length; + char fieldname[1000]; + int fieldname_size; + IssmDouble rtime; + int rstep; + int M,N; + + //fields that we retrive: + IssmDouble dfield; + char* sfield = NULL; + IssmDouble* dmatfield = NULL; + int* imatfield = NULL; + + //type of the returned field: + int type; + int found=0; + + while(1){ + + size_t ret_code = fread(&fieldname_size, sizeof(int), 1, fid); + if(ret_code != 1) break; //we are done. + + fread(fieldname, sizeof(char), fieldname_size, fid); + //_printf0_("fieldname: " << fieldname << "\n"); + + fread(&rtime, sizeof(IssmDouble), 1, fid); + fread(&rstep, sizeof(int), 1, fid); + + //check on field: + if ((step==rstep) && (strcmp(field,fieldname)==0)){ + + //ok, go read the result really: + fread(&type,sizeof(int),1,fid); + fread(&M,sizeof(int),1,fid); + if (type==1){ + fread(&dfield,sizeof(IssmDouble),1,fid); + } + else if (type==2){ + fread(&M,sizeof(int),1,fid); + sfield=xNew(M); + fread(sfield,sizeof(char),M,fid); + } + else if (type==3){ + fread(&N,sizeof(int),1,fid); + dmatfield=xNew(M*N); + fread(dmatfield,sizeof(IssmDouble),M*N,fid); + } + else if (type==4){ + fread(&N,sizeof(int),1,fid); + imatfield=xNew(M*N); + fread(imatfield,sizeof(int),M*N,fid); + } + else _error_("cannot read data of type " << type << "\n"); + found=1; + break; + } + else{ + //just skim to next results. + fread(&type,sizeof(int),1,fid); + fread(&M,sizeof(int),1,fid); + if (type==1){ + fseek(fid,sizeof(IssmDouble),SEEK_CUR); + } + else if(type==2){ + fseek(fid,M*sizeof(char),SEEK_CUR); + } + else if(type==3){ + fread(&N,sizeof(int),1,fid); + fseek(fid,M*N*sizeof(IssmDouble),SEEK_CUR); + } + else if(type==4){ + fread(&N,sizeof(int),1,fid); + fseek(fid,M*N*sizeof(int),SEEK_CUR); + } + else _error_("cannot read data of type " << type << "\n"); + } + } + if(found==0)_error_("cound not find " << field << " at step " << step << "\n"); + + /*assign output pointers:*/ + *pdoublemat=dmatfield; + *pdoublematsize=M*N; + *pdouble=dfield; + + /*return:*/ + return type; + +} +/*}}}*/ +int ComputeHistogram(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/ + + int nsamples; + char* directory=NULL; + char* model=NULL; + char** fields=NULL; + int* steps=NULL; + int nsteps; + int nfields; + int nbins; + int range,lower_row,upper_row; + int nfilesperdirectory; + + /*intermediary:*/ + IssmDouble* doublemat=NULL; + int doublematsize; + IssmDouble scalar; + + /*computation of average and variance itself:*/ + IssmDouble** maxxs = NULL; + IssmDouble** minxs = NULL; + int* xtype=NULL; + int* xsize=NULL; + + IssmDouble** maxmeans=NULL; + IssmDouble** minmeans=NULL; + int* meanxtype=NULL; + int* meanxsize=NULL; + + /*only work on the statistical communicator: */ + if (color==MPI_UNDEFINED)return 0; + + /*Retrieve parameters:*/ + parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum); + parameters->FindParam(&nsamples,QmuNsampleEnum); + parameters->FindParam(&directory,DirectoryNameEnum); + parameters->FindParam(&model,InputFileNameEnum); + parameters->FindParam(&fields,&nfields,FieldsEnum); + parameters->FindParam(&steps,&nsteps,StepsEnum); + parameters->FindParam(&nbins,NbinsEnum); + + /*Get rank from the stat comm communicator:*/ + IssmComm::SetComm(statcomm); + int my_rank=IssmComm::GetRank(); + + /*Open files and read them complelety, in a distributed way:*/ + range=DetermineLocalSize(nsamples,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm()); + + /*Initialize arrays:*/ + maxmeans=xNew(nfields); + minmeans=xNew(nfields); + meanxtype=xNew(nfields); + meanxsize=xNew(nfields); + + maxxs=xNew(nfields*nsteps); + minxs=xNew(nfields*nsteps); + xtype=xNew(nfields*nsteps); + xsize=xNew(nfields*nsteps); + + /*Start opening files:*/ + for(int i=(lower_row+1);i<=upper_row;i++){ + _printf0_("reading file #: " << i << "\n"); + char file[1000]; + long int length; + char* buffer=NULL; + + /*string:*/ + sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); + + /*open file: */ + _printf0_(" opening file: " << file << "\n"); + FILE* fid=fopen(file,"rb"); + if(fid==NULL)_error_("cound not open file: " << file << "\n"); + + /*figure out size of file, and read the whole thing:*/ + _printf0_(" reading file:\n"); + fseek (fid, 0, SEEK_END); + length = ftell (fid); + fseek (fid, 0, SEEK_SET); + buffer = xNew(length); + fread (buffer, sizeof(char), length, fid); + + /*close file:*/ + fclose (fid); + + /*create a memory stream with this buffer:*/ + _printf0_(" processing file:\n"); + fid=fmemopen(buffer, length, "rb"); + + /*start reading data from the buffer directly:*/ + for (int f=0;f(1); + minxs[counter]=xNew(1); + *maxxs[counter]=scalar; + *minxs[counter]=scalar; + xsize[counter]=1; + } + else if (xtype[counter]==3){ + maxxs[counter]=xNew(doublematsize); + xMemCpy(maxxs[counter],doublemat,doublematsize); + minxs[counter]=xNew(doublematsize); + xMemCpy(minxs[counter],doublemat,doublematsize); + xsize[counter]=doublematsize; + xDelete(doublemat); + } + else _error_("cannot carry out statistics on type " << xtype[counter]); + } + else{ + if(xtype[counter]==1){ + *maxxs[counter]=max(*maxxs[counter],scalar); + *minxs[counter]=min(*minxs[counter],scalar); + } + else if (xtype[counter]==3){ + IssmDouble* newmax=maxxs[counter]; + IssmDouble* newmin=minxs[counter]; + for(int k=0;knewmax[k])newmax[k]=doublemat[k]; + if(doublemat[k](doublemat); + } + else _error_("cannot carry out statistics on type " << xtype[counter]); + } + } + } + _printf0_(" average in time:\n"); + + /*Deal with average in time: */ + for (int f=0;f(1); + minmeans[f]=xNewZeroInit(1); + *maxmeans[f]=timemean; + *minmeans[f]=timemean; + } + else{ + *maxmeans[f]=max(*maxmeans[f],timemean); + *minmeans[f]=min(*minmeans[f],timemean); + } + } + else{ + meanxsize[f]=doublematsize; + fseek(fid,0,SEEK_SET); + IssmDouble* timemean=xNewZeroInit(doublematsize); + for (int j=0;j(doublemat); + } + + if(i==(lower_row+1)){ + maxmeans[f]=xNew(doublematsize); + xMemCpy(maxmeans[f],timemean,doublematsize); + minmeans[f]=xNew(doublematsize); + xMemCpy(minmeans[f],timemean,doublematsize); + } + else{ + IssmDouble* maxx=maxmeans[f]; + IssmDouble* minx=minmeans[f]; + + for(int k=0;k(buffer); + } + ISSM_MPI_Barrier(IssmComm::GetComm()); + _printf0_("Done reading files, now computing min and max.\n"); + + /*We have agregated minx and max across the cluster, now gather across the cluster onto + *cpu0 and then compute statistics:*/ + for (int f=0;f(size); + IssmDouble* allmin=xNew(size); + + ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); + ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); + + /*Store broadcasted value for later computation of histograms:*/ + maxxs[counter]=allmax; + minxs[counter]=allmin; + } + } /*}}}*/ + } + + /*Now do the same for the time mean fields:*/ + for (int f=0;f(size); + IssmDouble* allmin=xNew(size); + + ISSM_MPI_Allreduce(maxx,allmax,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); + ISSM_MPI_Allreduce(minx,allmin,size,ISSM_MPI_PDOUBLE,ISSM_MPI_MIN,IssmComm::GetComm()); + + + /*Store for later use in histogram computation:*/ + maxmeans[f]=allmax; + minmeans[f]=allmin; + + } /*}}}*/ + } + + /*Now that we have the min and max, we can start binning. First allocate + * histograms, then start filling them:*/ + IssmDouble** histogram=xNew(nfields*nsteps); + IssmDouble** timehistogram=xNew(nfields); + + _printf0_("Start reading files again, this time binning values in the histogram:\n"); + /*Start opening files:*/ + for (int i=(lower_row+1);i<=upper_row;i++){ + _printf0_("reading file #: " << i << "\n"); + char file[1000]; + long int length; + char* buffer=NULL; + + /*string:*/ + sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); + + /*open file: */ + _printf0_(" opening file:\n"); + FILE* fid=fopen(file,"rb"); + if(fid==NULL)_error_("cound not open file: " << file << "\n"); + + /*figure out size of file, and read the whole thing:*/ + _printf0_(" reading file:\n"); + fseek (fid, 0, SEEK_END); + length = ftell (fid); + fseek (fid, 0, SEEK_SET); + buffer = xNew(length); + fread (buffer, sizeof(char), length, fid); + + /*close file:*/ + fclose (fid); + + /*create a memory stream with this buffer:*/ + _printf0_(" processing file:\n"); + fid=fmemopen(buffer, length, "rb"); + + /*start reading data from the buffer directly:*/ + for (int f=0;f(nbins); + IssmDouble ma=*maxxs[counter]; + IssmDouble mi=*minxs[counter]; + int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index--; + if(ma==mi)index=0; + //_printf_( index << "|" << scalar << "|" << mi << "|" << ma << "|" << nbins << "\n"); + localhistogram[index]++; + histogram[counter]=localhistogram; + } + else if (xtype[counter]==3){ + IssmDouble* localhistogram=xNewZeroInit(doublematsize*nbins); + IssmDouble* ma=maxxs[counter]; + IssmDouble* mi=minxs[counter]; + for (int k=0;k=mi[k]); _assert_(index(doublemat); + } + else _error_("cannot carry out statistics on type " << xtype[counter]); + } + else{ + if(xtype[counter]==1){ + IssmDouble* localhistogram=histogram[counter]; + IssmDouble ma=*maxxs[counter]; + IssmDouble mi=*minxs[counter]; + int index=(scalar-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; + if(ma==mi)index=0; + localhistogram[index]++; + } + else if (xtype[counter]==3){ + IssmDouble* localhistogram=histogram[counter]; + IssmDouble* ma=maxxs[counter]; + IssmDouble* mi=minxs[counter]; + for (int k=0;k(doublemat); + } + else _error_("cannot carry out statistics on type " << xtype[counter]); + } + } + } + _printf0_(" average in time:\n"); + + /*Deal with average in time: */ + for (int f=0;f(nbins); + IssmDouble ma=*maxmeans[f]; + IssmDouble mi=*minmeans[f]; + int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; + if(ma==mi)index=0; + localhistogram[index]++; + timehistogram[f]=localhistogram; + } + else{ + IssmDouble* localhistogram=timehistogram[f]; + IssmDouble ma=*maxmeans[f]; + IssmDouble mi=*minmeans[f]; + int index=(timemean-mi)/(ma-mi)*nbins; if (index==nbins)index=nbins-1; + if(ma==mi)index=0; + localhistogram[index]++; + } + } + else{ + fseek(fid,0,SEEK_SET); + IssmDouble* timemean=xNewZeroInit(doublematsize); + for (int j=0;j(doublemat); + } + + if(i==(lower_row+1)){ + IssmDouble* localhistogram=xNewZeroInit(doublematsize*nbins); + IssmDouble* ma=maxmeans[f]; + IssmDouble* mi=minmeans[f]; + + for (int k=0;k(buffer); + } + _printf0_("Start aggregating histogram:\n"); + + /*We have agregated histograms across the cluster, now gather them across the cluster onto + *cpu0: */ + for (int f=0;f(nbins); + + ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + + /*add to results:*/ + if(my_rank==0){ + char fieldname[1000]; + + sprintf(fieldname,"%s%s",fields[f],"Histogram"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,1,nbins,j+1,0)); + + sprintf(fieldname,"%s%s",fields[f],"Max"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*maxxs[counter],j+1,0)); + sprintf(fieldname,"%s%s",fields[f],"Min"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*minxs[counter],j+1,0)); + } + } + } /*}}}*/ + else{ /*deal with arrays:{{{*/ + + int size=xsize[counter0]; + for (int j=0;j(size*nbins); + + ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + xDelete(histo); + + /*add to results:*/ + if(my_rank==0){ + char fieldname[1000]; + + sprintf(fieldname,"%s%s",fields[f],"Histogram"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,size,nbins,j+1,0)); + + sprintf(fieldname,"%s%s",fields[f],"Max"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,maxxs[counter],size,1,j+1,0)); + sprintf(fieldname,"%s%s",fields[f],"Min"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,minxs[counter],size,1,j+1,0)); + } + } + } /*}}}*/ + } + _printf0_("Start aggregating time mean histogram:\n"); + + /*Now do the same for the time mean fields:*/ + for (int f=0;f(nbins); + + ISSM_MPI_Allreduce(histo,allhisto,nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_MAX,IssmComm::GetComm()); + + /*add to results at time step 1:*/ + if(my_rank==0){ + char fieldname[1000]; + + sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,1,nbins,1,0)); + + sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*maxmeans[f],1,0)); + sprintf(fieldname,"%s%s",fields[f],"TimeMeaMin"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,*minmeans[f],1,0)); + } + } /*}}}*/ + else{ /*deal with arrays:{{{*/ + + int size=meanxsize[f]; + + /*we are broadcasting double arrays:*/ + IssmDouble* histo=timehistogram[f]; + IssmDouble* allhisto=xNewZeroInit(size*nbins); + + ISSM_MPI_Allreduce(histo,allhisto,size*nbins,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,IssmComm::GetComm()); + xDelete(histo); + /*add to results at step 1:*/ + if(my_rank==0){ + char fieldname[1000]; + + sprintf(fieldname,"%s%s",fields[f],"TimeMeanHistogram"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allhisto,size,nbins,1,0)); + sprintf(fieldname,"%s%s",fields[f],"TimeMeanMax"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,maxmeans[f],size,1,1,0)); + sprintf(fieldname,"%s%s",fields[f],"TimeMeanMin"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,minmeans[f],size,1,1,0)); + } + } /*}}}*/ + } + _printf0_("Done aggregating time mean histogram:\n"); + IssmComm::SetComm(ISSM_MPI_COMM_WORLD); +} +/*}}}*/ +int ComputeMeanVariance(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/ + + int nsamples; + char* directory=NULL; + char* model=NULL; + char** fields=NULL; + int* steps=NULL; + int nsteps; + int nfields; + int range,lower_row,upper_row; + int nfilesperdirectory; + + /*intermediary:*/ + IssmDouble* doublemat=NULL; + int doublematsize; + IssmDouble scalar; + + /*computation of average and variance itself:*/ + IssmDouble* x = NULL; + IssmDouble* x2 = NULL; + IssmDouble** xs = NULL; + IssmDouble** xs2 = NULL; + int* xtype=NULL; + int* xsize=NULL; + + IssmDouble** meanx=NULL; + IssmDouble** meanx2=NULL; + int* meantype=NULL; + int* meansize=NULL; + + /*only work on the statistical communicator: */ + if (color==MPI_UNDEFINED)return 0; + + /*Retrieve parameters:*/ + parameters->FindParam(&nfilesperdirectory,QmuNfilesPerDirectoryEnum); + parameters->FindParam(&nsamples,QmuNsampleEnum); + parameters->FindParam(&directory,DirectoryNameEnum); + parameters->FindParam(&model,InputFileNameEnum); + parameters->FindParam(&fields,&nfields,FieldsEnum); + parameters->FindParam(&steps,&nsteps,StepsEnum); + + /*Get rank from the stat comm communicator:*/ + IssmComm::SetComm(statcomm); + int my_rank=IssmComm::GetRank(); + + /*Open files and read them complelety, in a distributed way:*/ + range=DetermineLocalSize(nsamples,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm()); + + /*Initialize arrays:*/ + xs=xNew(nfields*nsteps); + xs2=xNew(nfields*nsteps); + xtype=xNew(nfields*nsteps); + xsize=xNew(nfields*nsteps); + + meantype=xNew(nfields); + meansize=xNew(nfields); + meanx=xNew(nfields); + meanx2=xNew(nfields); + + /*Start opening files:*/ + for (int i=(lower_row+1);i<=upper_row;i++){ + _printf0_("reading file #: " << i << "\n"); + char file[1000]; + long int length; + char* buffer=NULL; + + /*string:*/ + sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); + + /*open file: */ + _printf0_(" opening file: " << file << "\n"); + FILE* fid=fopen(file,"rb"); + if(fid==NULL) _error_(" could not open file: " << file << "\n"); + + /*figure out size of file, and read the whole thing:*/ + _printf0_(" reading file:\n"); + fseek (fid, 0, SEEK_END); + length = ftell (fid); + fseek (fid, 0, SEEK_SET); + buffer = xNew(length); + fread (buffer, sizeof(char), length, fid); + + /*close file:*/ + fclose (fid); + + /*create a memory stream with this buffer:*/ + _printf0_(" processing file:\n"); + fid=fmemopen(buffer, length, "rb"); + + /*start reading data from the buffer directly:*/ + for (int f=0;f(1); + xs2[counter]=xNew(1); + *xs[counter]=scalar; + *xs2[counter]=pow(scalar,2.0); + xsize[counter]=1; + } + else if (xtype[counter]==3){ + IssmDouble* doublemat2=xNew(doublematsize); + for(int k=0;k(1); + meanx2[f]=xNewZeroInit(1); + meansize[f]=1; + } + else{ + meanx[f]=xNewZeroInit(doublematsize); + meanx2[f]=xNewZeroInit(doublematsize); + meansize[f]=doublematsize; + } + } + fseek(fid,0,SEEK_SET); + if(meantype[f]==1){ + IssmDouble sc=0; + IssmDouble sc2=0; + for(int j=0;j(doublematsize); + IssmDouble* timemean2=xNewZeroInit(doublematsize); + + for(int j=0;j(buffer); + } + ISSM_MPI_Barrier(IssmComm::GetComm()); + _printf0_("Done reading files, now computing mean and variance.\n"); + + /*We have agregated x and x^2 across the cluster, now gather across the cluster onto + *cpu0 and then compute statistics:*/ + for (int f=0;fAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,j+1,0)); + sprintf(fieldname,"%s%s",fields[f],"Stddev"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,j+1,0)); + } + + } + } /*}}}*/ + else{ /*deal with arrays:{{{*/ + + int size=xsize[counter0]; + + IssmDouble* mean=xNew(size); + IssmDouble* stddev=xNew(size); + + for (int j=0;j(size); + IssmDouble* sumx2=xNew(size); + + ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); + ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); + + /*Build average and standard deviation. For standard deviation, use the + *following formula: sigma^2=E(x^2)-mu^2:*/ + for (int k=0;kAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,size,1,j+1,0)); + sprintf(fieldname,"%s%s",fields[f],"Stddev"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,size,1,j+1,0)); + } + } + } /*}}}*/ + } + /*Do the same but for the time mean:*/ + for (int f=0;fAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,1,0)); + sprintf(fieldname,"%s%s",fields[f],"TimeStddev"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,1,0)); + } + } /*}}}*/ + else{ /*deal with arrays:{{{*/ + + int size=meansize[f]; + IssmDouble* mean=xNew(size); + IssmDouble* stddev=xNew(size); + + /*we are broadcasting double arrays:*/ + x=meanx[f]; + x2=meanx2[f]; + + IssmDouble* sumx=xNew(size); + IssmDouble* sumx2=xNew(size); + + ISSM_MPI_Reduce(x,sumx,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); + ISSM_MPI_Reduce(x2,sumx2,size,ISSM_MPI_PDOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm()); + + /*Build average and standard deviation. For standard deviation, use the + *following formula: sigma^2=E(x^2)-mu^2:*/ + for (int k=0;kAddResult(new GenericExternalResult(results->Size()+1,fieldname,mean,size,1,1,0)); + sprintf(fieldname,"%s%s",fields[f],"TimeStddev"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,stddev,size,1,1,0)); + } + } /*}}}*/ + } + + _printf0_("Done with MeanVariance :\n"); + IssmComm::SetComm(ISSM_MPI_COMM_WORLD); + +} /*}}}*/ +int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm){ /*{{{*/ + + int nsamples; + char* directory=NULL; + char* model=NULL; + char** fields=NULL; + int* steps=NULL; + int nsteps; + int nfields; + int range,lower_row,upper_row; + int nfilesperdirectory; + int* indices=NULL; + int nindices; + + /*intermediary:*/ + IssmDouble* doublemat=NULL; + int doublematsize; + IssmDouble scalar; + + /*computation of average and variance itself:*/ + IssmDouble* x = NULL; + IssmDouble* allx=NULL; + IssmDouble** xs = NULL; + int* xtype=NULL; + int* xsize=NULL; + + /*only work on the statistical communicator: */ + if (color==MPI_UNDEFINED)return 0; + + /*Retrieve parameters:*/ + parameters->FindParam(&nsamples,QmuNsampleEnum); + parameters->FindParam(&directory,DirectoryNameEnum); + parameters->FindParam(&model,InputFileNameEnum); + parameters->FindParam(&fields,&nfields,FieldsEnum); + parameters->FindParam(&steps,&nsteps,StepsEnum); + parameters->FindParam(&indices,&nindices,IndicesEnum); + + /*Get rank from the stat comm communicator:*/ + IssmComm::SetComm(statcomm); + int my_rank=IssmComm::GetRank(); + + /*Open files and read them complelety, in a distributed way:*/ + range=DetermineLocalSize(nsamples,IssmComm::GetComm()); + GetOwnershipBoundariesFromRange(&lower_row,&upper_row,range,IssmComm::GetComm()); + + /*Initialize arrays:*/ + xs=xNew(nfields*nsteps); + xtype=xNew(nfields*nsteps); + xsize=xNew(nfields*nsteps); + + /*Start opening files:*/ + for (int i=(lower_row+1);i<=upper_row;i++){ + _printf0_("reading file #: " << i << "\n"); + char file[1000]; + long int length; + char* buffer=NULL; + + /*string:*/ + sprintf(file,"%s/%i/%s.outbin.%i",directory,my_rank+1,model,i); + + /*open file: */ + _printf0_(" opening file:\n"); + FILE* fid=fopen(file,"rb"); + + /*figure out size of file, and read the whole thing:*/ + _printf0_(" reading file:\n"); + fseek (fid, 0, SEEK_END); + length = ftell (fid); + fseek (fid, 0, SEEK_SET); + buffer = xNew(length); + fread (buffer, sizeof(char), length, fid); + + /*close file:*/ + fclose (fid); + + /*create a memory stream with this buffer:*/ + _printf0_(" processing file:\n"); + fid=fmemopen(buffer, length, "rb"); + + /*start reading data from the buffer directly:*/ + for (int f=0;f(range); + x[0]=scalar; + xs[counter]=x; + xsize[counter]=range; + } + else if (xtype[counter]==3){ + x=xNew(nindices*range); + for(int k=0;k(buffer); + } + ISSM_MPI_Barrier(IssmComm::GetComm()); + _printf0_("Done reading files, now assembling time series.\n"); + + for (int f=0;f(nsamples); + MPI_Gather(x, range, ISSM_MPI_PDOUBLE,allx, range, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm()); + /*add to results:*/ + if(my_rank==0){ + char fieldname[1000]; + + sprintf(fieldname,"%s%s",fields[f],"Samples"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allx,nsamples,1,j+1,0)); + } + } + else{ + /*we are broadcasting double arrays:*/ + x=xs[counter]; + allx=xNew(nsamples*nindices); + + MPI_Gather(x, range*nindices, ISSM_MPI_PDOUBLE,allx, range*nindices, ISSM_MPI_PDOUBLE, 0, IssmComm::GetComm()); + + /*add to results:*/ + if(my_rank==0){ + char fieldname[1000]; + sprintf(fieldname,"%s%s",fields[f],"Samples"); + results->AddResult(new GenericExternalResult(results->Size()+1,fieldname,allx,nsamples,nindices,j+1,0)); + } + } + } + } + _printf0_("Done with SampleSeries :\n"); + IssmComm::SetComm(ISSM_MPI_COMM_WORLD); + +} /*}}}*/ +int OutputStatistics(Parameters* parameters,Results* results,int color,ISSM_MPI_Comm statcomm){ /*{{{*/ + + char outputfilename[1000]; + char* directory=NULL; + char* model=NULL; + char* method=NULL; + int nsamples; + int* steps=NULL; + int nsteps; + + /*only work on the statistical communicator: */ + if (color==MPI_UNDEFINED)return 0; + + FemModel* femmodel=new FemModel(); + + /*Some parameters that will allow us to use the OutputResultsx module:*/ + parameters->AddObject(new BoolParam(QmuIsdakotaEnum,false)); + parameters->AddObject(new BoolParam(SettingsIoGatherEnum,true)); + + parameters->FindParam(&directory,DirectoryNameEnum); + parameters->FindParam(&model,InputFileNameEnum); + parameters->FindParam(&nsamples,QmuNsampleEnum); + parameters->FindParam(&steps,&nsteps,StepsEnum); + + sprintf(outputfilename,"%s/%s.stats",directory,model); + parameters->AddObject(new StringParam(OutputFileNameEnum,outputfilename)); + + /*Call OutputResults module:*/ + femmodel->parameters=parameters; + femmodel->results=results; + + OutputResultsx(femmodel); +} /*}}}*/ +bool DakotaDirStructure(int argc,char** argv){ /*{{{*/ + + char* input_file; + FILE* fid; + IoModel* iomodel=NULL; + int check; + + //qmu statistics + bool statistics = false; + int numdirectories = 0; + + /*First things first, set the communicator as a global variable: */ + IssmComm::SetComm(MPI_COMM_WORLD); + + /*Barrier:*/ + ISSM_MPI_Barrier(IssmComm::GetComm()); + _printf0_("Preparing directory structure for model outputs:" << "\n"); + + //open model input file for reading + input_file=xNew((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2)); + sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin"); + fid=fopen(input_file,"rb"); + if (fid==NULL) Cerr << "dirstructure error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl; + + //initialize IoModel, but light version, we just need it to fetch one constant: + iomodel=new IoModel(); + iomodel->fid=fid; + iomodel->FetchConstants(); + + //early return if statistics not requested: + iomodel->FindConstant(&statistics,"md.qmu.statistics"); + if(!statistics){ + delete iomodel; + fclose(fid); + return false; //important return value! + } + + iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories"); + + /*Ok, we have everything we need to create the directory structure:*/ + if(IssmComm::GetRank()==0){ + for (int i=0;i((strlen(argv[2])+strlen(argv[3])+strlen(".bin")+2)); + sprintf(input_file,"%s/%s%s",argv[2],argv[3],".bin"); + fid=fopen(input_file,"rb"); + if (fid==NULL) Cerr << "issm_dakota_statistics error message: could not open model " << input_file << " to retrieve qmu statistics parameters" << std::endl; + + //initialize IoModel, but light version, we'll need it to fetch constants: + iomodel=new IoModel(); + iomodel->fid=fid; + iomodel->FetchConstants(); + + //early return if statistics not requested: + iomodel->FindConstant(&statistics,"md.qmu.statistics"); + if(!statistics){ + delete iomodel; + fclose(fid); + return 0; + } + + //create parameters datasets with al the qmu statistics settings we need: + if(statistics){ + + /*Initialize parameters and results:*/ + results = new Results(); + parameters=new Parameters(); + + //solution type: + parameters->AddObject(new IntParam(SolutionTypeEnum,StatisticsSolutionEnum)); + + //root directory + directory=xNew(strlen(argv[2])+1); + xMemCpy(directory,argv[2],strlen(argv[2])+1); + parameters->AddObject(new StringParam(DirectoryNameEnum,directory)); + + //model name + model=xNew(strlen(argv[3])+1); + xMemCpy(model,argv[3],strlen(argv[3])+1); + parameters->AddObject(new StringParam(InputFileNameEnum,model)); + + //nsamples + iomodel->FindConstant(&nsamples,"md.qmu.method.params.samples"); + parameters->AddObject(new IntParam(QmuNsampleEnum,nsamples)); + + //ndirectories + iomodel->FindConstant(&numdirectories,"md.qmu.statistics.ndirectories"); + parameters->AddObject(new IntParam(QmuNdirectoriesEnum,numdirectories)); + + //nfiles per directory + iomodel->FindConstant(&nfilesperdirectory,"md.qmu.statistics.nfiles_per_directory"); + parameters->AddObject(new IntParam(QmuNfilesPerDirectoryEnum,nfilesperdirectory)); + + //At this point, we don't want to go forward any longer, we want to create an MPI + //communicator on which to carry out the computations: + if ((my_rank+1)*nfilesperdirectory>nsamples)color=MPI_UNDEFINED; + else color=0; + ISSM_MPI_Comm_split(ISSM_MPI_COMM_WORLD,color, my_rank, &statcomm); + + + iomodel->FindConstant(&numstatistics,"md.qmu.statistics.numstatistics"); + for (int i=1;i<=numstatistics;i++){ + + char* directory=NULL; + char* model=NULL; + int nsamples; + _printf0_("Dealing with qmu statistical computation #" << i << "\n"); + + sprintf(string,"md.qmu.statistics.method(%i).name",i); + iomodel->FindConstant(&name,string); + + sprintf(string,"md.qmu.statistics.method(%i).fields",i); + iomodel->FindConstant(&fields,&nfields,string); + parameters->AddObject(new StringArrayParam(FieldsEnum,fields,nfields)); + + sprintf(string,"md.qmu.statistics.method(%i).steps",i); + iomodel->FetchData(&steps,&dummy,&nsteps,string); + parameters->AddObject(new IntVecParam(StepsEnum,steps,nsteps)); + + if (strcmp(name,"Histogram")==0){ + /*fetch nbins: */ + sprintf(string,"md.qmu.statistics.method(%i).nbins",i); + iomodel->FindConstant(&nbins,string); + parameters->AddObject(new IntParam(NbinsEnum,nbins)); + ComputeHistogram(parameters,results,color,statcomm); + } + else if (strcmp(name,"SampleSeries")==0){ + /*fetch indices: */ + sprintf(string,"md.qmu.statistics.method(%i).indices",i); + iomodel->FetchData(&indices,&dummy,&nindices,string); + parameters->AddObject(new IntVecParam(IndicesEnum,indices,nindices)); + + ComputeSampleSeries(parameters,results,color,statcomm); + } + else if (strcmp(name,"MeanVariance")==0){ + ComputeMeanVariance(parameters,results,color,statcomm); + } + else _error_(" error creating qmu statistics methods parameters: unsupported method " << name); + } + } + //close model file: + fclose(fid); + + /*output results:*/ + OutputStatistics(parameters,results,color,statcomm); + + /*all meet here: */ + ISSM_MPI_Barrier(ISSM_MPI_COMM_WORLD); _printf0_("Output file.\n"); + + /*Delete ressources:*/ + delete parameters; + delete results; + + + +} /*}}}*/ Index: /issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.h =================================================================== --- /issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.h (revision 25627) +++ /issm/trunk-jpl/src/c/modules/QmuStatisticsx/QmuStatisticsx.h (revision 25627) @@ -0,0 +1,19 @@ +/*!\file: QmuStatisticsx.h + */ + +#ifndef _QMU_STATISTCSX_H_ +#define _QMU_STATISTCSX_H_ + +#include "../../classes/classes.h" + +int ComputeMeanVariance(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm); +int ComputeSampleSeries(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm); +int OutputStatistics(Parameters* parameters,Results* results,int color,ISSM_MPI_Comm statcomm); +int ComputeHistogram(Parameters* parameters,Results* results,int color, ISSM_MPI_Comm statcomm); +int readdata(IssmDouble** pdoublemat, int* pdoublematsize, IssmDouble* pdouble, FILE* fid,char* field,int step); +bool DakotaDirStructure(int argc,char** argv); +int DakotaStatistics(int argc,char** argv); + +/* local prototypes: */ + +#endif /* _QMU_STATISTCSX_H_ */ Index: /issm/trunk-jpl/src/c/modules/modules.h =================================================================== --- /issm/trunk-jpl/src/c/modules/modules.h (revision 25626) +++ /issm/trunk-jpl/src/c/modules/modules.h (revision 25627) @@ -81,4 +81,5 @@ #include "./PointCloudFindNeighborsx/PointCloudFindNeighborsx.h" #include "./PropagateFlagsFromConnectivityx/PropagateFlagsFromConnectivityx.h" +#include "./QmuStatisticsx/QmuStatisticsx.h" #include "./Reduceloadx/Reduceloadx.h" #include "./Reducevectorgtofx/Reducevectorgtofx.h" Index: /issm/trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25626) +++ /issm/trunk-jpl/src/c/shared/Enum/Enum.vim (revision 25627) @@ -112,4 +112,9 @@ syn keyword cConstant ControlInputSizeNEnum syn keyword cConstant ControlInputInterpolationEnum +syn keyword cConstant CumBslrEnum +syn keyword cConstant CumBslrIceEnum +syn keyword cConstant CumBslrHydroEnum +syn keyword cConstant CumGmtslrEnum +syn keyword cConstant CumGmslrEnum syn keyword cConstant DamageC1Enum syn keyword cConstant DamageC2Enum @@ -142,4 +147,5 @@ syn keyword cConstant ExtrapolationVariableEnum syn keyword cConstant FemModelCommEnum +syn keyword cConstant FieldsEnum syn keyword cConstant FlowequationFeFSEnum syn keyword cConstant FlowequationIsFSEnum @@ -211,4 +217,6 @@ syn keyword cConstant IndexEnum syn keyword cConstant InputFileNameEnum +syn keyword cConstant DirectoryNameEnum +syn keyword cConstant IndicesEnum syn keyword cConstant InputToDepthaverageInEnum syn keyword cConstant InputToDepthaverageOutEnum @@ -283,4 +291,5 @@ syn keyword cConstant MeshNumberofverticesEnum syn keyword cConstant ModelIdEnum +syn keyword cConstant NbinsEnum syn keyword cConstant NodesEnum syn keyword cConstant NumModelsEnum @@ -300,4 +309,5 @@ syn keyword cConstant QmuOutputEnum syn keyword cConstant QmuCurrEvalIdEnum +syn keyword cConstant QmuNsampleEnum syn keyword cConstant QmuResponsedescriptorsEnum syn keyword cConstant QmuVariableDescriptorsEnum @@ -307,11 +317,17 @@ syn keyword cConstant QmuResponsePartitionsEnum syn keyword cConstant QmuResponsePartitionsNpartEnum +syn keyword cConstant QmuStatisticsEnum +syn keyword cConstant QmuNumstatisticsEnum +syn keyword cConstant QmuNdirectoriesEnum +syn keyword cConstant QmuNfilesPerDirectoryEnum +syn keyword cConstant QmuStatisticsMethodEnum +syn keyword cConstant QmuMethodsEnum syn keyword cConstant RestartFileNameEnum syn keyword cConstant ResultsEnum syn keyword cConstant RootPathEnum +syn keyword cConstant ModelnameEnum syn keyword cConstant SaveResultsEnum syn keyword cConstant SolidearthPlanetRadiusEnum syn keyword cConstant SolidearthPlanetAreaEnum -syn keyword cConstant SealevelEustaticEnum syn keyword cConstant SolidearthSettingsAbstolEnum syn keyword cConstant RotationalAngularVelocityEnum @@ -412,4 +428,5 @@ syn keyword cConstant SteadystateRequestedOutputsEnum syn keyword cConstant StepEnum +syn keyword cConstant StepsEnum syn keyword cConstant StressbalanceAbstolEnum syn keyword cConstant StressbalanceFSreconditioningEnum @@ -688,4 +705,9 @@ syn keyword cConstant SealevelNEsaRateEnum syn keyword cConstant SealevelRSLEnum +syn keyword cConstant BslrEnum +syn keyword cConstant BslrIceEnum +syn keyword cConstant BslrHydroEnum +syn keyword cConstant BslrRateEnum +syn keyword cConstant GmtslrEnum syn keyword cConstant SealevelRSLEustaticEnum syn keyword cConstant SealevelRSLEustaticRateEnum @@ -1282,4 +1304,5 @@ syn keyword cConstant SpcTransientEnum syn keyword cConstant SsetEnum +syn keyword cConstant StatisticsSolutionEnum syn keyword cConstant SteadystateSolutionEnum syn keyword cConstant StressIntensityFactorEnum @@ -1360,5 +1383,4 @@ syn keyword cType Cfsurfacesquare syn keyword cType Channel -syn keyword cType classes syn keyword cType Constraint syn keyword cType Constraints @@ -1367,6 +1389,6 @@ syn keyword cType ControlInput syn keyword cType Covertree +syn keyword cType DataSetParam syn keyword cType DatasetInput -syn keyword cType DataSetParam syn keyword cType Definition syn keyword cType DependentObject @@ -1381,6 +1403,6 @@ syn keyword cType ElementInput syn keyword cType ElementMatrix +syn keyword cType ElementVector syn keyword cType Elements -syn keyword cType ElementVector syn keyword cType ExponentialVariogram syn keyword cType ExternalResult @@ -1389,10 +1411,9 @@ syn keyword cType Friction syn keyword cType Gauss -syn keyword cType GaussianVariogram -syn keyword cType gaussobjects syn keyword cType GaussPenta syn keyword cType GaussSeg syn keyword cType GaussTetra syn keyword cType GaussTria +syn keyword cType GaussianVariogram syn keyword cType GenericExternalResult syn keyword cType GenericOption @@ -1409,5 +1430,4 @@ syn keyword cType IssmDirectApplicInterface syn keyword cType IssmParallelDirectApplicInterface -syn keyword cType krigingobjects syn keyword cType Load syn keyword cType Loads @@ -1420,5 +1440,4 @@ syn keyword cType Matice syn keyword cType Matlitho -syn keyword cType matrixobjects syn keyword cType MatrixParam syn keyword cType Misfit @@ -1433,6 +1452,6 @@ syn keyword cType Observations syn keyword cType Option +syn keyword cType OptionUtilities syn keyword cType Options -syn keyword cType OptionUtilities syn keyword cType Param syn keyword cType Parameters @@ -1448,11 +1467,11 @@ syn keyword cType Regionaloutput syn keyword cType Results +syn keyword cType RiftStruct syn keyword cType Riftfront -syn keyword cType RiftStruct syn keyword cType SealevelMasks syn keyword cType Seg syn keyword cType SegInput +syn keyword cType SegRef syn keyword cType Segment -syn keyword cType SegRef syn keyword cType SpcDynamic syn keyword cType SpcStatic @@ -1473,4 +1492,8 @@ syn keyword cType Vertex syn keyword cType Vertices +syn keyword cType classes +syn keyword cType gaussobjects +syn keyword cType krigingobjects +syn keyword cType matrixobjects syn keyword cType AdjointBalancethickness2Analysis syn keyword cType AdjointBalancethicknessAnalysis @@ -1491,6 +1514,6 @@ syn keyword cType FreeSurfaceBaseAnalysis syn keyword cType FreeSurfaceTopAnalysis +syn keyword cType GLheightadvectionAnalysis syn keyword cType GiaAnalysis -syn keyword cType GLheightadvectionAnalysis syn keyword cType HydrologyDCEfficientAnalysis syn keyword cType HydrologyDCInefficientAnalysis Index: /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25626) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 25627) @@ -106,4 +106,9 @@ ControlInputSizeNEnum, ControlInputInterpolationEnum, + CumBslrEnum, + CumBslrIceEnum, + CumBslrHydroEnum, + CumGmtslrEnum, + CumGmslrEnum, DamageC1Enum, DamageC2Enum, @@ -136,4 +141,5 @@ ExtrapolationVariableEnum, FemModelCommEnum, + FieldsEnum, FlowequationFeFSEnum, FlowequationIsFSEnum, @@ -205,4 +211,6 @@ IndexEnum, InputFileNameEnum, + DirectoryNameEnum, + IndicesEnum, InputToDepthaverageInEnum, InputToDepthaverageOutEnum, @@ -277,4 +285,5 @@ MeshNumberofverticesEnum, ModelIdEnum, + NbinsEnum, NodesEnum, NumModelsEnum, @@ -294,4 +303,5 @@ QmuOutputEnum, QmuCurrEvalIdEnum, + QmuNsampleEnum, QmuResponsedescriptorsEnum, QmuVariableDescriptorsEnum, @@ -301,11 +311,17 @@ QmuResponsePartitionsEnum, QmuResponsePartitionsNpartEnum, + QmuStatisticsEnum, + QmuNumstatisticsEnum, + QmuNdirectoriesEnum, + QmuNfilesPerDirectoryEnum, + QmuStatisticsMethodEnum, + QmuMethodsEnum, RestartFileNameEnum, ResultsEnum, RootPathEnum, + ModelnameEnum, SaveResultsEnum, SolidearthPlanetRadiusEnum, SolidearthPlanetAreaEnum, - SealevelEustaticEnum, SolidearthSettingsAbstolEnum, RotationalAngularVelocityEnum, @@ -406,4 +422,5 @@ SteadystateRequestedOutputsEnum, StepEnum, + StepsEnum, StressbalanceAbstolEnum, StressbalanceFSreconditioningEnum, @@ -684,4 +701,9 @@ SealevelNEsaRateEnum, SealevelRSLEnum, + BslrEnum, + BslrIceEnum, + BslrHydroEnum, + BslrRateEnum, + GmtslrEnum, SealevelRSLEustaticEnum, SealevelRSLEustaticRateEnum, @@ -1281,4 +1303,5 @@ SpcTransientEnum, SsetEnum, + StatisticsSolutionEnum, SteadystateSolutionEnum, StressIntensityFactorEnum, Index: /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 25627) @@ -114,4 +114,9 @@ case ControlInputSizeNEnum : return "ControlInputSizeN"; case ControlInputInterpolationEnum : return "ControlInputInterpolation"; + case CumBslrEnum : return "CumBslr"; + case CumBslrIceEnum : return "CumBslrIce"; + case CumBslrHydroEnum : return "CumBslrHydro"; + case CumGmtslrEnum : return "CumGmtslr"; + case CumGmslrEnum : return "CumGmslr"; case DamageC1Enum : return "DamageC1"; case DamageC2Enum : return "DamageC2"; @@ -144,4 +149,5 @@ case ExtrapolationVariableEnum : return "ExtrapolationVariable"; case FemModelCommEnum : return "FemModelComm"; + case FieldsEnum : return "Fields"; case FlowequationFeFSEnum : return "FlowequationFeFS"; case FlowequationIsFSEnum : return "FlowequationIsFS"; @@ -213,4 +219,6 @@ case IndexEnum : return "Index"; case InputFileNameEnum : return "InputFileName"; + case DirectoryNameEnum : return "DirectoryName"; + case IndicesEnum : return "Indices"; case InputToDepthaverageInEnum : return "InputToDepthaverageIn"; case InputToDepthaverageOutEnum : return "InputToDepthaverageOut"; @@ -285,4 +293,5 @@ case MeshNumberofverticesEnum : return "MeshNumberofvertices"; case ModelIdEnum : return "ModelId"; + case NbinsEnum : return "Nbins"; case NodesEnum : return "Nodes"; case NumModelsEnum : return "NumModels"; @@ -302,4 +311,5 @@ case QmuOutputEnum : return "QmuOutput"; case QmuCurrEvalIdEnum : return "QmuCurrEvalId"; + case QmuNsampleEnum : return "QmuNsample"; case QmuResponsedescriptorsEnum : return "QmuResponsedescriptors"; case QmuVariableDescriptorsEnum : return "QmuVariableDescriptors"; @@ -309,11 +319,17 @@ case QmuResponsePartitionsEnum : return "QmuResponsePartitions"; case QmuResponsePartitionsNpartEnum : return "QmuResponsePartitionsNpart"; + case QmuStatisticsEnum : return "QmuStatistics"; + case QmuNumstatisticsEnum : return "QmuNumstatistics"; + case QmuNdirectoriesEnum : return "QmuNdirectories"; + case QmuNfilesPerDirectoryEnum : return "QmuNfilesPerDirectory"; + case QmuStatisticsMethodEnum : return "QmuStatisticsMethod"; + case QmuMethodsEnum : return "QmuMethods"; case RestartFileNameEnum : return "RestartFileName"; case ResultsEnum : return "Results"; case RootPathEnum : return "RootPath"; + case ModelnameEnum : return "Modelname"; case SaveResultsEnum : return "SaveResults"; case SolidearthPlanetRadiusEnum : return "SolidearthPlanetRadius"; case SolidearthPlanetAreaEnum : return "SolidearthPlanetArea"; - case SealevelEustaticEnum : return "SealevelEustatic"; case SolidearthSettingsAbstolEnum : return "SolidearthSettingsAbstol"; case RotationalAngularVelocityEnum : return "RotationalAngularVelocity"; @@ -414,4 +430,5 @@ case SteadystateRequestedOutputsEnum : return "SteadystateRequestedOutputs"; case StepEnum : return "Step"; + case StepsEnum : return "Steps"; case StressbalanceAbstolEnum : return "StressbalanceAbstol"; case StressbalanceFSreconditioningEnum : return "StressbalanceFSreconditioning"; @@ -690,4 +707,9 @@ case SealevelNEsaRateEnum : return "SealevelNEsaRate"; case SealevelRSLEnum : return "SealevelRSL"; + case BslrEnum : return "Bslr"; + case BslrIceEnum : return "BslrIce"; + case BslrHydroEnum : return "BslrHydro"; + case BslrRateEnum : return "BslrRate"; + case GmtslrEnum : return "Gmtslr"; case SealevelRSLEustaticEnum : return "SealevelRSLEustatic"; case SealevelRSLEustaticRateEnum : return "SealevelRSLEustaticRate"; @@ -1284,4 +1306,5 @@ case SpcTransientEnum : return "SpcTransient"; case SsetEnum : return "Sset"; + case StatisticsSolutionEnum : return "StatisticsSolution"; case SteadystateSolutionEnum : return "SteadystateSolution"; case StressIntensityFactorEnum : return "StressIntensityFactor"; Index: /issm/trunk-jpl/src/c/shared/Enum/README =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/README (revision 25626) +++ /issm/trunk-jpl/src/c/shared/Enum/README (revision 25627) @@ -1,3 +1,2 @@ - DO NOT MODIFY the following files: - StringToEnumx.cpp @@ -6,4 +5,7 @@ TO ADD AN ENUM: -add the corresponding lines to EnumDefinitions.h -and execute synchronizing script. +- add the desired, new Enums to EnumDefinitions.h +- execute `./Synchronize.sh` + +NOTE: Even when Synchronize.sh is run from this directory, ISSM_DIR must be + defined correctly Index: /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25626) +++ /issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 25627) @@ -114,4 +114,9 @@ else if (strcmp(name,"ControlInputSizeN")==0) return ControlInputSizeNEnum; else if (strcmp(name,"ControlInputInterpolation")==0) return ControlInputInterpolationEnum; + else if (strcmp(name,"CumBslr")==0) return CumBslrEnum; + else if (strcmp(name,"CumBslrIce")==0) return CumBslrIceEnum; + else if (strcmp(name,"CumBslrHydro")==0) return CumBslrHydroEnum; + else if (strcmp(name,"CumGmtslr")==0) return CumGmtslrEnum; + else if (strcmp(name,"CumGmslr")==0) return CumGmslrEnum; else if (strcmp(name,"DamageC1")==0) return DamageC1Enum; else if (strcmp(name,"DamageC2")==0) return DamageC2Enum; @@ -132,13 +137,13 @@ else if (strcmp(name,"DomainDimension")==0) return DomainDimensionEnum; else if (strcmp(name,"DomainType")==0) return DomainTypeEnum; - else if (strcmp(name,"DslModel")==0) return DslModelEnum; + else stage=2; + } + if(stage==2){ + if (strcmp(name,"DslModel")==0) return DslModelEnum; else if (strcmp(name,"DslModelid")==0) return DslModelidEnum; else if (strcmp(name,"DslNummodels")==0) return DslNummodelsEnum; else if (strcmp(name,"DslComputeFingerprints")==0) return DslComputeFingerprintsEnum; else if (strcmp(name,"EarthId")==0) return EarthIdEnum; - else stage=2; - } - if(stage==2){ - if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum; + else if (strcmp(name,"EplZigZagCounter")==0) return EplZigZagCounterEnum; else if (strcmp(name,"EsaHElastic")==0) return EsaHElasticEnum; else if (strcmp(name,"EsaHemisphere")==0) return EsaHemisphereEnum; @@ -147,4 +152,5 @@ else if (strcmp(name,"ExtrapolationVariable")==0) return ExtrapolationVariableEnum; else if (strcmp(name,"FemModelComm")==0) return FemModelCommEnum; + else if (strcmp(name,"Fields")==0) return FieldsEnum; else if (strcmp(name,"FlowequationFeFS")==0) return FlowequationFeFSEnum; else if (strcmp(name,"FlowequationIsFS")==0) return FlowequationIsFSEnum; @@ -216,4 +222,6 @@ else if (strcmp(name,"Index")==0) return IndexEnum; else if (strcmp(name,"InputFileName")==0) return InputFileNameEnum; + else if (strcmp(name,"DirectoryName")==0) return DirectoryNameEnum; + else if (strcmp(name,"Indices")==0) return IndicesEnum; else if (strcmp(name,"InputToDepthaverageIn")==0) return InputToDepthaverageInEnum; else if (strcmp(name,"InputToDepthaverageOut")==0) return InputToDepthaverageOutEnum; @@ -252,5 +260,8 @@ else if (strcmp(name,"LoveNfreq")==0) return LoveNfreqEnum; else if (strcmp(name,"LoveR0")==0) return LoveR0Enum; - else if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; + else stage=3; + } + if(stage==3){ + if (strcmp(name,"LoveShNmax")==0) return LoveShNmaxEnum; else if (strcmp(name,"LoveShNmin")==0) return LoveShNminEnum; else if (strcmp(name,"MassFluxSegments")==0) return MassFluxSegmentsEnum; @@ -260,8 +271,5 @@ else if (strcmp(name,"MasstransportMinThickness")==0) return MasstransportMinThicknessEnum; else if (strcmp(name,"MasstransportNumRequestedOutputs")==0) return MasstransportNumRequestedOutputsEnum; - else stage=3; - } - if(stage==3){ - if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; + else if (strcmp(name,"MasstransportPenaltyFactor")==0) return MasstransportPenaltyFactorEnum; else if (strcmp(name,"MasstransportRequestedOutputs")==0) return MasstransportRequestedOutputsEnum; else if (strcmp(name,"MasstransportStabilization")==0) return MasstransportStabilizationEnum; @@ -291,4 +299,5 @@ else if (strcmp(name,"MeshNumberofvertices")==0) return MeshNumberofverticesEnum; else if (strcmp(name,"ModelId")==0) return ModelIdEnum; + else if (strcmp(name,"Nbins")==0) return NbinsEnum; else if (strcmp(name,"Nodes")==0) return NodesEnum; else if (strcmp(name,"NumModels")==0) return NumModelsEnum; @@ -308,4 +317,5 @@ else if (strcmp(name,"QmuOutput")==0) return QmuOutputEnum; else if (strcmp(name,"QmuCurrEvalId")==0) return QmuCurrEvalIdEnum; + else if (strcmp(name,"QmuNsample")==0) return QmuNsampleEnum; else if (strcmp(name,"QmuResponsedescriptors")==0) return QmuResponsedescriptorsEnum; else if (strcmp(name,"QmuVariableDescriptors")==0) return QmuVariableDescriptorsEnum; @@ -315,11 +325,17 @@ else if (strcmp(name,"QmuResponsePartitions")==0) return QmuResponsePartitionsEnum; else if (strcmp(name,"QmuResponsePartitionsNpart")==0) return QmuResponsePartitionsNpartEnum; + else if (strcmp(name,"QmuStatistics")==0) return QmuStatisticsEnum; + else if (strcmp(name,"QmuNumstatistics")==0) return QmuNumstatisticsEnum; + else if (strcmp(name,"QmuNdirectories")==0) return QmuNdirectoriesEnum; + else if (strcmp(name,"QmuNfilesPerDirectory")==0) return QmuNfilesPerDirectoryEnum; + else if (strcmp(name,"QmuStatisticsMethod")==0) return QmuStatisticsMethodEnum; + else if (strcmp(name,"QmuMethods")==0) return QmuMethodsEnum; else if (strcmp(name,"RestartFileName")==0) return RestartFileNameEnum; else if (strcmp(name,"Results")==0) return ResultsEnum; else if (strcmp(name,"RootPath")==0) return RootPathEnum; + else if (strcmp(name,"Modelname")==0) return ModelnameEnum; else if (strcmp(name,"SaveResults")==0) return SaveResultsEnum; else if (strcmp(name,"SolidearthPlanetRadius")==0) return SolidearthPlanetRadiusEnum; else if (strcmp(name,"SolidearthPlanetArea")==0) return SolidearthPlanetAreaEnum; - else if (strcmp(name,"SealevelEustatic")==0) return SealevelEustaticEnum; else if (strcmp(name,"SolidearthSettingsAbstol")==0) return SolidearthSettingsAbstolEnum; else if (strcmp(name,"RotationalAngularVelocity")==0) return RotationalAngularVelocityEnum; @@ -367,5 +383,8 @@ else if (strcmp(name,"SmbAveraging")==0) return SmbAveragingEnum; else if (strcmp(name,"SmbDesfac")==0) return SmbDesfacEnum; - else if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; + else stage=4; + } + if(stage==4){ + if (strcmp(name,"SmbDpermil")==0) return SmbDpermilEnum; else if (strcmp(name,"SmbDsnowIdx")==0) return SmbDsnowIdxEnum; else if (strcmp(name,"SmbCldFrac")==0) return SmbCldFracEnum; @@ -383,8 +402,5 @@ else if (strcmp(name,"SmbIsd18opd")==0) return SmbIsd18opdEnum; else if (strcmp(name,"SmbIsdelta18o")==0) return SmbIsdelta18oEnum; - else stage=4; - } - if(stage==4){ - if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum; + else if (strcmp(name,"SmbIsdensification")==0) return SmbIsdensificationEnum; else if (strcmp(name,"SmbIsfirnwarming")==0) return SmbIsfirnwarmingEnum; else if (strcmp(name,"SmbIsgraingrowth")==0) return SmbIsgraingrowthEnum; @@ -423,4 +439,5 @@ else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum; else if (strcmp(name,"Step")==0) return StepEnum; + else if (strcmp(name,"Steps")==0) return StepsEnum; else if (strcmp(name,"StressbalanceAbstol")==0) return StressbalanceAbstolEnum; else if (strcmp(name,"StressbalanceFSreconditioning")==0) return StressbalanceFSreconditioningEnum; @@ -489,5 +506,8 @@ else if (strcmp(name,"Approximation")==0) return ApproximationEnum; else if (strcmp(name,"BalancethicknessMisfit")==0) return BalancethicknessMisfitEnum; - else if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum; + else stage=5; + } + if(stage==5){ + if (strcmp(name,"BalancethicknessOmega0")==0) return BalancethicknessOmega0Enum; else if (strcmp(name,"BalancethicknessOmega")==0) return BalancethicknessOmegaEnum; else if (strcmp(name,"BalancethicknessSpcthickness")==0) return BalancethicknessSpcthicknessEnum; @@ -506,8 +526,5 @@ else if (strcmp(name,"BasalforcingsPicoBasinId")==0) return BasalforcingsPicoBasinIdEnum; else if (strcmp(name,"BasalforcingsPicoBoxId")==0) return BasalforcingsPicoBoxIdEnum; - else stage=5; - } - if(stage==5){ - if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum; + else if (strcmp(name,"BasalforcingsPicoOverturningCoeff")==0) return BasalforcingsPicoOverturningCoeffEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanOverturning")==0) return BasalforcingsPicoSubShelfOceanOverturningEnum; else if (strcmp(name,"BasalforcingsPicoSubShelfOceanSalinity")==0) return BasalforcingsPicoSubShelfOceanSalinityEnum; @@ -612,5 +629,8 @@ else if (strcmp(name,"GiaLithosphereThickness")==0) return GiaLithosphereThicknessEnum; else if (strcmp(name,"GiaMantleViscosity")==0) return GiaMantleViscosityEnum; - else if (strcmp(name,"NGia")==0) return NGiaEnum; + else stage=6; + } + if(stage==6){ + if (strcmp(name,"NGia")==0) return NGiaEnum; else if (strcmp(name,"NGiaRate")==0) return NGiaRateEnum; else if (strcmp(name,"UGia")==0) return UGiaEnum; @@ -629,8 +649,5 @@ else if (strcmp(name,"HydrologydcEplThicknessTransient")==0) return HydrologydcEplThicknessTransientEnum; else if (strcmp(name,"HydrologydcMaskEplactiveElt")==0) return HydrologydcMaskEplactiveEltEnum; - else stage=6; - } - if(stage==6){ - if (strcmp(name,"HydrologydcMaskEplactiveNode")==0) return HydrologydcMaskEplactiveNodeEnum; + else if (strcmp(name,"HydrologydcMaskEplactiveNode")==0) return HydrologydcMaskEplactiveNodeEnum; else if (strcmp(name,"HydrologydcMaskThawedElt")==0) return HydrologydcMaskThawedEltEnum; else if (strcmp(name,"HydrologydcMaskThawedNode")==0) return HydrologydcMaskThawedNodeEnum; @@ -705,4 +722,9 @@ else if (strcmp(name,"SealevelNEsaRate")==0) return SealevelNEsaRateEnum; else if (strcmp(name,"SealevelRSL")==0) return SealevelRSLEnum; + else if (strcmp(name,"Bslr")==0) return BslrEnum; + else if (strcmp(name,"BslrIce")==0) return BslrIceEnum; + else if (strcmp(name,"BslrHydro")==0) return BslrHydroEnum; + else if (strcmp(name,"BslrRate")==0) return BslrRateEnum; + else if (strcmp(name,"Gmtslr")==0) return GmtslrEnum; else if (strcmp(name,"SealevelRSLEustatic")==0) return SealevelRSLEustaticEnum; else if (strcmp(name,"SealevelRSLEustaticRate")==0) return SealevelRSLEustaticRateEnum; @@ -730,5 +752,8 @@ else if (strcmp(name,"SigmaNN")==0) return SigmaNNEnum; else if (strcmp(name,"SigmaVM")==0) return SigmaVMEnum; - else if (strcmp(name,"SmbA")==0) return SmbAEnum; + else stage=7; + } + if(stage==7){ + if (strcmp(name,"SmbA")==0) return SmbAEnum; else if (strcmp(name,"SmbAValue")==0) return SmbAValueEnum; else if (strcmp(name,"SmbAccumulation")==0) return SmbAccumulationEnum; @@ -752,8 +777,5 @@ else if (strcmp(name,"SmbDlwrf")==0) return SmbDlwrfEnum; else if (strcmp(name,"SmbDswrf")==0) return SmbDswrfEnum; - else stage=7; - } - if(stage==7){ - if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum; + else if (strcmp(name,"SmbDzAdd")==0) return SmbDzAddEnum; else if (strcmp(name,"SmbDz")==0) return SmbDzEnum; else if (strcmp(name,"SmbDzMin")==0) return SmbDzMinEnum; @@ -853,5 +875,8 @@ else if (strcmp(name,"SurfaceOld")==0) return SurfaceOldEnum; else if (strcmp(name,"SurfaceLogVelMisfit")==0) return SurfaceLogVelMisfitEnum; - else if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; + else stage=8; + } + if(stage==8){ + if (strcmp(name,"SurfaceLogVxVyMisfit")==0) return SurfaceLogVxVyMisfitEnum; else if (strcmp(name,"SurfaceObservation")==0) return SurfaceObservationEnum; else if (strcmp(name,"SurfaceRelVelMisfit")==0) return SurfaceRelVelMisfitEnum; @@ -875,8 +900,5 @@ else if (strcmp(name,"Vx")==0) return VxEnum; else if (strcmp(name,"VxMesh")==0) return VxMeshEnum; - else stage=8; - } - if(stage==8){ - if (strcmp(name,"VxObs")==0) return VxObsEnum; + else if (strcmp(name,"VxObs")==0) return VxObsEnum; else if (strcmp(name,"VyAverage")==0) return VyAverageEnum; else if (strcmp(name,"Vy")==0) return VyEnum; @@ -976,5 +998,8 @@ else if (strcmp(name,"Outputdefinition82")==0) return Outputdefinition82Enum; else if (strcmp(name,"Outputdefinition83")==0) return Outputdefinition83Enum; - else if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; + else stage=9; + } + if(stage==9){ + if (strcmp(name,"Outputdefinition84")==0) return Outputdefinition84Enum; else if (strcmp(name,"Outputdefinition85")==0) return Outputdefinition85Enum; else if (strcmp(name,"Outputdefinition86")==0) return Outputdefinition86Enum; @@ -998,8 +1023,5 @@ else if (strcmp(name,"Absolute")==0) return AbsoluteEnum; else if (strcmp(name,"AdaptiveTimestepping")==0) return AdaptiveTimesteppingEnum; - else stage=9; - } - if(stage==9){ - if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum; + else if (strcmp(name,"AdjointBalancethickness2Analysis")==0) return AdjointBalancethickness2AnalysisEnum; else if (strcmp(name,"AdjointBalancethicknessAnalysis")==0) return AdjointBalancethicknessAnalysisEnum; else if (strcmp(name,"AdjointHorizAnalysis")==0) return AdjointHorizAnalysisEnum; @@ -1099,5 +1121,8 @@ else if (strcmp(name,"FreeSurfaceBaseAnalysis")==0) return FreeSurfaceBaseAnalysisEnum; else if (strcmp(name,"FreeSurfaceTopAnalysis")==0) return FreeSurfaceTopAnalysisEnum; - else if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; + else stage=10; + } + if(stage==10){ + if (strcmp(name,"FrontalForcingsDefault")==0) return FrontalForcingsDefaultEnum; else if (strcmp(name,"FrontalForcingsRignot")==0) return FrontalForcingsRignotEnum; else if (strcmp(name,"Fset")==0) return FsetEnum; @@ -1121,8 +1146,5 @@ else if (strcmp(name,"GroundingOnly")==0) return GroundingOnlyEnum; else if (strcmp(name,"GroundinglineMassFlux")==0) return GroundinglineMassFluxEnum; - else stage=10; - } - if(stage==10){ - if (strcmp(name,"Gset")==0) return GsetEnum; + else if (strcmp(name,"Gset")==0) return GsetEnum; else if (strcmp(name,"Gsl")==0) return GslEnum; else if (strcmp(name,"HOApproximation")==0) return HOApproximationEnum; @@ -1222,5 +1244,8 @@ else if (strcmp(name,"MpiDense")==0) return MpiDenseEnum; else if (strcmp(name,"Mpi")==0) return MpiEnum; - else if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; + else stage=11; + } + if(stage==11){ + if (strcmp(name,"MpiSparse")==0) return MpiSparseEnum; else if (strcmp(name,"Mumps")==0) return MumpsEnum; else if (strcmp(name,"NoFrictionOnPartiallyFloating")==0) return NoFrictionOnPartiallyFloatingEnum; @@ -1244,8 +1269,5 @@ else if (strcmp(name,"P0DG")==0) return P0DGEnum; else if (strcmp(name,"P1DG")==0) return P1DGEnum; - else stage=11; - } - if(stage==11){ - if (strcmp(name,"P1P1")==0) return P1P1Enum; + else if (strcmp(name,"P1P1")==0) return P1P1Enum; else if (strcmp(name,"P1P1GLS")==0) return P1P1GLSEnum; else if (strcmp(name,"P1bubble")==0) return P1bubbleEnum; @@ -1314,4 +1336,5 @@ else if (strcmp(name,"SpcTransient")==0) return SpcTransientEnum; else if (strcmp(name,"Sset")==0) return SsetEnum; + else if (strcmp(name,"StatisticsSolution")==0) return StatisticsSolutionEnum; else if (strcmp(name,"SteadystateSolution")==0) return SteadystateSolutionEnum; else if (strcmp(name,"StressIntensityFactor")==0) return StressIntensityFactorEnum; @@ -1344,5 +1367,8 @@ else if (strcmp(name,"TotalSmb")==0) return TotalSmbEnum; else if (strcmp(name,"TotalSmbScaled")==0) return TotalSmbScaledEnum; - else if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; + else stage=12; + } + if(stage==12){ + if (strcmp(name,"TransientArrayParam")==0) return TransientArrayParamEnum; else if (strcmp(name,"TransientInput")==0) return TransientInputEnum; else if (strcmp(name,"TransientParam")==0) return TransientParamEnum; @@ -1367,8 +1393,5 @@ else if (strcmp(name,"CalvingdevCoeff")==0) return CalvingdevCoeffEnum; else if (strcmp(name,"DeviatoricStress")==0) return DeviatoricStressEnum; - else stage=12; - } - if(stage==12){ - if (strcmp(name,"EtaAbsGradient")==0) return EtaAbsGradientEnum; + else if (strcmp(name,"EtaAbsGradient")==0) return EtaAbsGradientEnum; else if (strcmp(name,"MeshZ")==0) return MeshZEnum; else if (strcmp(name,"NearestInterp")==0) return NearestInterpEnum; Index: /issm/trunk-jpl/src/c/shared/Enum/Synchronize.sh =================================================================== --- /issm/trunk-jpl/src/c/shared/Enum/Synchronize.sh (revision 25626) +++ /issm/trunk-jpl/src/c/shared/Enum/Synchronize.sh (revision 25627) @@ -1,6 +1,10 @@ #!/bin/bash -#Synchronize EnumToStringx.cpp and StringToEnumx.cpp and matlab Enums - -#Get all lines of EnumDefinitions2.h which hold Enum | remove all commas > put everything in file temp +################################################################################ +# Synchronize EnumToStringx.cpp and StringToEnumx.cpp and MATLAB Enums +# +# Get all lines of EnumDefinitions.h which hold Enum | remove all commas > put everything in file temp +# +# NOTE: Even when run from this directory, ISSM_DIR must be defined correctly. +# cat EnumDefinitions.h | grep -e "[0-9]Enum," -e "[a-zA-Z]Enum," -e "MaximumNumberOfDefinitionsEnum" | grep -v include | sed -e "s/,/ /g" | awk '{print $1}' > temp Index: /issm/trunk-jpl/src/m/classes/clusters/local.m =================================================================== --- /issm/trunk-jpl/src/m/classes/clusters/local.m (revision 25627) +++ /issm/trunk-jpl/src/m/classes/clusters/local.m (revision 25627) @@ -0,0 +1,77 @@ +%LOCAL cluster class definition +% +% Usage: +% cluster=local('name','astrid','np',3); +% cluster=local('name',oshostname(),'np',3,'login','username'); + +classdef local + properties (SetAccess=public) + % {{{ + name = ''; + np = 1; + codepath = [IssmConfig('ISSM_PREFIX') '/bin']; + etcpath = [issmdir() '/etc']; + executionpath = [issmdir() '/execution']; + verbose = 1; + shell = '/bin/sh'; + %}}} + end + methods + function cluster=local(varargin) % {{{ + + %use provided options to change fields + options=pairoptions(varargin{:}); + + %get name + cluster.name=getfieldvalue(options,'name',oshostname()); + + %OK get other fields + cluster=AssignObjectFields(pairoptions(varargin{:}),cluster); + end + %}}} + function disp(cluster) % {{{ + % display the object + disp(sprintf('class ''%s'' object ''%s'' = ',class(cluster),inputname(1))); + disp(sprintf(' name: %s',cluster.name)); + disp(sprintf(' np: %i',cluster.np)); + disp(sprintf(' codepath: %s',cluster.codepath)); + disp(sprintf(' executionpath: %s',cluster.executionpath)); + disp(sprintf(' etcpath: %s',cluster.etcpath)); + disp(sprintf(' verbose: %s',cluster.verbose)); + disp(sprintf(' shell: %s',cluster.shell)); + end + %}}} + function md = checkconsistency(cluster,md,solution,analyses) % {{{ + if cluster.np<1 + md = checkmessage(md,['number of processors should be at least 1']); + end + if isnan(cluster.np), + md = checkmessage(md,'number of processors should not be NaN!'); + end + end + %}}} + function BuildQueueScript(cluster,dirname,modelname,solution,io_gather,isvalgrind,isgprof,isdakota,isoceancoupling) % {{{ + % Which executable are we calling? + executable='issm.exe'; % default + + if isdakota, + executable='issm_dakota.exe'; + end + + fid=fopen([modelname '.queue'],'w'); + fprintf(fid,'#!%s\n',cluster.shell); + fprintf(fid,'mpiexec -np %i %s/%s %s %s %s \n',cluster.np,cluster.codepath,executable,solution,'./',modelname); + fclose(fid); + end + %}}} + function UploadQueueJob(cluster,modelname,dirname,filelist)% {{{ + %do nothing really. + end %}}} + function LaunchQueueJob(cluster,modelname,dirname,filelist,restart,batch)% {{{ + system(['source ' modelname '.queue']); + + end %}}} + function Download(cluster,dirname,filelist)% {{{ + end %}}} + end +end Index: /issm/trunk-jpl/src/m/classes/qmu.m =================================================================== --- /issm/trunk-jpl/src/m/classes/qmu.m (revision 25626) +++ /issm/trunk-jpl/src/m/classes/qmu.m (revision 25627) @@ -13,4 +13,5 @@ method = struct(); params = struct(); + statistics = qmustatistics(); results = struct(); numberofresponses = 0; @@ -225,4 +226,5 @@ return; end + WriteData(fid,prefix,'data',self.method.params.samples,'name','md.qmu.method.params.samples','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','numberofresponses','format','Integer'); WriteData(fid,prefix,'object',self,'fieldname','variabledescriptors','format','StringArray'); @@ -240,4 +242,6 @@ end WriteData(fid,prefix,'data',flag,'name','md.qmu.mass_flux_segments_present','format','Boolean'); + self.statistics.marshall(prefix,md,fid); + end % }}} function savemodeljs(self,fid,modelname) % {{{ Index: /issm/trunk-jpl/src/m/classes/qmustatistics.m =================================================================== --- /issm/trunk-jpl/src/m/classes/qmustatistics.m (revision 25627) +++ /issm/trunk-jpl/src/m/classes/qmustatistics.m (revision 25627) @@ -0,0 +1,162 @@ +%QMUSTATISTICS class definition +% +% Usage: +% stats=qmustatistics(); + +classdef qmustatistics + properties (SetAccess=public) + + %generic information: + nfiles_per_directory = 5; %number of files per output directory. + ndirectories = 50; %number of output directories. %ndirectories should be < numcpus + + method = struct('name',{'None'},'fields',{[]},'steps',{[]},'nbins',{NaN},'indices',{[]}); + + %name: name of method, one of 'None','Histogram','SampleSeries' or 'MeanVariance' + %fields: fields for the statistics being requested, ex: 'Sealevel','BslrIce','BslrHydro' + %steps: time steps at which each field statistic is computed, ex: [1:2,5,20] or 1:100 + %nbins: number of bins for 'Histogram' statistics + %indices: vertex indices at which to retrieve samples + + end + methods + function self = qmustatistics(varargin) % {{{ + switch nargin + case 0 + self=setdefaultparameters(self); + case 1 + inputstruct=varargin{1}; + list1 = properties('qmustatistics'); + list2 = fieldnames(inputstruct); + for i=1:length(list1) + fieldname = list1{i}; + if ismember(fieldname,list2), + self.(fieldname) = inputstruct.(fieldname); + end + end + otherwise + error('constructor not supported'); + end + end % }}} + function self = setdefaultparameters(self) % {{{ + + self.method.name='None'; + self.nfiles_per_directory = 5; %number of files per output directory. + self.ndirectories = 50; %number of output directories. %ndirectories should be < numcpus + + end % }}} + function md = checkconsistency(self,md,solution,analyses) % {{{ + + %Early return: + if md.qmu.isdakota==0; return; end; + if strcmpi(self.method(1).name,'None'),return;end; + + %Checks: + md = checkfield(md,'fieldname','qmu.statistics.nfiles_per_directory','>=',1); + if self.ndirectories>md.cluster.np, + error('qmustatistics consistency check: number of cluster cpus should be > to number of output directories'); + end + if (self.ndirectories*self.nfiles_per_directory)~=md.qmu.method.params.samples, + error('qmustatistics consistency check: number of directories x number of files per directory should be == to number of samples requested!'); + end + for i=1:length(self.method), + m=self.method(i); + if strcmpi(m.name,'Histogram'), + md = checkfield(md,'fieldname',sprintf('qmu.statistics.method(%i).nbins',i),'>=',1,'<=',md.qmu.method.params.samples); + end + for f=1:length(m.fields), + if ~ischar(m.fields{f}), + error(sprintf('qmustatistics consistency check error: qmu.statistics.method(%i).fields{%i} is not a string!',i,f)); + end + end + for s=1:length(m.steps), + if m.steps(s)<=0, + error(sprintf('qmustatistics consistency check error: qmu.statistics.method(%i).steps(%i) should be >0!',i,s)); + end + if m.steps(s)> md.mesh.numberofvertices + error(sprintf('qmustatistics consistency check error: qmu.statistics.method(%i).steps(%i) should be < md.mesh.numberofvertices!',i,s)); + end + end + + end + + end % }}} + function list=defaultoutputs(self,md) % {{{ + end % }}} + function disp(self) % {{{ + + disp(sprintf('qmustatistics: post-dakota run processing of QMU statistics:')); + + if strcmpi(self.method(1).name,'None'), return; end; + + %generic information: + fielddisplay(self,'nfiles_per_directory','number of files per output directory'); + fielddisplay(self,'ndirectories','number of output directories. %ndirectories should be < numcpus'); + + for i=1:length(self.method), + disp([' method #' num2str(i)]); + disp(self.method(i)); + end + + end % }}} + function marshall(self,prefix,md,fid) % {{{ + + if strcmpi(self.method(1).name,'None'), + WriteData(fid,prefix,'name','md.qmu.statistics','data',0,'format','Boolean'); + statistics=0; + return; + else + WriteData(fid,prefix,'name','md.qmu.statistics','data',1,'format','Boolean'); + statistics=1; + end + + if statistics, + WriteData(fid,prefix,'name','md.qmu.statistics.nfiles_per_directory','data',self.nfiles_per_directory,'format','Integer'); + WriteData(fid,prefix,'name','md.qmu.statistics.ndirectories','data',self.ndirectories,'format','Integer'); + WriteData(fid,prefix,'name','md.qmu.statistics.numstatistics','data',length(self.method),'format','Integer'); + for i=1:length(self.method), + m=self.method(i); + WriteData(fid,prefix,'name',sprintf('md.qmu.statistics.method(%i).name',i),'data',m.name,'Format','String'); + WriteData(fid,prefix,'data',m.fields,'name',sprintf('md.qmu.statistics.method(%i).fields',i),'format','StringArray'); + WriteData(fid,prefix,'data',m.steps,'name',sprintf('md.qmu.statistics.method(%i).steps',i),'format','IntMat','mattype',3); + + if strcmpi(m.name,'Histogram'), + WriteData(fid,prefix,'name',sprintf('md.qmu.statistics.method(%i).nbins',i),'data',m.nbins,'format','Integer'); + elseif strcmpi(m.name,'MeanVariance'), + %do nothing + elseif strcmpi(m.name,'SampleSeries'), + WriteData(fid,prefix,'data',m.indices,'name',sprintf('md.qmu.statistics.method(%i).indices',i),'format','IntMat','mattype',3); + else + error(sprintf('qmustatistics marshall error message: unknown type ''%s'' for qmu.statistics.method(%i)',m.name,i)); + end + + end + end + + end % }}} + function self = extrude(self,md) % {{{ + + end % }}} + function savemodeljs(self,fid,modelname) % {{{ + + writejs1Darray(fid,[modelname '.stressbalance.spcvx'],self.spcvx); + writejs1Darray(fid,[modelname '.stressbalance.spcvy'],self.spcvy); + writejs1Darray(fid,[modelname '.stressbalance.spcvz'],self.spcvz); + writejsdouble(fid,[modelname '.stressbalance.restol'],self.restol); + writejsdouble(fid,[modelname '.stressbalance.reltol'],self.reltol); + writejsdouble(fid,[modelname '.stressbalance.abstol'],self.abstol); + writejsdouble(fid,[modelname '.stressbalance.isnewton'],self.isnewton); + writejsdouble(fid,[modelname '.stressbalance.FSreconditioning'],self.FSreconditioning); + writejsdouble(fid,[modelname '.stressbalance.maxiter'],self.maxiter); + writejsdouble(fid,[modelname '.stressbalance.shelf_dampening'],self.shelf_dampening); + writejs1Darray(fid,[modelname '.stressbalance.vertex_pairing'],self.vertex_pairing); + writejsdouble(fid,[modelname '.stressbalance.penalty_factor'],self.penalty_factor); + writejsdouble(fid,[modelname '.stressbalance.rift_penalty_lock'],self.rift_penalty_lock); + writejsdouble(fid,[modelname '.stressbalance.rift_penalty_threshold'],self.rift_penalty_threshold); + writejs2Darray(fid,[modelname '.stressbalance.referential'],self.referential); + writejs2Darray(fid,[modelname '.stressbalance.loadingforce'],self.loadingforce); + writejscellstring(fid,[modelname '.stressbalance.requested_outputs'],self.requested_outputs); + + end % }}} + end +end Index: /issm/trunk-jpl/src/m/contrib/larour/confidenceintervals.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/confidenceintervals.m (revision 25627) +++ /issm/trunk-jpl/src/m/contrib/larour/confidenceintervals.m (revision 25627) @@ -0,0 +1,36 @@ +function level=confidenceintervals(x,h,threshold) +%CONFIDENCEINTERVALS: build confidence interval from histogram. + + nods=size(h,1); nbins=size(h,2); + + %make sure h is one column smaller than x + if size(x,2) ~= (size(h,2)+1), + error('x should be one column larger than the histogram'); + end + + %First build cdf: + c=[zeros(nods,1) cumsum(h,2)]; + + %In percentage: + threshold=threshold/100; + + flags=zeros(nods,nbins+1); + pos=find(c=xmin,1)-offset); +id2x=min(numel(xdata),posx(end)+offset); + +ymin=min(Y(:)); ymax=max(Y(:)); +posy=find(ydata>=ymin); +id1y=max(1,find(ydata<=ymax,1)-offset); +id2y=min(numel(ydata),posy(end)+offset); + +disp(sprintf(' -- BedMachine %s : loading %s',continent,string)); +data = double(ncread(morlighemnc,string,[id1x id1y],[id2x-id1x+1 id2y-id1y+1],[1 1]))'; +xdata=xdata(id1x:id2x); +ydata=ydata(id1y:id2y); +data(find(data==-999999))=NaN; + +disp(sprintf(' -- BedMachine %s : interpolating %s',continent,string)); +if strcmp(string,'mask') | strcmp(string,'source'), + %Need nearest neighbor to avoid interpolation between 0 and 2 + if strcmpi(continent,'greenland'), ydata=flipud(ydata); data=flipud(data); end + output = InterpFromGridToMesh(xdata,ydata,data,double(X),double(Y),0); +else + ydata=flipud(ydata); data=flipud(data); + output = InterpFromGridToMesh(xdata,ydata,data,double(X),double(Y),0); +end + +output(find(output==-999999))=NaN; Index: /issm/trunk-jpl/src/m/contrib/larour/interpBedmap2.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/interpBedmap2.m (revision 25627) +++ /issm/trunk-jpl/src/m/contrib/larour/interpBedmap2.m (revision 25627) @@ -0,0 +1,47 @@ +function [dataout] = interpBedmap2(X,Y,string), +%INTERPBEDMAP2 - interpolate bedmap2 data +% +% Available data: +% 1. bed is bed height +% 2. surface is surface height +% 3. thickness is ice thickness +% 4. icemask_grounded_and_shelves is a mask file showing the grounding line and the extent of the floating ice shelves +% 5. rockmask is a mask file showing rock outcrops +% 6. lakemask_vostok is a mask file showing the extent of the lake cavity of Lake Vostok +% 7. bed_uncertainty is the bed uncertainty grid shown in figure 12 of the manuscript +% 8. thickness_uncertainty_5km is the thickness uncertainty grid shown in figure 11 of the manuscript +% 9. data_coverage is a binary grid showing the dis tribution of ice thickness data used in the grid of ice thickness +% +% Usage: +% [dataout] = interpBedmap2(X,Y,string) + +%reqad data +%path = '/u/astrid-r1b/ModelData/BedMap2/bedmap2_bin/'; +%path = '~/issm-jpl/proj-group/ModelData/BedMap2/bedmap2_bin/'; +path = '/Users/larour/ModelData/BedMap2/bedmap2_bin/'; +fid=fopen([path '/bedmap2_' string '.flt'],'r','l'); +data=fread(fid,[6667,6667],'float32'); +fclose(fid); + +% define grid +if strcmp(string,'thickness_uncertainty_5km'), + ncols =1361; + nrows =1361; + xll =-3401000; + yll =-3402000; + gridsize =5000; +else + ncols =6667; + nrows =6667; + xll =-3333000; + yll =-3333000; + gridsize =1000; +end +x_m=xll+(0:1:ncols-1)'*gridsize; +y_m=yll+(0:1:nrows-1)'*gridsize; + +%Change default to NaN +data(find(data==-9999))=NaN; + +%Interpolate +dataout = InterpFromGridToMesh(x_m,y_m,flipud(data'),double(X),double(Y),0); Index: /issm/trunk-jpl/src/m/contrib/larour/mme_autotime_correlation_matrix.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/mme_autotime_correlation_matrix.m (revision 25627) +++ /issm/trunk-jpl/src/m/contrib/larour/mme_autotime_correlation_matrix.m (revision 25627) @@ -0,0 +1,53 @@ +function matrix=mme_time_correlation_matrix(mme,varargin); + + if nargin==2, + type=varargin{1}; + elseif nargin==3, + mme2=varargin{1}; + type=varargin{2}; + else + error('mme_time_correlation_matrix usage error: 2 or 3 arguments only allowed!'); + end + + if nargin==2, + + %Out of a multi model ensemble (nsamples x nsteps) of runs, build + %a temporal correlation matrix (of size nsteps x nsteps) + + nsamples=size(mme,1); + nsteps=size(mme,2); + + %initialize with 1 in the diagonal: + matrix=eye(nsteps,nsteps); + + %go through time steps, and fill up the top part. + for i=1:nsteps, + for j=i+1:nsteps, + matrix(i,j)=corr(mme(:,i),mme(:,j),'Type',type); + matrix(j,i)=matrix(i,j); + end + end + else + + %Same kind of computations, except it's not autocorrelation: + nsamples=size(mme,1); nsamples2=size(mme2,1); + nsteps=size(mme,2); nsteps2=size(mme2,2); + + if nsteps2~=nsteps, + error('number of time steps from both sample matrices should be identical!'); + end + if nsamples2~=nsamples, + error('number of samples from both sample matrices should be identical!'); + end + + %initialize with 1 in the diagonal: + matrix=zeros(nsteps,nsteps); + + %go through time steps, and fill up the top part. + for i=1:nsteps, + for j=i:nsteps, + matrix(i,j)=corr(mme(:,i),mme2(:,j),'Type',type); + matrix(j,i)=matrix(i,j); + end + end + end Index: sm/trunk-jpl/src/m/contrib/larour/mme_time_correlation_matrix.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/mme_time_correlation_matrix.m (revision 25626) +++ (revision ) @@ -1,18 +1,0 @@ -function matrix=mme_time_correlation_matrix(mme,type); - - %Out of a multi model ensemble (nsamples x nsteps) of runs, build - %a temporal correlation matrix (of size nsteps x nsteps) - - nsamples=size(mme,1); - nsteps=size(mme,2); - - %initialize with 1 in the diagonal: - matrix=eye(nsteps,nsteps); - - %go through time steps, and fill up the top part. - for i=1:nsteps, - for j=i+1:nsteps, - matrix(i,j)=corr(mme(:,i),mme(:,j),'Type',type); - matrix(j,i)=matrix(i,j); - end - end Index: /issm/trunk-jpl/src/m/contrib/larour/oceanmip.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/larour/oceanmip.m (revision 25626) +++ /issm/trunk-jpl/src/m/contrib/larour/oceanmip.m (revision 25627) @@ -46,5 +46,6 @@ pos=find(diff(time)<0); if isempty(pos), - pos=(length(time)-12*100+1):length(time); + %pos=(length(time)-12*100+1):length(time); + pos=1:length(time); else pos=[(pos+1):length(time)]; @@ -75,4 +76,6 @@ zostogam(year)=mean(zostoga(i-11:i)); end + %control against 2006: + zostogam=zostogam-zostogam(1); self.zostoga{end+1}=zostogam; clear zostoga; @@ -109,6 +112,16 @@ disp(['interpolating model ' self.model{i} ' onto model mesh']); + %If we have 182 cells in long, trim by 1 on each side: + if size(self.pbo{i},1)==182, + self.pbo{i}=self.pbo{i}(2:181,:,:); + self.zos{i}=self.zos{i}(2:181,:,:); + self.long{i}=self.long{i}(2:181,:,:); + self.lat{i}=self.lat{i}(2:181,:,:); + end + + %triangulation: - long=double(self.long{i}); long=long(:); + long=double(self.long{i}); long=long(:); + pos=find(long<0); long(pos)=long(pos)+360; lat=double(self.lat{i}); lat=lat(:); [newl,uniqpos]=unique([lat,long],'rows','stable'); @@ -122,4 +135,7 @@ ind=find(areas<1e-8); index(ind,:)=[]; + %fix if we have orphans + [index long lat dummy newpos]=FixMesh(index,long,lat,1:length(long)); + time=self.time{i}; %retrieve fields: @@ -133,4 +149,5 @@ parfor j=1:length(time), + %for j=1:length(time), if mod(j,10)==0, s=sprintf(' progress: %.2g ',j/length(time)*100); @@ -138,9 +155,14 @@ end - pbo=omip_pbo(:,:,j); pbo=pbo(:); pbo=pbo(uniqpos); + pbo=omip_pbo(:,:,j); pbo=pbo(:); pbo=pbo(uniqpos); zos=omip_zos(:,:,j); zos=zos(:); zos=zos(uniqpos); - mesh_pbo(:,j)=InterpFromMeshToMesh2d(index,long,lat,pbo,meshlong,meshlat); - mesh_zos(:,j)=InterpFromMeshToMesh2d(index,long,lat,zos,meshlong,meshlat); + pboj= InterpFromMeshToMesh2d(index,long,lat,pbo(newpos),meshlong,meshlat); + pos=find(abs(pboj)>1e4); pboj(pos)=0; + mesh_pbo(:,j)=pboj; + + zosj=InterpFromMeshToMesh2d(index,long,lat,zos(newpos),meshlong,meshlat); + pos=find(abs(zosj)>1e5); zosj(pos)=0; + mesh_zos(:,j)=zosj; end self.mesh_pbo{end+1}=mesh_pbo; @@ -200,6 +222,5 @@ pbo=self.mesh_pbo{i}; pbo=pbo/1000; %in meters end - time=self.time{i}; - array=[pbo;time]; + array=pbo; break; end Index: /issm/trunk-jpl/src/m/plot/plot_coastlines.m =================================================================== --- /issm/trunk-jpl/src/m/plot/plot_coastlines.m (revision 25626) +++ /issm/trunk-jpl/src/m/plot/plot_coastlines.m (revision 25627) @@ -2,5 +2,6 @@ %define coastline: {{{ -coast=[-83.83 -180 ;... +coast=[-90 -180; ... + -83.83 -180 ;... -84.33 -178 ;... -84.5 -174 ;... @@ -614,4 +615,7 @@ -83.5 176 ;... -83.83 180 ;... +-85 180 ;... +-87 180 ;... +-90 180 ;... NaN NaN ;... -77.58 166.33 ;... @@ -9866,8 +9870,9 @@ 19.77 -156.05 ;... 19.95 -155.87 ;... -20.27 -155.88 ]; +20.27 -155.88 ;... +NaN NaN]; -coastlat=coast(:,1); -coastlon=coast(:,2); +coastlat=[coast(:,1); coast(:,1)]; +coastlon=[coast(:,2); coast(:,2)-360]; %}}} @@ -9902,4 +9907,9 @@ p=polyshape(coastlon,coastlat); plot(p,'FaceColor','w','FaceAlpha',1); + if strcmpi(getfieldvalue(options,'coordcent','atlantic'),'pacific'), + xlim(getfieldvalue(options,'xlim',[-360 0])); + else + xlim(getfieldvalue(options,'xlim',[-180 180])); + end %p=plot(coastlon,coastlat,'k-'); end Index: /issm/trunk-jpl/src/m/plot/processdata.m =================================================================== --- /issm/trunk-jpl/src/m/plot/processdata.m (revision 25626) +++ /issm/trunk-jpl/src/m/plot/processdata.m (revision 25627) @@ -42,4 +42,10 @@ end +%special case for coord latlong: +if strcmpi(getfieldvalue(options,'coord','xy'),'latlon') | strcmpi(getfieldvalue(options,'coord','xy'),'latlong'), + [data datatype] = processdatalatlong(md,data,options); + return; +end + %needed later on if isprop(md.mesh,'numberofvertices2d'), Index: /issm/trunk-jpl/src/m/plot/processdatalatlong.m =================================================================== --- /issm/trunk-jpl/src/m/plot/processdatalatlong.m (revision 25627) +++ /issm/trunk-jpl/src/m/plot/processdatalatlong.m (revision 25627) @@ -0,0 +1,29 @@ +function [data datatype] = processdatalatlong(md,data,options); + + % datatype = 1 -> elements + % datatype = 2 -> nodes + %what is the mesh we are using: + x0=md.mesh.long; + y0=md.mesh.lat; + + %add row at lat=90 and lat=-90 + add=[(-180:.1:-1e-5)';(1e-5:.1:180)']; + nadd=length(add); + xextra=[add;add]; + yextra=[90*ones(nadd,1); -90*ones(nadd,1)]; + x=[x0;xextra]; + y=[y0;yextra]; + elements=delaunay(x,y); + + %with this mesh, interpolate data: + if length(data)==length(md.mesh.long), + datatype=2; + + %interpolate data: + extradata=griddata(x0,y0,data,xextra,yextra,'nearest'); + + data=[data; extradata]; + elseif length(data)==length(md.mesh.elements), + error('processdatalatlong error message: coord ''latlong'' case not covered for element data '); + datatype=1; + end Index: /issm/trunk-jpl/src/m/plot/processmesh.m =================================================================== --- /issm/trunk-jpl/src/m/plot/processmesh.m (revision 25626) +++ /issm/trunk-jpl/src/m/plot/processmesh.m (revision 25627) @@ -24,9 +24,20 @@ if strcmp(domaintype(md.mesh),'3Dsurface'), [x y z elements is2d isplanet] = processmesh(md.mesh,options); - if strcmpi(getfieldvalue(options,'coord','xy'),'latlon'), - x=md.mesh.long; - y=md.mesh.lat; + if strcmpi(getfieldvalue(options,'coord','xy'),'latlon') | strcmpi(getfieldvalue(options,'coord','xy'),'latlong'), + x0=md.mesh.long; + y0=md.mesh.lat; + %add row at lat=90 and lat=-90 + add=[(-180:.1:-1e-5)';(1e-5:.1:180)']; + nadd=length(add); + xextra=[add;add]; + yextra=[90*ones(nadd,1); -90*ones(nadd,1)]; + x=[x0;xextra]; + y=[y0;yextra]; + + if strcmpi(getfieldvalue(options,'coordcent','atlantic'),'pacific'), + pos=find(x>0); x(pos)=-360+x(pos); + end elements=delaunay(x,y); - z=md.mesh.lat; z(:)=0; + z=x; z(:)=0; end return; @@ -42,5 +53,5 @@ else elements=md.mesh.elements; - if ~strcmpi(getfieldvalue(options,'coord','xy'),'latlon') , + if ~strcmpi(getfieldvalue(options,'coord','xy'),'latlon') & ~strcmpi(getfieldvalue(options,'coord','xy'),'latlong') , x=md.mesh.x; if isprop(md.mesh,'x2d'), x2d=md.mesh.x2d; end @@ -50,4 +61,7 @@ x=md.mesh.long; y=md.mesh.lat; + if strcmpi(getfieldvalue(options,'coordcent','atlantic'),'pacific'), + pos=find(x>0); x(pos)-360+x(pos); + end end end Index: /issm/trunk-jpl/src/m/plot/subplotmodel.m =================================================================== --- /issm/trunk-jpl/src/m/plot/subplotmodel.m (revision 25626) +++ /issm/trunk-jpl/src/m/plot/subplotmodel.m (revision 25627) @@ -25,4 +25,11 @@ ymin = 1-vmargin(2)-height; +box=getfieldvalue(options,'box','off'); +if strcmpi(box,'on'), + visible='on'; +else + visible='off'; +end + for i = 1:nlines xmin = hmargin(1); @@ -30,5 +37,6 @@ if(((i-1)*ncols+j)==num) ha = axes('Units','normalized', ... - 'Position',[xmin ymin width height],'XTickLabel','','YTickLabel','','Visible','off'); + 'Position',[xmin ymin width height],'XTickLabel','','YTickLabel','',... + 'Visible',visible,'box',box); return end Index: /issm/trunk-jpl/src/m/qmu/postqmu.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/postqmu.m (revision 25626) +++ /issm/trunk-jpl/src/m/qmu/postqmu.m (revision 25627) @@ -36,5 +36,5 @@ end -if md.qmu.output, +if md.qmu.output & strcmpi(md.qmu.statistics.method(1).name,'None'), if strcmpi(md.qmu.method.method,'nond_sampling'), dakotaresults.modelresults={}; @@ -46,4 +46,9 @@ end end +end +if ~strcmpi(md.qmu.statistics.method(1).name,'None'), + md.qmu.isdakota=0; + md=loadresultsfromdisk(md,[md.miscellaneous.name,'.stats']); + md.qmu.isdakota=1; end Index: /issm/trunk-jpl/src/m/solve/loadresultsfromcluster.m =================================================================== --- /issm/trunk-jpl/src/m/solve/loadresultsfromcluster.m (revision 25626) +++ /issm/trunk-jpl/src/m/solve/loadresultsfromcluster.m (revision 25627) @@ -30,5 +30,5 @@ end end - if md.qmu.output, + if md.qmu.output & strcmpi(md.qmu.statistics.method(1).name,'None'), if strcmpi(md.qmu.method.method,'nond_sampling'), for i=1:md.qmu.method.params.samples @@ -36,4 +36,7 @@ end end + end + if ~strcmpi(md.qmu.statistics.method(1).name,'None'), + filelist{end+1}=[md.miscellaneous.name '.stats']; end else Index: /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.m (revision 25626) +++ /issm/trunk-jpl/src/m/solve/loadresultsfromdisk.m (revision 25627) @@ -35,4 +35,12 @@ error(['No result found in binary file ' filename '. Check for solution crash.']); end + if isempty(structure(1).SolutionType), + if ~isempty(structure(end).SolutionType), + structure(1).SolutionType=structure(end).SolutionType; + else + warning(['Cannot find a solution type in the results! Ascribing one: ''NoneSolution''.']); + structure(1).SolutionType='NoneSolution'; + end + end md.results.(structure(1).SolutionType)=structure; @@ -40,5 +48,5 @@ md.private.solution=structure(1).SolutionType; - %read log files onto fields + %read log files onto fields if exist([md.miscellaneous.name '.errlog'],'file'), md.results.(structure(1).SolutionType)(1).errlog=char(textread([md.miscellaneous.name '.errlog'],'%s','delimiter','\n')); Index: /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m =================================================================== --- /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 25626) +++ /issm/trunk-jpl/src/m/solve/parseresultsfromdisk.m (revision 25627) @@ -5,6 +5,55 @@ else results=parseresultsfromdiskioserial(md,filename); -end - + %results=parseresultsfromdiskioserialsequential(md,filename); +end + +function results=parseresultsfromdiskioserialsequential(md,filename) % {{{ + +%Open file +fid=fopen(filename,'rb'); +if(fid==-1), + error(['loadresultsfromdisk error message: could not open ',filename,' for binary reading']); +end + +%first pass to figure out the steps we have: +steps=[]; +while 1, + result = ReadDataDimensions(fid); + if isempty(result), + break; + end + if result.step~=-9999, + steps=[steps result.step]; + end +end + +steps=unique(steps); + +%create structure: +results=struct('step',num2cell(steps)); + +%second pass to fill the steps we have: +fseek(fid,0,-1); %rewind +while 1, + result = ReadData(fid,md); + if isempty(result), + break; + end + if result.step==-9999, + result.step=1; + result.time=0; + end + %find where we are putting this step: + ind=find(steps==result.step); + if isempty(ind), + error('could not find a step in our pre-structure! Something is very off!'); + end + + %plug: + results(ind).time=result.time; + results(ind).(result.fieldname)=result.field; +end +fclose(fid); +% }}} function results=parseresultsfromdiskioserial(md,filename) % {{{ Index: /issm/trunk-jpl/test/NightlyRun/test2004.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 25626) +++ /issm/trunk-jpl/test/NightlyRun/test2004.m (revision 25627) @@ -394,23 +394,31 @@ %eustatic run: -md.solidearth.settings.rigid=0; md.solidearth.settings.elastic=0;md.solidearth.settings.rotation=0; +md.solidearth.settings.rigid=0; +md.solidearth.settings.elastic=0; +md.solidearth.settings.rotation=0; md.solidearth.requested_outputs= {'default',... - 'SurfaceloadIceThicknessChange','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... - 'SealevelNEsaRate', 'SealevelUEsaRate', 'NGiaRate', 'UGiaRate','SealevelEustaticMask','SealevelEustaticOceanMask'}; + 'SurfaceloadIceThicknessChange','Sealevel','SealevelRSLRate','SealevelriseCumDeltathickness',... + 'SealevelNEsaRate', 'SealevelUEsaRate', 'NGiaRate', 'UGiaRate','SealevelEustaticMask','SealevelEustaticOceanMask'}; md=solve(md,'Sealevelrise'); Seustatic=md.results.SealevelriseSolution.Sealevel; %eustatic + rigid run: -md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=0;md.solidearth.settings.rotation=0; +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=0; +md.solidearth.settings.rotation=0; md=solve(md,'Sealevelrise'); Srigid=md.results.SealevelriseSolution.Sealevel; %eustatic + rigid + elastic run: -md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1;md.solidearth.settings.rotation=0; +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=1; +md.solidearth.settings.rotation=0; md=solve(md,'Sealevelrise'); Selastic=md.results.SealevelriseSolution.Sealevel; %eustatic + rigid + elastic + rotation run: -md.solidearth.settings.rigid=1; md.solidearth.settings.elastic=1; md.solidearth.settings.rotation=1; +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=1; +md.solidearth.settings.rotation=1; md=solve(md,'Sealevelrise'); Srotation=md.results.SealevelriseSolution.Sealevel; Index: /issm/trunk-jpl/test/NightlyRun/test2005.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2005.m (revision 25627) +++ /issm/trunk-jpl/test/NightlyRun/test2005.m (revision 25627) @@ -0,0 +1,107 @@ +%Test Name: EarthSlr + +%mesh earth: +md=model; +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh + +%parameterize solidearth solution: +%solidearth loading: {{{ +md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); +md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[1;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +%antarctica +late=sum(md.mesh.lat(md.mesh.elements),2)/3; +longe=sum(md.mesh.long(md.mesh.elements),2)/3; +pos=find(late < -80); +md.solidearth.surfaceload.icethicknesschange(pos)=-100; +%greenland +pos=find(late>70 & late<80 & longe>-60 & longe<-30); +md.solidearth.surfaceload.icethicknesschange(pos)=-100; + +%elastic loading from love numbers: +md.solidearth.lovenumbers=lovenumbers('maxdeg',100); + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); +icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); +icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(md.solidearth.surfaceload.icethicknesschange); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; + +%make sure wherever there is an ice load, that the mask is set to ice: +pos=find(md.solidearth.surfaceload.icethicknesschange); +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +% }}} + +md.solidearth.settings.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + + +%Miscellaneous +md.miscellaneous.name='test2005'; + +%Solution parameters +md.solidearth.settings.reltol=NaN; +md.solidearth.settings.abstol=1e-3; +md.solidearth.settings.computesealevelchange=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.solidearth.settings.maxiter=10; + +%eustatic + rigid + elastic + rotation run: +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=1; +md.solidearth.settings.rotation=1; + +%transient settings: +md.timestepping.start_time=0; +md.timestepping.final_time=10; +md.timestepping.time_step=1; +md.transient.isslr=1; +md.transient.issmb=0; +md.transient.isgia=1; +md.transient.ismasstransport=0; +md.transient.isstressbalance=0; +md.transient.isthermal=0; +dh=md.solidearth.surfaceload.icethicknesschange; +deltathickness=zeros(md.mesh.numberofelements+1,10); +for i=1:10, + deltathickness(1:end-1,i)=dh*i; +end +deltathickness(end,:)=0:1:9; +md.solidearth.surfaceload.icethicknesschange=deltathickness; + +%hack: +md.geometry.surface=zeros(md.mesh.numberofvertices,1); +md.geometry.thickness=ones(md.mesh.numberofvertices,1); +md.geometry.base=-ones(md.mesh.numberofvertices,1); +md.geometry.bed=md.geometry.base; + +%run transient solution: +md=solve(md,'Transient'); + +S1=md.results.TransientSolution(1).Sealevel; +S5=md.results.TransientSolution(5).Sealevel; +S10=md.results.TransientSolution(10).Sealevel; +Seus10=md.results.TransientSolution(10).Bslr; + +%Fields and tolerances to track changes +field_names={'Sealevel1','Sealevel5','Sealevel10','Seustatic10'}; +field_tolerances={1e-13,1e-13,1e-13,1e-13,1e-13}; +field_values={S1,S5,S10,Seus10}; + Index: /issm/trunk-jpl/test/NightlyRun/test2006.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2006.m (revision 25627) +++ /issm/trunk-jpl/test/NightlyRun/test2006.m (revision 25627) @@ -0,0 +1,197 @@ +%Test Name: EarthSlr Dakota Sampling glaciers. + +%mesh earth: +md=model; +md.cluster=generic('name',oshostname(),'np',5); +md.mesh=gmshplanet('radius',6.371012*10^3,'resolution',700.); %700 km resolution mesh + +%parameterize solidearth solution: +%solidearth loading: {{{ +md.solidearth.surfaceload.icethicknesschange=zeros(md.mesh.numberofelements,1); +md.solidearth.sealevel=zeros(md.mesh.numberofvertices,1); +md.dsl.global_average_thermosteric_sea_level_change=[1;0]; +md.dsl.sea_surface_height_change_above_geoid=zeros(md.mesh.numberofvertices+1,1); +md.dsl.sea_water_pressure_change_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); +%antarctica +late=sum(md.mesh.lat(md.mesh.elements),2)/3; +longe=sum(md.mesh.long(md.mesh.elements),2)/3; +pos=find(late < -80); +md.solidearth.surfaceload.icethicknesschange(pos)=-100; +%greenland +pos=find(late>70 & late<80 & longe>-60 & longe<-30); +md.solidearth.surfaceload.icethicknesschange(pos)=-100; + +%elastic loading from love numbers: +md.solidearth.lovenumbers=lovenumbers('maxdeg',100); + +%}}} +%mask: {{{ +mask=gmtmask(md.mesh.lat,md.mesh.long); +icemask=ones(md.mesh.numberofvertices,1); +pos=find(mask==0); +icemask(pos)=-1; +pos=find(sum(mask(md.mesh.elements),2)<3); +icemask(md.mesh.elements(pos,:))=-1; +md.mask.ice_levelset=icemask; +md.mask.ocean_levelset=-icemask; + +%make sure that the elements that have loads are fully grounded: +pos=find(md.solidearth.surfaceload.icethicknesschange); +md.mask.ocean_levelset(md.mesh.elements(pos,:))=1; + +%make sure wherever there is an ice load, that the mask is set to ice: +pos=find(md.solidearth.surfaceload.icethicknesschange); +md.mask.ice_levelset(md.mesh.elements(pos,:))=-1; +% }}} + +md.solidearth.settings.ocean_area_scaling=0; + +%Geometry for the bed, arbitrary: +md.geometry.bed=-ones(md.mesh.numberofvertices,1); + +%Materials: +md.materials=materials('hydro'); + + +%Miscellaneous +md.miscellaneous.name='test2006'; + +%Solution parameters +md.solidearth.settings.reltol=NaN; +md.solidearth.settings.abstol=1e-3; +md.solidearth.settings.computesealevelchange=1; + +% max number of iteration reverted back to 10 (i.e., the original default value) +md.solidearth.settings.maxiter=10; + +%eustatic + rigid + elastic + rotation run: +md.solidearth.settings.rigid=1; +md.solidearth.settings.elastic=1; +md.solidearth.settings.rotation=1; + +%transient settings: +md.timestepping.start_time=0; +md.timestepping.final_time=10; +md.timestepping.time_step=1; +md.transient.isslr=1; +md.transient.issmb=0; +md.transient.isgia=1; +md.transient.ismasstransport=0; +md.transient.isstressbalance=0; +md.transient.isthermal=0; +dh=md.solidearth.surfaceload.icethicknesschange; +deltathickness=zeros(md.mesh.numberofelements+1,10); +for i=1:10, + deltathickness(1:end-1,i)=dh*i; +end +deltathickness(end,:)=0:1:9; +md.solidearth.surfaceload.icethicknesschange=deltathickness; + +%hack: +md.geometry.surface=zeros(md.mesh.numberofvertices,1); +md.geometry.thickness=ones(md.mesh.numberofvertices,1); +md.geometry.base=-ones(md.mesh.numberofvertices,1); +md.geometry.bed=md.geometry.base; + + +%Uncertainty Quantification +%ice sheets {{{ +npart=1; nt=1; +partition=-ones(md.mesh.numberofelements,1); +pos=find(late < -80); partition(pos)=0; +pos=find(late>70 & late<80 & longe>-60 & longe<-30); partition(pos)=0; + +%variables: +qmuvar.surfaceload=normal_uncertain('descriptor','scaled_SurfaceloadIceThicknessChange',... + 'mean',1*ones(npart,nt),... + 'stddev',.1*ones(npart,nt),... %10% standard deviation + 'partition',partition,... + 'transient','on',... + 'nsteps',nt); +%}}} + +%correlation: +md.qmu.correlation_matrix=[]; + +%variables final declaration: +md.qmu.variables=struct(); +md.qmu.variables.surfaceload=qmuvar.surfaceload; + +locations=[1 5 10 15 20]; +%responses % {{{ +md.qmu.responses.sealevel1=response_function('descriptor','Outputdefinition1'); +md.qmu.responses.sealevel2=response_function('descriptor','Outputdefinition2'); +md.qmu.responses.sealevel3=response_function('descriptor','Outputdefinition3'); +md.qmu.responses.sealevel4=response_function('descriptor','Outputdefinition4'); +md.qmu.responses.sealevel5=response_function('descriptor','Outputdefinition5'); + +%output definitions: +for i=1:length(locations), + ind=locations(i); + if i==1, + md.outputdefinition.definitions={nodalvalue('name','SNode','definitionstring','Outputdefinition1', ... + 'model_string','Sealevel','node',ind)}; + else + md.outputdefinition.definitions{end+1}=nodalvalue('name','SNode','definitionstring',['Outputdefinition' num2str(i)], ... + 'model_string','Sealevel','node',ind); + end +end +% }}} + %algorithm: % {{{ + md.qmu.method =dakota_method('nond_samp'); + md.qmu.method(end)=dmeth_params_set(md.qmu.method(end),... + 'seed',1234,... + 'samples',10,... + 'sample_type','random'); + md.qmu.output=1; + %}}} + %parameters % {{{ + md.qmu.params.direct=true; + md.qmu.params.interval_type='forward'; + md.qmu.params.analysis_driver='matlab'; + md.qmu.params.evaluation_scheduling='master'; + md.qmu.params.processors_per_evaluation=2; + md.qmu.params.tabular_graphics_data=true; + md.qmu.isdakota=1; + md.verbose=verbose(0); md.verbose.qmu=1; + % }}} + %qmu statistics %{{{ + md.qmu.statistics.nfiles_per_directory=2; + md.qmu.statistics.ndirectories=5; + + md.qmu.statistics.method(1).name='Histogram'; + md.qmu.statistics.method(1).fields={'Sealevel','BslrIce'}; + md.qmu.statistics.method(1).steps=1:10; + md.qmu.statistics.method(1).nbins=20; + + md.qmu.statistics.method(2).name='MeanVariance'; + md.qmu.statistics.method(2).fields={'Sealevel','BslrIce'}; + md.qmu.statistics.method(2).steps=[1:10]; + + md.qmu.statistics.method(3).name='SampleSeries'; + md.qmu.statistics.method(3).fields={'Sealevel','BslrIce'}; + md.qmu.statistics.method(3).steps=[1:10]; + md.qmu.statistics.method(3).indices=locations; + %}}} + +%run transient dakota solution: +mds=solve(md,'Transient'); + +%run without statistics computations: +md.qmu.statistics.method(1).name='None'; +md=solve(md,'Transient'); + +%compare statistics with our own here: +svalues=mds.results.StatisticsSolution(end).SealevelSamples; %all values at locations. + +dvalues=zeros(md.qmu.method.params.samples,length(locations)); +for i=1:md.qmu.method.params.samples, + dvalues(i,:)=md.results.dakota.modelresults{i}.TransientSolution(end).Sealevel(locations); +end + +samplesnorm=norm(dvalues-svalues,'fro'); + +%Fields and tolerances to track changes +field_names={'Samples Norm'}; +field_tolerances={1e-13}; +field_values={samplesnorm}; Index: /issm/trunk-jpl/test/NightlyRun/test2010.m =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test2010.m (revision 25626) +++ /issm/trunk-jpl/test/NightlyRun/test2010.m (revision 25627) @@ -75,5 +75,5 @@ % uncomment following 2 lines for md=solve(md,'Sealevelrise'); -eus=md.results.SealevelriseSolution.SealevelRSLEustatic; +eus=md.results.SealevelriseSolution.Bslr; slr=md.results.SealevelriseSolution.Sealevel; moixz=md.results.SealevelriseSolution.SealevelInertiaTensorXZ;