source: issm/oecreview/Archive/19101-20495/ISSM-20135-20136.diff

../trunk-jpl/src/c/analyses/SealevelriseAnalysis.cpp

@@ -51 +51 @@
         int     numoutputs;
         char**  requestedoutputs = NULL;
 
+        /*transition vectors: */
+        IssmDouble **transitions    = NULL;
+        int         *transitions_M    = NULL;
+        int         *transitions_N    = NULL;
+        int          ntransitions;
+
         /*some constant parameters: */
         parameters->AddObject(iomodel->CopyConstantObject(SealevelriseReltolEnum));
         parameters->AddObject(iomodel->CopyConstantObject(SealevelriseAbstolEnum));
@@ -125 +131 @@
                 xDelete<IssmDouble>(G_elastic_local);
         }
         
+        /*Transitions: */
+        iomodel->FetchData(&transitions,&transitions_M,&transitions_N,&ntransitions,SealevelriseTransitionsEnum);
+        if(transitions){
+                parameters->AddObject(new DoubleMatArrayParam(SealevelriseTransitionsEnum,transitions,ntransitions,transitions_M,transitions_N));
+
+                for(int i=0;i<ntransitions;i++){
+                        IssmDouble* transition=transitions[i];
+                        xDelete<IssmDouble>(transition);
+                }
+                xDelete<IssmDouble*>(transitions);
+                xDelete<IssmDouble*>(transitions_M);
+                xDelete<IssmDouble*>(transitions_N);
+        }
+
         /*Requested outputs*/
         iomodel->FetchData(&requestedoutputs,&numoutputs,SealevelriseRequestedOutputsEnum);
         parameters->AddObject(new IntParam(SealevelriseNumRequestedOutputsEnum,numoutputs));

../trunk-jpl/src/c/cores/transient_core.cpp

@@ -20 +20 @@
 
         /*parameters: */
         IssmDouble starttime,finaltime,dt,yts;
-        bool       isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,islevelset,isdamageevolution,ishydrology,iscalving;
+        bool       isstressbalance,ismasstransport,issmb,isFS,isthermal,isgroundingline,isgia,isslr,iscoupler,islevelset,isdamageevolution,ishydrology,iscalving;
         bool       save_results,dakota_analysis;
         bool       time_adapt;
         int        output_frequency;
@@ -50 +50 @@
         femmodel->parameters->FindParam(&issmb,TransientIssmbEnum);
         femmodel->parameters->FindParam(&isthermal,TransientIsthermalEnum);
         femmodel->parameters->FindParam(&isgia,TransientIsgiaEnum);
+        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
+        femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
         femmodel->parameters->FindParam(&isgroundingline,TransientIsgroundinglineEnum);
         femmodel->parameters->FindParam(&islevelset,TransientIslevelsetEnum);
         femmodel->parameters->FindParam(&isdamageevolution,TransientIsdamageevolutionEnum);
@@ -140 +142 @@
                         #endif
                 }
 
+                /*Sea level rise: */
+                if(isslr | iscoupler) sealevelrise_core(femmodel);
+
                 /*unload results*/
                 if(VerboseSolution()) _printf0_("   computing requested outputs\n");
                 femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs,save_results);

../trunk-jpl/src/c/cores/cores.h

@@ -49 +49 @@
 void smb_core(FemModel* femmodel);
 void damage_core(FemModel* femmodel);
 void sealevelrise_core(FemModel* femmodel);
-Vector<IssmDouble> * sealevelrise_core_eustatic(FemModel* femmodel);
-Vector<IssmDouble> * sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic);
+Vector<IssmDouble>* sealevelrise_core_eustatic(FemModel* femmodel);
+Vector<IssmDouble>* sealevelrise_core_noneustatic(FemModel* femmodel,Vector<IssmDouble>* Sg_eustatic);
 IssmDouble objectivefunction(IssmDouble search_scalar,FemModel* femmodel);
 
 //optimization
@@ -61 +61 @@
 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);
 
 //solution configuration
 void CorePointerFromSolutionEnum(void (**psolutioncore)(FemModel*),Parameters* parameters,int solutiontype);

../trunk-jpl/src/c/cores/sealevelrise_core.cpp

@@ -9 +9 @@
 #include "../modules/modules.h"
 #include "../solutionsequences/solutionsequences.h"
 
