Index: ../trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h =================================================================== --- ../trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h (revision 26109) +++ ../trunk-jpl/src/c/toolkits/issm/IssmAbsVec.h (revision 26110) @@ -51,6 +51,7 @@ virtual void PointwiseDivide(IssmAbsVec* x,IssmAbsVec* y)=0; virtual void PointwiseMult(IssmAbsVec* x,IssmAbsVec* y)=0; virtual void Pow(doubletype scale_factor)=0; + virtual void Sum(doubletype* pvalue)=0; }; #endif //#ifndef _ISSM_ABS_VEC_H_ Index: ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h =================================================================== --- ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 26109) +++ ../trunk-jpl/src/c/toolkits/issm/IssmMpiVec.h (revision 26110) @@ -593,6 +593,10 @@ } /*}}}*/ + void Sum(doubletype* pvalue){/*{{{*/ + _error_("not support yet!"); + } + /*}}}*/ void BucketsBuildScatterBuffers(int** pnumvalues_forcpu,int** prow_indices_forcpu,doubletype** pvalues_forcpu,int** pmodes_forcpu,DataSet** bucketsforcpu,int num_procs){/*{{{*/ /*intermediary: */ Index: ../trunk-jpl/src/c/toolkits/issm/IssmVec.h =================================================================== --- ../trunk-jpl/src/c/toolkits/issm/IssmVec.h (revision 26109) +++ ../trunk-jpl/src/c/toolkits/issm/IssmVec.h (revision 26110) @@ -219,6 +219,10 @@ vector->Pow(scale_factor); } /*}}}*/ + void Sum(doubletype* pvalue){/*{{{*/ + vector->Sum(pvalue); + } + /*}}}*/ }; #endif //#ifndef _ISSMVEC_H_ Index: ../trunk-jpl/src/c/toolkits/objects/Vector.h =================================================================== --- ../trunk-jpl/src/c/toolkits/objects/Vector.h (revision 26109) +++ ../trunk-jpl/src/c/toolkits/objects/Vector.h (revision 26110) @@ -405,5 +405,16 @@ else this->ivector->Pow(scale_factor); } /*}}}*/ -}; +void Sum(doubletype* pvalue){ /*{{{*/ + _assert_(this);/*{{{*/ + + if(type==PetscVecType){ + #ifdef _HAVE_PETSC_ + this->pvector->Sum(pvalue); + #endif + } + else this->ivector->Sum(pvalue); +} +/*}}}*/ +}; /*}}}*/ #endif //#ifndef _VECTOR_H_ Index: ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h =================================================================== --- ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h (revision 26109) +++ ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.h (revision 26110) @@ -57,6 +57,7 @@ IssmDouble Max(void); void Scale(IssmDouble scale_factor); void Pow(IssmDouble scale_factor); + void Sum(IssmDouble* pvalue); void PointwiseDivide(PetscVec* x,PetscVec* y); void PointwiseMult(PetscVec* x,PetscVec* y); IssmDouble Dot(PetscVec* vector); Index: ../trunk-jpl/src/c/Makefile.am =================================================================== --- ../trunk-jpl/src/c/Makefile.am (revision 26109) +++ ../trunk-jpl/src/c/Makefile.am (revision 26110) @@ -96,6 +96,7 @@ ./classes/Vertex.cpp \ ./classes/Hook.cpp \ ./classes/Radar.cpp \ + ./classes/BarystaticContributions.cpp \ ./classes/ExternalResults/Results.cpp \ ./classes/Elements/Element.cpp \ ./classes/Elements/Elements.cpp \ Index: ../trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp =================================================================== --- ../trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26109) +++ ../trunk-jpl/src/c/analyses/SealevelchangeAnalysis.cpp (revision 26110) @@ -182,7 +182,9 @@ parameters->AddObject(new DoubleMatParam(CumBslcOceanPartitionEnum,bslcocean_partition,npartocean,1)); xDelete(partitionocean); } - + /*New optimized code:*/ + BarystaticContributions* barystaticcontributions=new BarystaticContributions(iomodel); + parameters->AddObject(new GenericParam(barystaticcontributions,BarystaticContributionsEnum)); /*Deal with external multi-model ensembles: {{{*/ if(isexternal){ Index: ../trunk-jpl/src/c/classes/BarystaticContributions.cpp =================================================================== --- ../trunk-jpl/src/c/classes/BarystaticContributions.cpp (nonexistent) +++ ../trunk-jpl/src/c/classes/BarystaticContributions.cpp (revision 26110) @@ -0,0 +1,149 @@ +/*!\file BarystaticContributions.c + * \brief: implementation of the BarystaticContributions object + */ + +/*Include files: {{{*/ +#ifdef HAVE_CONFIG_H + #include +#else +#error "Cannot compile with HAVE_CONFIG_H symbol! run configure first!" +#endif + +#include "./BarystaticContributions.h" +#include "../toolkits/toolkits.h" +#include "./classes.h" +/*}}}*/ + +/*Constructors and destructors:*/ +BarystaticContributions::BarystaticContributions(IoModel* iomodel ){ /*{{{*/ + + int nel; + + iomodel->FetchData(&nice,"md.solidearth.npartice"); + if(nice){ + iomodel->FetchData(&pice,&nel,NULL,"md.solidearth.partitionice"); + ice=new Vector(nice); + cumice=new Vector(nice); cumice->Set(0); cumice->Assemble(); + } + + iomodel->FetchData(&nhydro,"md.solidearth.nparthydro"); + if(nhydro){ + iomodel->FetchData(&phydro,&nel,NULL,"md.solidearth.partitionhydro"); + hydro=new Vector(nhydro); + cumhydro=new Vector(nhydro); cumhydro->Set(0); cumhydro->Assemble(); + } + iomodel->FetchData(&nocean,"md.solidearth.npartocean"); + if(nocean){ + iomodel->FetchData(&pocean,&nel,NULL,"md.solidearth.partitionocean"); + ocean=new Vector(nocean); + cumocean=new Vector(nocean); cumocean->Set(0); cumocean->Assemble(); + } + +} /*}}}*/ +BarystaticContributions::~BarystaticContributions(){ /*{{{*/ + delete ice; delete cumice; + delete hydro; delete cumhydro; + delete ocean; delete cumocean; + if(nice)xDelete(pice); + if(nhydro)xDelete(phydro); + if(nocean)xDelete(pocean); +}; /*}}}*/ + +/*Support routines:*/ +IssmDouble BarystaticContributions::Total(){ /*{{{*/ + + IssmDouble sumice,sumhydro,sumocean; + + ice->Assemble(); + hydro->Assemble(); + ocean->Assemble(); + + ice->Sum(&sumice); + hydro->Sum(&sumhydro); + ocean->Sum(&sumocean); + + return sumice+sumhydro+sumocean; + +} /*}}}*/ +IssmDouble BarystaticContributions::CumTotal(){ /*{{{*/ + + IssmDouble sumice,sumhydro,sumocean; + + cumice->Assemble(); + cumhydro->Assemble(); + cumocean->Assemble(); + + cumice->Sum(&sumice); + cumhydro->Sum(&sumhydro); + cumocean->Sum(&sumocean); + + + return sumice+sumhydro+sumocean; + +} /*}}}*/ +void BarystaticContributions::Cumulate(Parameters* parameters){ /*{{{*/ + + cumice->AXPY(ice,1); + cumocean->AXPY(ocean,1); + cumhydro->AXPY(hydro,1); + + +} /*}}}*/ +void BarystaticContributions::Save(Results* results, Parameters* parameters, IssmDouble oceanarea){ /*{{{*/ + + int step; + IssmDouble time; + IssmDouble rho_water; + + IssmDouble* cumice_serial=NULL; + IssmDouble* cumhydro_serial=NULL; + IssmDouble* cumocean_serial=NULL; + + IssmDouble sumice,sumhydro,sumocean; + + parameters->FindParam(&step,StepEnum); + parameters->FindParam(&time,TimeEnum); + parameters->FindParam(&rho_water,TimeEnum); + + ice->Sum(&sumice); hydro->Sum(&sumhydro); ocean->Sum(&sumocean); + results->AddResult(new GenericExternalResult(results->Size()+1,BslcEnum,-this->Total()/oceanarea/rho_water,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,BslcIceEnum,-sumice/oceanarea/rho_water,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,BslcHydroEnum,-sumice/oceanarea/rho_water,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,BslcOceanEnum,-sumocean/oceanarea/rho_water,step,time)); + + cumice->Sum(&sumice); cumhydro->Sum(&sumhydro); cumocean->Sum(&sumocean); + results->AddResult(new GenericExternalResult(results->Size()+1,CumBslcEnum,this->CumTotal()/oceanarea/rho_water,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,CumBslcIceEnum,sumice/oceanarea/rho_water,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,CumBslcHydroEnum,sumhydro/oceanarea/rho_water,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,CumBslcOceanEnum,sumocean/oceanarea/rho_water,step,time)); + + cumice_serial=this->cumice->ToMPISerial0(); for (int i=0;icumhydro->ToMPISerial0(); for (int i=0;icumocean->ToMPISerial0(); for (int i=0;iAddResult(new GenericExternalResult(results->Size()+1,CumBslcIcePartitionEnum,cumice_serial,nice,1,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,CumBslcHydroPartitionEnum,cumhydro_serial,nhydro,1,step,time)); + results->AddResult(new GenericExternalResult(results->Size()+1,CumBslcOceanPartitionEnum,cumocean_serial,nocean,1,step,time)); + + if(IssmComm::GetRank()==0){ + xDelete(cumice_serial); + xDelete(cumhydro_serial); + xDelete(cumocean_serial); + } + return; + +} /*}}}*/ +void BarystaticContributions::Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue){ /*{{{*/ + + int id; + + id=reCast(pice[eid]); + ice->SetValue(id,icevalue,ADD_VAL); + + id=reCast(phydro[eid]); + hydro->SetValue(id,hydrovalue,ADD_VAL); + + id=reCast(pocean[eid]); + ocean->SetValue(id,oceanvalue,ADD_VAL); + +} /*}}}*/ Index: ../trunk-jpl/src/c/classes/BarystaticContributions.h =================================================================== --- ../trunk-jpl/src/c/classes/BarystaticContributions.h (nonexistent) +++ ../trunk-jpl/src/c/classes/BarystaticContributions.h (revision 26110) @@ -0,0 +1,46 @@ +/*!\file BarystaticContributions.h + * \brief: header file for barystatic contribution object + */ + +#ifndef _BARYSTATICCONTRIBUTIONS_H_ +#define _BARYSTATICCONTRIBUTIONS_H_ + +/*Headers:*/ +class IoModel; +class Parameters; +class Results; +template class Vector; +#include "../shared/shared.h" + +class BarystaticContributions { + + public: + + Vector* ice; //contributions to every ice partition + Vector* cumice; //cumulated contributions to every ice partition + int nice; //number of ice partitions + IssmDouble* pice; //ice partition + + Vector* hydro; //contributions to every hydro partition + Vector* cumhydro; //cumulated contributions to every hydro partition + int nhydro; //number of hydro partitions + IssmDouble* phydro; //hydro partition + + Vector* ocean; //contributions to every ocean partition + Vector* cumocean; //cumulated contributions to every ocean partition + int nocean; //number of ocean partitions + IssmDouble* pocean; //ocean partition + + /*BarystaticContributions constructors, destructors :*/ + BarystaticContributions(IoModel* iomodel ); + ~BarystaticContributions(); + + /*routines:*/ + IssmDouble Total(); + IssmDouble CumTotal(); + void Cumulate(Parameters* parameters); + void Save(Results* results, Parameters* parameters, IssmDouble oceanarea); + void Set(int eid, IssmDouble icevalue, IssmDouble hydrovalue, IssmDouble oceanvalue); + +}; +#endif /* _BARYSTATICCONTRIBUTIONS_H_ */ Index: ../trunk-jpl/src/c/classes/Elements/Element.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Element.h (revision 26109) +++ ../trunk-jpl/src/c/classes/Elements/Element.h (revision 26110) @@ -37,6 +37,7 @@ template class Vector; class ElementMatrix; class ElementVector; +class BarystaticContributions; /*}}}*/ class Element: public Object{ @@ -388,6 +389,12 @@ virtual void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old,SealevelMasks* mask)=0; virtual void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks)=0; virtual void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y)=0; + + virtual void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze)=0; + virtual void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks)=0; + virtual void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks)=0; + virtual void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks)=0; + virtual void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks)=0; #endif }; Index: ../trunk-jpl/src/c/classes/Elements/Penta.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 26109) +++ ../trunk-jpl/src/c/classes/Elements/Penta.h (revision 26110) @@ -226,6 +226,12 @@ void SealevelchangeSal(IssmDouble* Sgo, IssmDouble* Sg_old, SealevelMasks* masks){_error_("not implemented yet!");}; void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; + + void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");}; + void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; + void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");}; #endif /*}}}*/ Index: ../trunk-jpl/src/c/classes/Elements/Seg.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 26109) +++ ../trunk-jpl/src/c/classes/Elements/Seg.h (revision 26110) @@ -181,6 +181,13 @@ void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; + + void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");}; + void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; + void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");}; + #endif #ifdef _HAVE_DAKOTA_ Index: ../trunk-jpl/src/c/classes/Elements/Tetra.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26109) +++ ../trunk-jpl/src/c/classes/Elements/Tetra.h (revision 26110) @@ -188,6 +188,12 @@ void DeformationFromSurfaceLoads(IssmDouble* Up ,IssmDouble* North, IssmDouble* East, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet!");}; void GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x,IssmDouble* y){_error_("not implemented yet");}; + void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){_error_("not implemented yet");}; + void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; + void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){_error_("not implemented yet");}; + void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){_error_("not implemented yet");}; + #endif #ifdef _HAVE_DAKOTA_ Index: ../trunk-jpl/src/c/classes/Elements/Tria.cpp =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26109) +++ ../trunk-jpl/src/c/classes/Elements/Tria.cpp (revision 26110) @@ -5431,6 +5431,7 @@ /*}}}*/ #endif #ifdef _HAVE_SEALEVELCHANGE_ +//old code IssmDouble Tria::OceanAverage(IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/ if(masks->isoceanin[this->lid]){ @@ -5450,7 +5451,7 @@ } /*}}}*/ -void Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/ +void Tria::SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old, SealevelMasks* masks){/*{{{*/ /*early return if we are not on an ice cap OR ocean:*/ if(!masks->isiceonly[this->lid] && !masks->isoceanin[this->lid]){ dI_list[0] = 0.0; // this is important!!! @@ -5544,7 +5545,7 @@ return; }/*}}}*/ -void Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/ +void Tria::SetSealevelMasks(SealevelMasks* masks){ /*{{{*/ masks->isiceonly[this->lid]=this->IsIceOnlyInElement(); masks->isoceanin[this->lid]=this->IsOceanInElement(); @@ -5557,10 +5558,9 @@ /*are we not fully grounded: */ if ((gr_input->GetInputMin())<0) masks->notfullygrounded[this->lid]=true; else masks->notfullygrounded[this->lid]=false; - } /*}}}*/ -void Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/ +void Tria::SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/ /*diverse:*/ int gsize; @@ -5950,7 +5950,7 @@ return bslchydro; } /*}}}*/ -void Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/ +void Tria::SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks){ /*{{{*/ /*diverse:*/ int gsize; @@ -6000,7 +6000,7 @@ return; } /*}}}*/ -void Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/ +void Tria::SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old, SealevelMasks* masks){ /*{{{*/ /*diverse:*/ int gsize,dummy; @@ -6042,7 +6042,7 @@ return; } /*}}}*/ -void Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/ +void Tria::DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North ,IssmDouble* East,IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/ /*diverse:*/ int gsize; @@ -6131,7 +6131,7 @@ return; } /*}}}*/ -void Tria::GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/ +void Tria::GiaDeflection(Vector* wg,Vector* dwgdt, Matlitho* litho, IssmDouble* x, IssmDouble* y){/*{{{*/ IssmDouble xyz_list[NUMVERTICES][3]; @@ -6237,8 +6237,8 @@ } /*}}}*/ -//new logic -void Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/ +//new code +void Tria::SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze){ /*{{{*/ /*diverse:*/ int nel; @@ -6427,7 +6427,7 @@ return; } /*}}}*/ -void Tria::SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea){ /*{{{*/ +void Tria::SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks){ /*{{{*/ /*diverse:*/ IssmDouble area; @@ -6435,7 +6435,6 @@ IssmDouble phi_water=1.0; //WARNING: do not touch this, default is entire elemnt contributes barystatic IssmDouble I,W,BP; //change in ice thickness, water column or bottom pressure (Farrel and Clarke, Equ. 4) bool notfullygrounded=false; - bool scaleoceanarea= false; bool computerigid= false; int glfraction=1; int npartice,nparthydro,npartocean; @@ -6468,9 +6467,6 @@ #ifdef _ISSM_DEBUG_ constant=0; this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); #endif - /*barystatic_contribution[0]+=0; - barystatic_contribution[1]+=0; - barystatic_contribution[2]+=0;*/ return; } } @@ -6482,9 +6478,6 @@ #ifdef _ISSM_DEBUG_ this->AddInput(SealevelBarystaticMaskEnum,&constant,P0Enum); #endif - /*barystatic_contribution[0]+=0; - barystatic_contribution[1]+=0; - barystatic_contribution[2]+=0;*/ return; } } @@ -6493,9 +6486,6 @@ * hydrology:*/ if(!isice && !ishydro){ if(!masks->isoceanin[this->lid]){ - /*barystatic_contribution[0]+=0; - barystatic_contribution[1]+=0; - barystatic_contribution[2]+=0;*/ return; } @@ -6508,7 +6498,6 @@ rho_ice=FindParam(MaterialsRhoIceEnum); rho_water=FindParam(MaterialsRhoSeawaterEnum); this->parameters->FindParam(&computerigid,SolidearthSettingsRigidEnum); - this->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); this->parameters->FindParam(&glfraction,SolidearthSettingsGlfractionEnum); this->parameters->FindParam(&npartice,SolidearthNpartIceEnum); this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum); @@ -6517,9 +6506,9 @@ this->Element::GetInputValue(&area,AreaEnum); /*Deal with ice loads if we are on grounded ice:*/ - if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){ /*{{{*/ + if(masks->isiceonly[this->lid] && !masks->isfullyfloating[this->lid]){ - /*Compute fraction of the element that is grounded: */ + /*Compute fraction of the element that is grounded: {{{*/ if(notfullygrounded){ IssmDouble xyz_list[NUMVERTICES][3]; ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); @@ -6531,6 +6520,7 @@ #endif } else phi_ice=1.0; + /*}}}*/ /*Inform mask: */ constant+=100; //1 for ice. @@ -6568,15 +6558,13 @@ } /*}}}*/ - /*Compute barystatic contribution:*/ - _assert_(oceanarea>0.); - if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 - bslcice = rho_ice*area*I/(oceanarea*rho_water); + /*Compute barystatic contribution in kg:*/ + bslcice = rho_ice*area*I; _assert_(!xIsNan(bslcice)); /*Transfer thickness change into kg/m^2:*/ I=I*rho_ice*phi_ice; - } /*}}}*/ + } /*Deal with water loads if we are on ground:*/ if(!masks->isfullyfloating[this->lid]){ @@ -6586,10 +6574,8 @@ if (!deltathickness_input)_error_("DeltaTwsEnum input needed to compute sea level change!"); deltathickness_input->GetInputAverage(&W); - /*Compute barystatic component:*/ - _assert_(oceanarea>0.); - if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 - bslchydro = rho_freshwater*area*phi_water*W/(oceanarea*rho_water); + /*Compute barystatic component in kg:*/ + bslchydro = rho_freshwater*area*phi_water*W; _assert_(!xIsNan(bslchydro)); /*convert from m to kg/m^2:*/ @@ -6604,8 +6590,8 @@ if (!bottompressure_change_input)_error_("DeltaBottomPressureEnum pressure input needed to compute sea level change fingerprint!"); bottompressure_change_input->GetInputAverage(&BP); - /*Compute barystatic component:*/ - bslcbp = rho_water*area*BP/(oceanarea*rho_water); + /*Compute barystatic component in kg:*/ + bslcbp = rho_water*area*BP; /*convert from m to kg/m^2:*/ BP=BP*rho_water; @@ -6612,28 +6598,19 @@ } /*Plug all loads into total load vector:*/ - localloads[this->lid]+=I+W+BP; + loads->SetValue(this->sid,I+W+BP,INS_VAL); - /*Plug bslcice into barystatic contribution vector:*/ - if(barystatic_contribution_onpartition){ - int idi=reCast(partition[this->Sid()])+1; - int idj=0; - idj=0;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcice,ADD_VAL); - idj=1;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslchydro,ADD_VAL); - idj=2;barystatic_contribution_onpartition->SetValues(1,&idi,1,&idj,&bslcbp,ADD_VAL); - } + /*Keep track of barystatic contributions:*/ + barycontrib->Set(this->Sid(),bslcice,bslchydro,bslcbp); - barystatic_contribution[0]+=bslcice; - barystatic_contribution[1]+=bslchydro; - barystatic_contribution[2]+=bslcbp; - } /*}}}*/ -IssmDouble Tria::SealevelchangeConvolution(IssmDouble* loads){ /*{{{*/ +void Tria::SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks){ /*{{{*/ /*sal green function:*/ IssmDouble* G=NULL; - IssmDouble Sealevel[NUMVERTICES]={0,0,0}; + IssmDouble SealevelGRD[NUMVERTICES]={0,0,0}; + IssmDouble oceanaverage,oceanarea=0; bool sal = false; int size; @@ -6646,15 +6623,115 @@ for(int i=0;iAddInput(SealevelRSLEnum,Sealevel,P1Enum); + /*compute ocean average over element:*/ + OceanAverageOptim(&oceanaverage,&oceanarea,SealevelGRD,masks); + + /*add ocean average in the global sealevelloads vector:*/ + sealevelloads->SetValue(this->sid,oceanaverage,INS_VAL); + oceanareas->SetValue(this->sid,oceanarea,INS_VAL); + + return; +} /*}}}*/ +void Tria::SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks){ /*{{{*/ + + bool converged=false; + IssmDouble SealevelGRD[3]={0,0,0}; + IssmDouble oceanaverage,oceanarea=0; + int nel; + bool sal = false; + IssmDouble* G=NULL; + int size; + + this->parameters->FindParam(&nel,MeshNumberofelementsEnum); + this->parameters->FindParam(&sal,SolidearthSettingsRigidEnum); + + if(sal){ + this->inputs->GetArrayPtr(SealevelchangeGEnum,this->lid,&G,&size); + + for(int i=0;iSetValue(this->sid,oceanaverage,INS_VAL); } /*}}}*/ +void Tria::OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks){ /*{{{*/ + IssmDouble phi=1.0; + bool iscoastline=false; + IssmDouble area; + IssmDouble Sg_avg=0; //output + + /*Do we have an ocean?:*/ + if(!masks->isoceanin[this->lid]){ + *poceanarea=0; + *poceanaverage=0; + } + + /*Do we have a coastline?:*/ + if(!masks->isfullyfloating[this->lid])iscoastline=true; + + if(iscoastline){ + IssmDouble xyz_list[NUMVERTICES][3]; + ::GetVerticesCoordinates(&xyz_list[0][0],vertices,NUMVERTICES); + phi=1.0-this->GetGroundedPortion(&xyz_list[0][0]); + } + + /*Get area of element:*/ + this->Element::GetInputValue(&area,AreaEnum); + + /*Average over ocean if there is no coastline, area of the ocean + *is the areaa of the element:*/ + if(!iscoastline){ + + /*Average Sg over vertices:*/ + for(int i=0;iGetGroundedPart(&point1,&fraction1,&fraction2,&mainlyfloating); + //!mainlyfloating so that the integration is done on the ocean (and not the grounded) part + Gauss* gauss = this->NewGauss(point1,fraction1,fraction2,!mainlyfloating,2); + + /* Start looping on the number of gaussian points and average over these gaussian points: */ + total_weight=0; + Sg_avg=0; + while(gauss->next()){ + IssmDouble Sg_gauss=0; + TriaRef::GetInputValue(&Sg_gauss, Sg, gauss,P1Enum); + Sg_avg+=Sg_gauss*gauss->weight; + total_weight+=gauss->weight; + } + Sg_avg=Sg_avg/total_weight; + delete gauss; + + *poceanaverage=Sg_avg; + *poceanarea=area; + return; + +} +/*}}}*/ #endif #ifdef _HAVE_DAKOTA_ Index: ../trunk-jpl/src/c/classes/Elements/Tria.h =================================================================== --- ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 26109) +++ ../trunk-jpl/src/c/classes/Elements/Tria.h (revision 26110) @@ -161,20 +161,22 @@ void EsaGeodetic3D(Vector* pUp,Vector* pNorth,Vector* pEast,IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius,IssmDouble* xx,IssmDouble* yy,IssmDouble* zz); #endif #ifdef _HAVE_SEALEVELCHANGE_ + void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks); + void SealevelchangeGeometry(IssmDouble* lat, IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze); + IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea); + IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea); + void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks); + void SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks); + void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks); + void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); IssmDouble OceanAverage(IssmDouble* Sg, SealevelMasks* masks); - void SealevelchangeMomentOfInertia(IssmDouble* dI_list,IssmDouble* Sg_old,SealevelMasks* masks); - void SetSealevelMasks(SealevelMasks* masks); - void SealevelchangeGeometry(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze); - void SealevelchangeGeometryOptim(IssmDouble* latitude,IssmDouble* longitude,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze); - - IssmDouble SealevelchangeBarystaticIce(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition, IssmDouble oceanarea); - IssmDouble SealevelchangeBarystaticHydro(IssmDouble* Sgi, SealevelMasks* masks, Vector* barystatic_contribution, IssmDouble* partition,IssmDouble oceanarea); - void SealevelchangeBarystaticBottomPressure(IssmDouble* Sgi,SealevelMasks* masks); - void SealevelchangeSal(IssmDouble* Sgo,IssmDouble* Sg_old,SealevelMasks* masks); - void DeformationFromSurfaceLoads(IssmDouble* Up, IssmDouble* North, IssmDouble* East, IssmDouble* Sg,SealevelMasks* masks); - void GiaDeflection(Vector* wg,Vector* dwgdt,Matlitho* litho, IssmDouble* x,IssmDouble* y); - void SealevelchangeBarystaticLoads(IssmDouble* barystatic_contribution,IssmDouble* localloads, SealevelMasks* masks, Matrix* barystatic_contribution_onpartition, IssmDouble* partition, IssmDouble oceanarea); - IssmDouble SealevelchangeConvolution(IssmDouble* loads); + void SetSealevelMasks(SealevelMasks* masks); + + void SealevelchangeGeometryOptim(IssmDouble* lat,IssmDouble* longi,IssmDouble* radius, IssmDouble* xx, IssmDouble* yy, IssmDouble* zz, IssmDouble* xxe, IssmDouble* yye, IssmDouble* zze); + void SealevelchangeBarystaticLoads(Vector* loads, BarystaticContributions* barycontrib, SealevelMasks* masks); + void SealevelchangeInitialConvolution(Vector* sealevelloads,Vector* oceanareas,IssmDouble* loads, SealevelMasks* masks); + void SealevelchangeOceanConvolution(Vector* newsealevelloads, IssmDouble* sealevelloads, IssmDouble* loads, SealevelMasks* masks); + void OceanAverageOptim(IssmDouble* poceanaverage, IssmDouble* poceanarea, IssmDouble* Sg, SealevelMasks* masks); #endif /*}}}*/ /*Tria specific routines:{{{*/ Index: ../trunk-jpl/src/c/classes/FemModel.cpp =================================================================== --- ../trunk-jpl/src/c/classes/FemModel.cpp (revision 26109) +++ ../trunk-jpl/src/c/classes/FemModel.cpp (revision 26110) @@ -4748,164 +4748,6 @@ /*}}}*/ #endif #ifdef _HAVE_SEALEVELCHANGE_ -void FemModel::SealevelchangeBarystatic(Vector* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition,SealevelMasks* masks) { /*{{{*/ - - /*serialized vectors:*/ - IssmDouble bslcice = 0.; - IssmDouble bslcice_cpu = 0.; - IssmDouble bslchydro = 0.; - IssmDouble bslchydro_cpu = 0.; - IssmDouble area = 0.; - IssmDouble oceanarea = 0.; - IssmDouble oceanarea_cpu = 0.; - int bp_compute_fingerprints= 0; - bool isoceantransport=false; - - Vector* bslcice_partition=NULL; - IssmDouble* bslcice_partition_serial=NULL; - IssmDouble* partitionice=NULL; - int npartice,nel; - - Vector* bslchydro_partition=NULL; - IssmDouble* bslchydro_partition_serial=NULL; - IssmDouble* partitionhydro=NULL; - bool istws=0; - int nparthydro; - - int npartocean; - Vector* bslcocean_partition=NULL; - IssmDouble* bslcocean_partition_serial=NULL; - IssmDouble* partitionocean=NULL; - - /*Initialize temporary vector that will be used to sum barystatic components - * on all local elements, prior to assembly:*/ - int gsize = this->nodes->NumberOfDofs(GsetEnum); - IssmDouble* RSLgi=xNewZeroInit(gsize); - int* indices=xNew(gsize); - for(int i=0;ielements->objects){ - i +=1; - Element* element = xDynamicCast(object); - element->GetInputValue(&area,AreaEnum); - if (masks->isoceanin[i]) oceanarea_cpu += area; - } - ISSM_MPI_Reduce (&oceanarea_cpu,&oceanarea,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&oceanarea,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - _assert_(oceanarea>0.); - - /*Initialize partition vectors to retrieve barystatic contributions: */ - this->parameters->FindParam(&npartice,SolidearthNpartIceEnum); - if(npartice){ - this->parameters->FindParam(&partitionice,&nel,NULL,SolidearthPartitionIceEnum); - bslcice_partition= new Vector(npartice); - } - - this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum); - if(nparthydro){ - this->parameters->FindParam(&partitionhydro,&nel,NULL,SolidearthPartitionHydroEnum); - bslchydro_partition= new Vector(nparthydro); - } - - this->parameters->FindParam(&npartocean,SolidearthNpartOceanEnum); - if(npartocean){ - this->parameters->FindParam(&partitionocean,&nel,NULL,SolidearthPartitionOceanEnum); - bslchydro_partition= new Vector(npartocean); - } - /*For later: - npartbarystatic=npartice; - if(nparthydro>npartbarystatic)npartbarystatic=nparthydro; - if(npartocean>npartbarystatic)npartbarystatic=npartocean; - bslc_partition=new Matrix(IssmDouble>(npartbarystatic,3); - - bslc_cpu[0]=0; bslc_cpu[1]=0; bslc_cpu[2]=0; - for(Object* & object : this->elements->objects){ - Element* element = xDynamicCast(object); - element->SealevelchangeBarystaticLoads(&bslc_cpu[0], localloads,masks, bslcice_partition,partitionice,oceanarea); - } - MPI Bcast localloads -> loads - - for(Object* & object : this->elements->objects){ - Element* element = xDynamicCast(object); - element->SealevelchangeConvolution(loads); - } - */ - - - - - - /*Call the barystatic sea level change core for ice : */ - bslcice_cpu=0; - for(Object* & object : this->elements->objects){ - Element* element = xDynamicCast(object); - bslcice_cpu+=element->SealevelchangeBarystaticIce(RSLgi,masks, bslcice_partition,partitionice,oceanarea); - } - - /*Call the barystatic sea level change core for hydro: */ - bslchydro_cpu=0; //make sure to initialize this, so we have a total barystatic contribution computed at 0. - this->parameters->FindParam(&istws,TransientIshydrologyEnum); - if(istws){ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - bslchydro_cpu+=element->SealevelchangeBarystaticHydro(RSLgi,masks, bslchydro_partition,partitionhydro,oceanarea); - } - } - - /*Call the barystatic sea level change core for bottom pressures: */ - this->parameters->FindParam(&bp_compute_fingerprints,SolidearthSettingsComputeBpGrdEnum); - this->parameters->FindParam(&isoceantransport,TransientIsoceantransportEnum); - if(bp_compute_fingerprints && isoceantransport){ - for(int i=0;iSize();i++){ - Element* element=xDynamicCast(elements->GetObjectByOffset(i)); - element->SealevelchangeBarystaticBottomPressure(RSLgi,masks); - } - } - - /*Plug values once and assemble: */ - pRSLgi->SetValues(gsize,indices,RSLgi,ADD_VAL); - pRSLgi->Assemble(); - - /*Sum all barystatic components from all cpus:*/ - ISSM_MPI_Reduce (&bslcice_cpu,&bslcice,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&bslcice,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - _assert_(!xIsNan(bslcice)); - - ISSM_MPI_Reduce (&bslchydro_cpu,&bslchydro,1,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,0,IssmComm::GetComm() ); - ISSM_MPI_Bcast(&bslchydro,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm()); - _assert_(!xIsNan(bslchydro)); - - /*Take care of partition vectors:*/ - if(bslcice_partition){ - bslcice_partition->Assemble(); - bslcice_partition_serial=bslcice_partition->ToMPISerial(); - } - if(bslchydro_partition){ - bslchydro_partition->Assemble(); - bslchydro_partition_serial=bslchydro_partition->ToMPISerial(); - } - - - /*Free ressources:*/ - xDelete(indices); - xDelete(RSLgi); - if(bslchydro_partition)delete bslchydro_partition; - if(bslcice_partition)delete bslcice_partition; - if(partitionhydro)xDelete(partitionhydro); - if(partitionice)xDelete(partitionice); - - /*Assign output pointers:*/ - *poceanarea = oceanarea; - *pbslcice = bslcice; - *pbslchydro = bslchydro; - *pbslc=bslchydro+bslcice; - *pbslcice_partition=bslcice_partition_serial; - *pbslchydro_partition=bslchydro_partition_serial; - -} -/*}}}*/ void FemModel::SealevelchangeSal(Vector* pRSLgo, Vector* pRSLg_old, SealevelMasks* masks, bool verboseconvolution){/*{{{*/ /*serialized vectors:*/ Index: ../trunk-jpl/src/c/classes/classes.h =================================================================== --- ../trunk-jpl/src/c/classes/classes.h (revision 26109) +++ ../trunk-jpl/src/c/classes/classes.h (revision 26110) @@ -18,6 +18,7 @@ #include "./Massfluxatgate.h" #include "./Misfit.h" #include "./SealevelMasks.h" +#include "./BarystaticContributions.h" #include "./Nodalvalue.h" #include "./Numberedcostfunction.h" #include "./Cfsurfacesquare.h" Index: ../trunk-jpl/src/c/cores/cores.h =================================================================== --- ../trunk-jpl/src/c/cores/cores.h (revision 26109) +++ ../trunk-jpl/src/c/cores/cores.h (revision 26110) @@ -80,10 +80,7 @@ void WriteLockFile(char* filename); void ResetBoundaryConditions(FemModel* femmodel, int analysis_type); void PrintBanner(void); -void TransferForcing(FemModel* femmodel,int forcingenum); -void TransferSealevel(FemModel* femmodel,int forcingenum); void EarthMassTransport(FemModel* femmodel); -void slcconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs); //solution configuration void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype); Index: ../trunk-jpl/src/c/cores/sealevelchange_core.cpp =================================================================== --- ../trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26109) +++ ../trunk-jpl/src/c/cores/sealevelchange_core.cpp (revision 26110) @@ -11,6 +11,12 @@ #include "../shared/shared.h" #include "../modules/modules.h" #include "../solutionsequences/solutionsequences.h" +/*support routines local definitions:{{{*/ +void TransferForcing(FemModel* femmodel,int forcingenum); +void TransferSealevel(FemModel* femmodel,int forcingenum); +void slcconvergence(bool* pconverged, Vector* RSLg,Vector* RSLg_old,IssmDouble eps_rel,IssmDouble eps_abs); +IssmDouble SealevelloadsOceanAverage(Vector* sealevelloads,Vector* oceanareas, IssmDouble oceanarea); +/*}}}*/ /*main cores:*/ void sealevelchange_core(FemModel* femmodel){ /*{{{*/ @@ -383,6 +389,107 @@ } }; /*}}}*/ +void grd_core_optim(FemModel* femmodel,SealevelMasks* masks) { /*{{{*/ + + /*variables:{{{*/ + int nel; + BarystaticContributions* barycontrib=NULL; + GenericParam* barycontribparam=NULL; + IssmDouble barystatic; + + Vector* loads=NULL; + IssmDouble* allloads=NULL; + Vector* sealevelloads=NULL; + Vector* oldsealevelloads=NULL; + IssmDouble sealevelloadsaverage; + IssmDouble* allsealevelloads=NULL; + Vector* oceanareas=NULL; + IssmDouble oceanarea; + bool scaleoceanarea=false; + IssmDouble rho_water; + + IssmDouble eps_rel; + IssmDouble eps_abs; + int step; + IssmDouble time; + + IssmDouble cumbslc; + IssmDouble cumbslcice; + IssmDouble cumbslchydro; + /*}}}*/ + + /*retrieve parameters: */ + femmodel->parameters->FindParam(&scaleoceanarea,SolidearthSettingsOceanAreaScalingEnum); + barycontribparam = xDynamicCast*>(femmodel->parameters->FindParamObject(BarystaticContributionsEnum)); + barycontrib=barycontribparam->GetParameterValue(); + femmodel->parameters->FindParam(&rho_water,MaterialsRhoSeawaterEnum); + femmodel->parameters->FindParam(&eps_rel,SolidearthSettingsReltolEnum); + femmodel->parameters->FindParam(&eps_abs,SolidearthSettingsAbstolEnum); + + /*initialize matrices and vectors:*/ + femmodel->parameters->FindParam(&nel,MeshNumberofelementsEnum); + loads=new Vector(nel); + sealevelloads=new Vector(nel); + oceanareas=new Vector(nel); + + /*buildup loads: */ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeBarystaticLoads(loads, barycontrib,masks); + } + + //Communicate loads from local to global: + loads->Assemble(); allloads=loads->ToMPISerial(); + + /*convolve loads:*/ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeInitialConvolution(sealevelloads,oceanareas,allloads,masks); + } + + //Get ocean area: + oceanareas->Assemble(); oceanareas->Sum(&oceanarea); _assert_(oceanarea>0.); + if(scaleoceanarea) oceanarea=3.619e+14; // use true ocean area, m^2 + + //Get sea level loads ocean average: + sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea); + + //substract ocean average and barystatic contributionfrom sea level loads: + barystatic=barycontrib->Total()/oceanarea/rho_water; + sealevelloads->Shift(-sealevelloadsaverage+barystatic); + allsealevelloads=sealevelloads->ToMPISerial(); + + bool converged=false; + for(;;){ + + oldsealevelloads=sealevelloads->Duplicate(); + + /*convolve load and sealevel loads on oceans:*/ + for(Object* & object : femmodel->elements->objects){ + Element* element = xDynamicCast(object); + element->SealevelchangeOceanConvolution(sealevelloads, allsealevelloads, allloads,masks); + } + sealevelloads->Assemble(); + + //substract ocean average and barystatic contribution + sealevelloadsaverage=SealevelloadsOceanAverage(sealevelloads,oceanareas,oceanarea); + sealevelloads->Shift(-sealevelloadsaverage+barystatic); + + //broadcast loads + allsealevelloads=sealevelloads->ToMPISerial(); + + //convergence? + slcconvergence(&converged,sealevelloads,oldsealevelloads,eps_rel,eps_abs); + if (converged)break; + } + + /*cumulate barystatic contributions and save to results: */ + barycontrib->Cumulate(femmodel->parameters); + barycontrib->Save(femmodel->results,femmodel->parameters,oceanarea); + +} +/*}}}*/ + //Geometry: SealevelMasks* sealevel_masks(FemModel* femmodel) { /*{{{*/ @@ -1061,3 +1168,11 @@ *pconverged=converged; } /*}}}*/ +IssmDouble SealevelloadsOceanAverage(Vector* sealevelloads,Vector* oceanareas, IssmDouble oceanarea){ /*{{{*/ + IssmDouble sealevelloadsaverage; + Vector* sealevelloadsvolume=sealevelloads->Duplicate(); + sealevelloadsvolume->PointwiseMult(sealevelloads,oceanareas); + sealevelloadsvolume->Sum(&sealevelloadsaverage); + delete sealevelloadsvolume; + return sealevelloadsaverage/oceanarea; +} /*}}}*/ Index: ../trunk-jpl/src/c/shared/Enum/Enum.vim =================================================================== --- ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26109) +++ ../trunk-jpl/src/c/shared/Enum/Enum.vim (revision 26110) @@ -746,6 +746,7 @@ syn keyword cConstant BslcEnum syn keyword cConstant BslcIceEnum syn keyword cConstant BslcHydroEnum +syn keyword cConstant BslcOceanEnum syn keyword cConstant BslcRateEnum syn keyword cConstant GmtslcEnum syn keyword cConstant SealevelRSLBarystaticEnum @@ -1070,6 +1071,7 @@ syn keyword cConstant BalancevelocitySolutionEnum syn keyword cConstant BasalforcingsIsmip6Enum syn keyword cConstant BasalforcingsPicoEnum +syn keyword cConstant BarystaticContributionsEnum syn keyword cConstant BeckmannGoosseFloatingMeltRateEnum syn keyword cConstant BedSlopeSolutionEnum syn keyword cConstant BoolExternalResultEnum @@ -1432,6 +1434,7 @@ syn keyword cType AmrBamg syn keyword cType AmrNeopz syn keyword cType ArrayInput +syn keyword cType BarystaticContributions syn keyword cType BoolInput syn keyword cType BoolParam syn keyword cType Cfdragcoeffabsgrad Index: ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26109) +++ ../trunk-jpl/src/c/shared/Enum/EnumDefinitions.h (revision 26110) @@ -742,6 +742,7 @@ BslcEnum, BslcIceEnum, BslcHydroEnum, + BslcOceanEnum, BslcRateEnum, GmtslcEnum, SealevelRSLBarystaticEnum, @@ -1069,6 +1070,7 @@ BalancevelocitySolutionEnum, BasalforcingsIsmip6Enum, BasalforcingsPicoEnum, + BarystaticContributionsEnum, BeckmannGoosseFloatingMeltRateEnum, BedSlopeSolutionEnum, BoolExternalResultEnum, Index: ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26109) +++ ../trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp (revision 26110) @@ -748,6 +748,7 @@ case BslcEnum : return "Bslc"; case BslcIceEnum : return "BslcIce"; case BslcHydroEnum : return "BslcHydro"; + case BslcOceanEnum : return "BslcOcean"; case BslcRateEnum : return "BslcRate"; case GmtslcEnum : return "Gmtslc"; case SealevelRSLBarystaticEnum : return "SealevelRSLBarystatic"; @@ -1072,6 +1073,7 @@ case BalancevelocitySolutionEnum : return "BalancevelocitySolution"; case BasalforcingsIsmip6Enum : return "BasalforcingsIsmip6"; case BasalforcingsPicoEnum : return "BasalforcingsPico"; + case BarystaticContributionsEnum : return "BarystaticContributions"; case BeckmannGoosseFloatingMeltRateEnum : return "BeckmannGoosseFloatingMeltRate"; case BedSlopeSolutionEnum : return "BedSlopeSolution"; case BoolExternalResultEnum : return "BoolExternalResult"; Index: ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp =================================================================== --- ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26109) +++ ../trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp (revision 26110) @@ -766,6 +766,7 @@ else if (strcmp(name,"Bslc")==0) return BslcEnum; else if (strcmp(name,"BslcIce")==0) return BslcIceEnum; else if (strcmp(name,"BslcHydro")==0) return BslcHydroEnum; + else if (strcmp(name,"BslcOcean")==0) return BslcOceanEnum; else if (strcmp(name,"BslcRate")==0) return BslcRateEnum; else if (strcmp(name,"Gmtslc")==0) return GmtslcEnum; else if (strcmp(name,"SealevelRSLBarystatic")==0) return SealevelRSLBarystaticEnum; @@ -873,11 +874,11 @@ else if (strcmp(name,"SmbSmbCorr")==0) return SmbSmbCorrEnum; else if (strcmp(name,"SmbSmbref")==0) return SmbSmbrefEnum; else if (strcmp(name,"SmbSzaValue")==0) return SmbSzaValueEnum; - else if (strcmp(name,"SmbT")==0) return SmbTEnum; else stage=8; } if(stage==8){ - if (strcmp(name,"SmbTa")==0) return SmbTaEnum; + if (strcmp(name,"SmbT")==0) return SmbTEnum; + else if (strcmp(name,"SmbTa")==0) return SmbTaEnum; else if (strcmp(name,"SmbTeValue")==0) return SmbTeValueEnum; else if (strcmp(name,"SmbTemperaturesAnomaly")==0) return SmbTemperaturesAnomalyEnum; else if (strcmp(name,"SmbTemperaturesLgm")==0) return SmbTemperaturesLgmEnum; @@ -996,11 +997,11 @@ else if (strcmp(name,"Outputdefinition31")==0) return Outputdefinition31Enum; else if (strcmp(name,"Outputdefinition32")==0) return Outputdefinition32Enum; else if (strcmp(name,"Outputdefinition33")==0) return Outputdefinition33Enum; - else if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum; else stage=9; } if(stage==9){ - if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; + if (strcmp(name,"Outputdefinition34")==0) return Outputdefinition34Enum; + else if (strcmp(name,"Outputdefinition35")==0) return Outputdefinition35Enum; else if (strcmp(name,"Outputdefinition36")==0) return Outputdefinition36Enum; else if (strcmp(name,"Outputdefinition37")==0) return Outputdefinition37Enum; else if (strcmp(name,"Outputdefinition38")==0) return Outputdefinition38Enum; @@ -1096,6 +1097,7 @@ else if (strcmp(name,"BalancevelocitySolution")==0) return BalancevelocitySolutionEnum; else if (strcmp(name,"BasalforcingsIsmip6")==0) return BasalforcingsIsmip6Enum; else if (strcmp(name,"BasalforcingsPico")==0) return BasalforcingsPicoEnum; + else if (strcmp(name,"BarystaticContributions")==0) return BarystaticContributionsEnum; else if (strcmp(name,"BeckmannGoosseFloatingMeltRate")==0) return BeckmannGoosseFloatingMeltRateEnum; else if (strcmp(name,"BedSlopeSolution")==0) return BedSlopeSolutionEnum; else if (strcmp(name,"BoolExternalResult")==0) return BoolExternalResultEnum; @@ -1118,12 +1120,12 @@ else if (strcmp(name,"ChannelDischarge")==0) return ChannelDischargeEnum; else if (strcmp(name,"Closed")==0) return ClosedEnum; else if (strcmp(name,"Colinear")==0) return ColinearEnum; - else if (strcmp(name,"Constraints")==0) return ConstraintsEnum; - else if (strcmp(name,"Contact")==0) return ContactEnum; else stage=10; } if(stage==10){ - if (strcmp(name,"Contour")==0) return ContourEnum; + if (strcmp(name,"Constraints")==0) return ConstraintsEnum; + else if (strcmp(name,"Contact")==0) return ContactEnum; + else if (strcmp(name,"Contour")==0) return ContourEnum; else if (strcmp(name,"Contours")==0) return ContoursEnum; else if (strcmp(name,"ControlInput")==0) return ControlInputEnum; else if (strcmp(name,"ControlInputGrad")==0) return ControlInputGradEnum; @@ -1241,12 +1243,12 @@ else if (strcmp(name,"L2ProjectionBaseAnalysis")==0) return L2ProjectionBaseAnalysisEnum; else if (strcmp(name,"L2ProjectionEPLAnalysis")==0) return L2ProjectionEPLAnalysisEnum; else if (strcmp(name,"LACrouzeixRaviart")==0) return LACrouzeixRaviartEnum; - else if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; - else if (strcmp(name,"LambdaS")==0) return LambdaSEnum; else stage=11; } if(stage==11){ - if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; + if (strcmp(name,"LATaylorHood")==0) return LATaylorHoodEnum; + else if (strcmp(name,"LambdaS")==0) return LambdaSEnum; + else if (strcmp(name,"LevelsetAnalysis")==0) return LevelsetAnalysisEnum; else if (strcmp(name,"LevelsetfunctionPicard")==0) return LevelsetfunctionPicardEnum; else if (strcmp(name,"LinearFloatingMeltRate")==0) return LinearFloatingMeltRateEnum; else if (strcmp(name,"LliboutryDuval")==0) return LliboutryDuvalEnum; @@ -1364,12 +1366,12 @@ else if (strcmp(name,"SMBpdd")==0) return SMBpddEnum; else if (strcmp(name,"SMBpddSicopolis")==0) return SMBpddSicopolisEnum; else if (strcmp(name,"SMBsemic")==0) return SMBsemicEnum; - else if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; - else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; else stage=12; } if(stage==12){ - if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; + if (strcmp(name,"SSAApproximation")==0) return SSAApproximationEnum; + else if (strcmp(name,"SSAFSApproximation")==0) return SSAFSApproximationEnum; + else if (strcmp(name,"SSAHOApproximation")==0) return SSAHOApproximationEnum; else if (strcmp(name,"Scaled")==0) return ScaledEnum; else if (strcmp(name,"SealevelAbsolute")==0) return SealevelAbsoluteEnum; else if (strcmp(name,"SealevelEmotion")==0) return SealevelEmotionEnum; Index: ../trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h =================================================================== --- ../trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h (revision 26109) +++ ../trunk-jpl/src/c/toolkits/issm/IssmSeqVec.h (revision 26110) @@ -323,6 +323,14 @@ } /*}}}*/ + void Sum(doubletype* pvalue){/*{{{*/ + + doubletype value=0; + int i; + for(i=0;iM;i++)value+=this->vector[i]; + *pvalue=value; + } + /*}}}*/ }; #endif //#ifndef _ISSM_SEQ_VEC_H_ Index: ../trunk-jpl/src/c/toolkits/objects/Matrix.h =================================================================== --- ../trunk-jpl/src/c/toolkits/objects/Matrix.h (revision 26109) +++ ../trunk-jpl/src/c/toolkits/objects/Matrix.h (revision 26110) @@ -281,6 +281,22 @@ return output; } /*}}}*/ + doubletype* ToMPISerial(void){/*{{{*/ + + doubletype* output=NULL; + + if(type==PetscMatType){ + #ifdef _HAVE_PETSC_ + output=this->pmatrix->ToMPISerial(); + #endif + } + else{ + _error_("not implemented yet!"); + } + + return output; + } + /*}}}*/ void SetValues(int m,int* idxm,int n,int* idxn,IssmDouble* values,InsMode mode){/*{{{*/ if(type==PetscMatType){ @@ -306,10 +322,8 @@ } /*}}}*/ - /* - * sets all values to 0 but keeps the structure of a sparse matrix - */ void SetZero(void) {/*{{{*/ + // sets all values to 0 but keeps the structure of a sparse matrix if(type==PetscMatType){ #ifdef _HAVE_PETSC_ this->pmatrix->SetZero(); Index: ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp =================================================================== --- ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp (revision 26109) +++ ../trunk-jpl/src/c/toolkits/petsc/objects/PetscVec.cpp (revision 26110) @@ -248,6 +248,13 @@ } /*}}}*/ +void PetscVec::Sum(IssmDouble* pvalue){/*{{{*/ + + _assert_(this->vector); + VecSum(this->vector,pvalue); + +} +/*}}}*/ IssmDouble PetscVec::Dot(PetscVec* input){/*{{{*/ IssmDouble dot; Index: ../trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h =================================================================== --- ../trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h (revision 26109) +++ ../trunk-jpl/src/c/toolkits/petsc/patches/petscpatches.h (revision 26110) @@ -31,6 +31,7 @@ void PetscOptionsDetermineSolverType(int* psolver_type); void MatMultPatch(Mat A,Vec X, Vec AX,ISSM_MPI_Comm comm); void MatToSerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm); +void MatToMPISerial(double** poutmatrix,Mat matrix,ISSM_MPI_Comm comm); Vec SerialToVec(double* vector,int vector_size); InsertMode ISSMToPetscInsertMode(InsMode mode); NormType ISSMToPetscNormMode(NormMode mode); Index: ../trunk-jpl/test/NightlyRun/test2004.m =================================================================== --- ../trunk-jpl/test/NightlyRun/test2004.m (revision 26109) +++ ../trunk-jpl/test/NightlyRun/test2004.m (revision 26110) @@ -133,6 +133,10 @@ disp(' reading bedrock'); md.geometry.bed=-ones(md.mesh.numberofvertices,1); + md.geometry.base=md.geometry.bed; + md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1); + md.geometry.surface=md.geometry.bed+md.geometry.thickness; + end % }}} %Slc: {{{ if bas.iscontinentany('antarctica'), @@ -166,11 +170,11 @@ md.masstransport.spcthickness=mean(delHAIS(md.mesh.elements),2)/100; end - md.solidearth.initialsealevel=zeros(md.mesh.numberofvertices,1); + md.initialization.sealevel=zeros(md.mesh.numberofvertices,1); - md.dsl.global_average_thermosteric_sea_level_change=[0;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); + md.dsl.global_average_thermosteric_sea_level=[0;0]; + md.dsl.sea_surface_height_above_geoid=zeros(md.mesh.numberofvertices+1,1); + md.dsl.sea_water_pressure_at_sea_floor=zeros(md.mesh.numberofvertices+1,1); end %}}} % material properties: {{{ @@ -312,6 +316,9 @@ %geometry: {{{ di=md.materials.rho_ice/md.materials.rho_water; md.geometry.bed=-ones(md.mesh.numberofvertices,1); + md.geometry.base=md.geometry.bed; + md.geometry.thickness=1000*ones(md.mesh.numberofvertices,1); + md.geometry.surface=md.geometry.bed+md.geometry.thickness; % }}} %materials: {{{ md.materials=materials('hydro'); @@ -345,6 +352,9 @@ sl.transfer('mask.ice_levelset'); sl.transfer('mask.ocean_levelset'); sl.transfer('geometry.bed'); +sl.transfer('geometry.surface'); +sl.transfer('geometry.thickness'); +sl.transfer('geometry.base'); sl.transfer('mesh.lat'); sl.transfer('mesh.long'); sl.transfer('masstransport.spcthickness'); @@ -399,6 +409,16 @@ md.transient.ismasstransport=1; md.transient.isslc=1; +%Initializations: +md.basalforcings.groundedice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.basalforcings.floatingice_melting_rate=zeros(md.mesh.numberofvertices,1); +md.initialization.vx=zeros(md.mesh.numberofvertices,1); +md.initialization.vy=zeros(md.mesh.numberofvertices,1); +md.initialization.sealevel=zeros(md.mesh.numberofvertices,1); +md.initialization.bottompressure=zeros(md.mesh.numberofvertices,1); +md.initialization.dsl=zeros(md.mesh.numberofvertices,1); +md.initialization.str=0; +md.smb.mass_balance=zeros(md.mesh.numberofvertices,1); %max number of iterations reverted back to 10 (i.e. the original default value) md.solidearth.settings.maxiter=10;