Changeset 26114

Timestamp:

03/18/21 20:40:04 (4 years ago)

Author:

Eric.Larour

Message:

CHG: had taken out these routines which are needed!

Location:

issm/trunk-jpl/src/c/classes

Files:

• FemModel.cpp (modified) (1 diff)
• FemModel.h (modified) (1 diff)

issm/trunk-jpl/src/c/classes/FemModel.cpp

@@ -4749 +4749 @@
 #endif
 #ifdef _HAVE_SEALEVELCHANGE_
+void FemModel::SealevelchangeBarystatic(Vector<IssmDouble>* 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<IssmDouble>* bslcice_partition=NULL;
+        IssmDouble* bslcice_partition_serial=NULL;
+        IssmDouble* partitionice=NULL;
+        int npartice,nel;
+
+        Vector<IssmDouble>* bslchydro_partition=NULL;
+        IssmDouble* bslchydro_partition_serial=NULL;
+        IssmDouble* partitionhydro=NULL;
+        bool istws=0;
+        int nparthydro;
+                
+        int npartocean;
+        Vector<IssmDouble>* 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<IssmDouble>(gsize);
+        int* indices=xNew<int>(gsize);
+   for(int i=0;i<gsize;i++) indices[i]=i;
+
+        /*First, figure out the area of the ocean, which is needed to compute the barystatic component: */
+        int i = -1;
+        for(Object* & object : this->elements->objects){
+                i +=1;
+                Element* element = xDynamicCast<Element*>(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<IssmDouble>(npartice);
+        }
+
+        this->parameters->FindParam(&nparthydro,SolidearthNpartHydroEnum);
+        if(nparthydro){
+                this->parameters->FindParam(&partitionhydro,&nel,NULL,SolidearthPartitionHydroEnum);
+                bslchydro_partition= new Vector<IssmDouble>(nparthydro);
+        }
+
+        this->parameters->FindParam(&npartocean,SolidearthNpartOceanEnum);
+        if(npartocean){
+                this->parameters->FindParam(&partitionocean,&nel,NULL,SolidearthPartitionOceanEnum);
+                bslchydro_partition= new Vector<IssmDouble>(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<Element*>(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<Element*>(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<Element*>(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;i<elements->Size();i++){
+                        Element* element=xDynamicCast<Element*>(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;i<elements->Size();i++){
+                        Element* element=xDynamicCast<Element*>(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<IssmDouble>(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<IssmDouble>(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<int>(indices);
+        xDelete<IssmDouble>(RSLgi);
+        if(bslchydro_partition)delete bslchydro_partition;
+        if(bslcice_partition)delete bslcice_partition;
+        if(partitionhydro)xDelete<IssmDouble>(partitionhydro);
+        if(partitionice)xDelete<IssmDouble>(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<IssmDouble>* pRSLgo, Vector<IssmDouble>* pRSLg_old,  SealevelMasks* masks, bool verboseconvolution){/*{{{*/
 

issm/trunk-jpl/src/c/classes/FemModel.h

@@ -164 +164 @@
                 #endif
                 #ifdef _HAVE_SEALEVELCHANGE_
+                void SealevelchangeBarystatic(Vector<IssmDouble>* pRSLgi, IssmDouble* poceanarea, IssmDouble* pbslc,IssmDouble* pbslcice, IssmDouble* pbslchydro, IssmDouble** pbslcice_partition,IssmDouble** pbslchydro_partition, SealevelMasks* masks); 
                 void SealevelchangeSal(Vector<IssmDouble>* pSgo, Vector<IssmDouble>* pSg_old,  SealevelMasks* masks,bool verboseconvolution);
                 void SealevelchangeRotationalFeedback(Vector<IssmDouble>* pRSLgo_rot, Vector<IssmDouble>* pRSLg_old, IssmDouble* pIxz, IssmDouble* pIyz, IssmDouble* pIzz, SealevelMasks* masks);