-void sealevelrise_core(FemModel* femmodel){
+void sealevelrise_core(FemModel* femmodel){ /*{{{*/
 
         Vector<IssmDouble> *Sg    = NULL;
         Vector<IssmDouble> *Sg_eustatic    = NULL; 
-        bool save_results;
+        bool save_results,isslr,iscoupler;
         int configuration_type;
+        int solution_type;
         int        numoutputs        = 0;
         char     **requested_outputs = NULL;
-
-        if(VerboseSolution()) _printf0_("   computing sea level rise\n");
-
+        
         /*Recover some parameters: */
         femmodel->parameters->FindParam(&configuration_type,ConfigurationTypeEnum);
+        femmodel->parameters->FindParam(&solution_type,SolutionTypeEnum);
         femmodel->parameters->FindParam(&save_results,SaveResultsEnum);
-        femmodel->parameters->FindParam(&numoutputs,SealevelriseNumRequestedOutputsEnum);
-        if(numoutputs) femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
+        femmodel->parameters->FindParam(&isslr,TransientIsslrEnum);
+        femmodel->parameters->FindParam(&iscoupler,TransientIscouplerEnum);
 
+        /*several cases here, depending on value of iscoupler and isslr: 
+        ( !iscoupler & !isslr)       we are not interested in being here :) 
+        ( !iscoupler & isslr)        we are running in uncoupled mode
+        ( iscoupler & isslr)         we are running in coupled mode, and better be earth
+        ( iscoupler & !isslr)        we are running in coupled mode, and better be an ice cap
+        */
+
+        /*early return: */
+        if( !iscoupler & !isslr) return;  //we are not interested in being here :) 
+
+        /*In what follows we assume we are all running slr, either in coupled, or uncoupled mode:*/
+        if(VerboseSolution()) _printf0_("   computing sea level rise\n");
+
         /*set configuration: */
-        femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
+        if(isslr)femmodel->SetCurrentConfiguration(SealevelriseAnalysisEnum);
 
-        /*call two sub cores:*/
-        Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS.
+        /*transfer deltathickness forcing from ice caps to earth model: */
+        if(iscoupler) TransferForcing(femmodel,SealevelriseDeltathicknessEnum);
 
-        Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //sea-level rise dependent terms (2nd and 5th terms on the RHS)
+        /*call sea-level rise sub cores:*/
+        if(isslr){
+                Sg_eustatic=sealevelrise_core_eustatic(femmodel); //generalized eustatic (Farrel and Clark, Eq 4, 1st, 3rd and 4rd terms on the RHS.
 
-        /*get results out:*/
-        InputUpdateFromVectorx(femmodel,Sg,SealevelriseSEnum,VertexSIdEnum);
+                Sg=sealevelrise_core_noneustatic(femmodel,Sg_eustatic); //sea-level rise dependent terms (2nd and 5th terms on the RHS)
 
-        if(save_results){
-                if(VerboseSolution()) _printf0_("   saving results\n");
-                femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
+                /*get results into elements:*/
+                InputUpdateFromVectorx(femmodel,Sg,SealevelriseSEnum,VertexSIdEnum);
+
+                if(save_results){
+                        if(VerboseSolution()) _printf0_("   saving results\n");
+                        femmodel->parameters->FindParam(&requested_outputs,&numoutputs,SealevelriseRequestedOutputsEnum);
+                        femmodel->RequestedOutputsx(&femmodel->results,requested_outputs,numoutputs);
+                }
+
+                if(solution_type==SealevelriseSolutionEnum)femmodel->RequestedDependentsx();
+
+                /*Free ressources:*/    
+                delete Sg;
+                delete Sg_eustatic;
+                if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
         }
+
+        /*transfer sea-level back to ice caps: */
+        if(iscoupler)TransferSealevel(femmodel,SealevelriseSEnum);
+} 
+/*}}}*/
+void TransferForcing(FemModel* femmodel,int forcingenum){ /*{{{*/
+
+        /*forcing being transferred from models to earth: */
+        IssmDouble** forcings=NULL;
+        IssmDouble*  forcing=NULL; 
+        Vector<IssmDouble>* forcingglobal=NULL; 
+        int*         nvs=NULL;
         
-        /*Free ressources:*/    
-        if(numoutputs){for(int i=0;i<numoutputs;i++){xDelete<char>(requested_outputs[i]);} xDelete<char*>(requested_outputs);}
-        delete Sg;
-        delete Sg_eustatic;
+        /*transition vectors:*/
+        IssmDouble** transitions=NULL;
+        int          ntransitions; 
+        int*         transitions_m=NULL;
+        int*         transitions_n=NULL;
+        int          nv;
+        
+        /*communicators:*/
+        ISSM_MPI_Comm tocomm;
+        ISSM_MPI_Comm* fromcomms=NULL;
+        ISSM_MPI_Status status;
+        int         my_rank;
+        int         modelid,earthid;
+        int         nummodels;
 
+        /*Recover some parameters: */
+        femmodel->parameters->FindParam(&modelid,ModelIdEnum);
+        femmodel->parameters->FindParam(&earthid,EarthIdEnum);
+        femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
+        my_rank=IssmComm::GetRank();
+        
+        /*retrieve the inter communicators that will be used to send data from each ice cap to the earth: */
+        if(modelid==earthid)femmodel->parameters->FindParam(&fromcomms,&nv,IcecapToEarthCommEnum);
+        else femmodel->parameters->FindParam(&tocomm, IcecapToEarthCommEnum);
 
-}
+        /*For each icecap, retrieve the forcing vector that will be sent to the earth model: */
+        if(modelid!=earthid){
+                nv=femmodel->vertices->NumberOfVertices();
+                GetVectorFromInputsx(&forcing,femmodel,forcingenum,VertexSIdEnum);
+        }
+
+        /*Send the forcing to the earth model:{{{*/
+        if(my_rank==0){
+                if(modelid==earthid){
+                        forcings=xNew<IssmDouble*>(nummodels-1);
+                        nvs=xNew<int>(nummodels-1);
+                        for(int i=0;i<earthid;i++){
+                                MPI_Recv(nvs+i, 1, ISSM_MPI_INT, 0,i, fromcomms[i], &status);
+                                forcings[i]=xNew<IssmDouble>(nvs[i]);
+                                MPI_Recv(forcings[i], nvs[i], ISSM_MPI_DOUBLE, 0,i, fromcomms[i], &status);
+                        }
+                        
+                }
+                else{
+                        ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, modelid, tocomm);
+                        ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, tocomm);
+                }
+        }
+        /*}}}*/
+
+        /*On the earth model, consolidate all the forcings into one, and update the elements dataset accordingly: {{{*/
+        if(modelid==earthid){
+                
+                /*Out of all the delta thicknesses, build one delta thickness vector made of all the ice cap contributions. 
+                 *First, build the global delta thickness vector in the earth model: */
+                nv=femmodel->vertices->NumberOfVertices();
+                forcingglobal= new Vector<IssmDouble>(nv);
+
+                /*Retrieve transition vectors, used to plug from each ice cap into the global forcing:*/
+                femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
+
+                if(ntransitions!=earthid)_error_("TransferForcing error message: number of transition vectors is not equal to the number of icecaps!");
+
+                /*Go through all the delta thicknesses coming from each ice cap: */
+                if(my_rank==0){
+                        for(int i=0;i<earthid;i++){
+
+                                IssmDouble* forcingfromcap= forcings[i]; //careful, this only exists on rank 0 of the earth model!
+                                IssmDouble* transition=transitions[i];
+                                int         M=transitions_m[i];
+
+                                /*build index to plug values: */
+                                int*        index=xNew<int>(M); for(int i=0;i<M;i++)index[i]=reCast<int>(transition[i])-1; //matlab indexing!
+
+
+                                /*We are going to plug this vector into the earth model, at the right vertices corresponding to this particular 
+                                 * ice cap: */
+                                forcingglobal->SetValues(M,index,forcingfromcap,INS_VAL);
+                                xDelete<int>(index);
+                        }
+                }
+
+                /*Assemble vector:*/
+                forcingglobal->Assemble();
+                
+                /*Plug into elements:*/
+                InputUpdateFromVectorx(femmodel,forcingglobal,forcingenum,VertexSIdEnum);
+        } 
+        /*}}}*/
+
+        /*Free ressources:{{{*/
+        if(forcings){
+                for(int i=0;i<nummodels-1;i++){
+                        IssmDouble* temp=forcings[i]; xDelete<IssmDouble>(temp);
+                }
+                xDelete<IssmDouble*>(forcings);
+        }
+        if(forcing)xDelete<IssmDouble>(forcing);
+        if(forcingglobal)delete forcingglobal;
+        if(transitions){
+                for(int i=0;i<earthid;i++){
+                        IssmDouble* temp=transitions[i];
+                        xDelete<IssmDouble>(temp);
+                }
+                xDelete<IssmDouble*>(transitions);
+                xDelete<int>(transitions_m);
+                xDelete<int>(transitions_n);
+        }
+        if(nvs)xDelete<int>(nvs);
+        /*}}}*/
+
+} /*}}}*/
+void TransferSealevel(FemModel* femmodel,int forcingenum){ /*{{{*/
+
+        /*forcing being transferred from earth to ice caps: */
+        IssmDouble*  forcing=NULL; 
+        IssmDouble*  forcingglobal=NULL; 
+        
+        /*transition vectors:*/
+        IssmDouble** transitions=NULL;
+        int          ntransitions; 
+        int*         transitions_m=NULL;
+        int*         transitions_n=NULL;
+        int          nv;
+        
+        /*communicators:*/
+        ISSM_MPI_Comm fromcomm;
+        ISSM_MPI_Comm* tocomms=NULL;
+        ISSM_MPI_Status status;
+        int         my_rank;
+        int         modelid,earthid;
+        int         nummodels;
+        int         numcoms;
+
+        /*Recover some parameters: */
+        femmodel->parameters->FindParam(&modelid,ModelIdEnum);
+        femmodel->parameters->FindParam(&earthid,EarthIdEnum);
+        femmodel->parameters->FindParam(&nummodels,NumModelsEnum);
+        my_rank=IssmComm::GetRank();
+        
+        /*retrieve the inter communicators that will be used to send data from earth to ice caps:*/
+        if(modelid==earthid)femmodel->parameters->FindParam(&tocomms,&numcoms,IcecapToEarthCommEnum);
+        else femmodel->parameters->FindParam(&fromcomm, IcecapToEarthCommEnum);
+
+
+        /*Retrieve sea-level on earth model: */
+        if(modelid==earthid){
+                nv=femmodel->vertices->NumberOfVertices();
+                GetVectorFromInputsx(&forcingglobal,femmodel,forcingenum,VertexSIdEnum);
+        }
+
+        /*Send the forcing to the ice caps:{{{*/
+        if(my_rank==0){
+                
+                if(modelid==earthid){
+                        
+                        /*Retrieve transition vectors, used to figure out global forcing contribution to each ice cap's own elements: */
+                        femmodel->parameters->FindParam(&transitions,&ntransitions,&transitions_m,&transitions_n,SealevelriseTransitionsEnum);
+                        
+                        if(ntransitions!=earthid)_error_("TransferSeaLevel error message: number of transition vectors is not equal to the number of icecaps!");
+
+                        for(int i=0;i<earthid;i++){
+                                nv=transitions_m[i];
+                                forcing=xNew<IssmDouble>(nv);
+                                IssmDouble* transition=transitions[i];
+                                for(int j=0;j<nv;j++){
+                                        forcing[j]=forcingglobal[reCast<int>(transition[j])-1];
+                                }
+                                ISSM_MPI_Send(&nv, 1, ISSM_MPI_INT, 0, i, tocomms[i]);
+                                ISSM_MPI_Send(forcing, nv, ISSM_MPI_DOUBLE, 0, i, tocomms[i]);
+                        }
+                }
+                else{
+                        ISSM_MPI_Recv(&nv, 1, ISSM_MPI_INT, 0, modelid, fromcomm, &status);
+                        forcing=xNew<IssmDouble>(nv);
+                        ISSM_MPI_Recv(forcing, nv, ISSM_MPI_DOUBLE, 0, modelid, fromcomm, &status);
+                }
+        }
+        /*}}}*/
+
+        /*On each ice cap, spread the forcing across cpus, and update the elements dataset accordingly: {{{*/
+        if(modelid!=earthid){
+
+                ISSM_MPI_Bcast(&nv,1,ISSM_MPI_INT,0,IssmComm::GetComm());
+                if(my_rank!=0)forcing=xNew<IssmDouble>(nv);
+                ISSM_MPI_Bcast(forcing,nv,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
+
+                /*Plug into elements:*/
+                InputUpdateFromVectorx(femmodel,forcing,forcingenum,VertexSIdEnum);
+        } 
+        /*}}}*/
+
+        /*Free ressources:{{{*/
+        if(forcingglobal)xDelete<IssmDouble>(forcingglobal);
+        if(forcing)xDelete<IssmDouble>(forcing);
+        if(transitions){
+                for(int i=0;i<ntransitions;i++){
+                        IssmDouble* temp=transitions[i];
+                        xDelete<IssmDouble>(temp);
+                }
+                xDelete<IssmDouble*>(transitions);
+                xDelete<int>(transitions_m);
+                xDelete<int>(transitions_n);
+        }
+        /*}}}*/
+
+} /*}}}*/