Changeset 23717


Ignore:
Timestamp:
02/12/19 16:16:35 (6 years ago)
Author:
seroussi
Message:

CHG: finished cleaning up ocean coupling

Location:
issm/trunk-jpl/src/c
Files:
3 edited

Legend:

Unmodified
Added
Removed

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

    r23715 r23717  
    8080
    8181        #if defined(_HAVE_OCEAN_ )
    82         if(isoceancoupling) OceanExchangeDatax(femmodel);
     82        if(isoceancoupling) OceanExchangeDatax(femmodel,true);
    8383        #endif
    8484
    … …  
    126126                femmodel->parameters->SetParam(save_results,SaveResultsEnum);
    127127
    128                 if(isoceancoupling){ /*{{{*/
    129 
    130                         #ifndef _HAVE_AD_
    131                         if(VerboseSolution()) _printf0_("   ocean coupling: exchanging information\n");
    132                         int my_rank;
    133                         ISSM_MPI_Comm tomitgcmcomm;
    134                         ISSM_MPI_Status status;
    135 
    136                         my_rank=IssmComm::GetRank();
    137                         GenericParam<ISSM_MPI_Comm>* parcom = dynamic_cast<GenericParam<ISSM_MPI_Comm>*>(femmodel->parameters->FindParamObject(ToMITgcmCommEnum));
    138                         if(!parcom)_error_("TransferForcing error message: could not find ToMITgcmCommEnum communicator");
    139                         tomitgcmcomm=parcom->GetParameterValue();
    140                         int ngrids_ocean, nels_ocean;
    141                         IssmDouble oceantime;
    142                         IssmDouble rho_ice;
    143                         IssmDouble *oceanmelt         = NULL;
    144                         IssmDouble *icebase_oceangrid = NULL;
    145                         IssmDouble *icemask_oceangrid = NULL;
    146                         IssmDouble *oceangridx        = NULL;
    147                         IssmDouble *oceangridy        = NULL;
    148                         IssmDouble *x_ice             = NULL;
    149                         IssmDouble *y_ice             = NULL;
    150                         IssmDouble *lat_ice           = NULL;
    151                         IssmDouble *lon_ice           = NULL;
    152                         IssmDouble *icebase           = NULL;
    153                         IssmDouble *icemask           = NULL;
    154                         IssmDouble *melt_mesh         = NULL;
    155                         int        *index_ice         = NULL;
    156                         int        *index_ocean       = NULL;
    157                         int         ngrids_ice=femmodel->vertices->NumberOfVertices();
    158                         int         nels_ice=femmodel->elements->NumberOfElements();
    159 
    160                         /*Recover mesh and inputs needed*/
    161                         femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
    162                         femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&index_ice);
    163                         femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum);
    164                         femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum);
    165                         BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean);
    166 
    167                         femmodel->vertices->LatLonList(&lat_ice,&lon_ice);
    168 
    169                         /*Interpolate ice base and mask onto ocean grid*/
    170                         GetVectorFromInputsx(&icebase,femmodel,BaseEnum,VertexSIdEnum);
    171                         Options* options = new Options();
    172                         GenericOption<double> *odouble = new GenericOption<double>();
    173                         const char* name = "default";
    174                         odouble->name =xNew<char>(strlen(name)+1);
    175                         memcpy(odouble->name,name,(strlen(name)+1)*sizeof(char));
    176                         odouble->value=+9999.;
    177                         odouble->size[0]=1;
    178                         odouble->size[1]=1;
    179                         options->AddOption(odouble);
    180                         InterpFromMeshToMesh2dx(&icebase_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
    181                                                 icebase,ngrids_ice,1,oceangridx,oceangridy,ngrids_ocean,options);
    182                         delete options;
    183                         xDelete<IssmDouble>(icebase);
    184 
    185                         GetVectorFromInputsx(&icemask,femmodel,MaskIceLevelsetEnum,VertexSIdEnum);
    186                         Options* options2 = new Options();
    187                         GenericOption<double> *odouble2 = new GenericOption<double>();
    188                         const char* name2 = "default";
    189                         odouble2->name =xNew<char>(strlen(name2)+1);
    190                         memcpy(odouble2->name,name2,(strlen(name2)+1)*sizeof(char));
    191                         odouble2->value=+1.;
    192                         odouble2->size[0]=1;
    193                         odouble2->size[1]=1;
    194                         options2->AddOption(odouble2);
    195                         InterpFromMeshToMesh2dx(&icemask_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
    196                                                 icemask,ngrids_ice,1,oceangridx,oceangridy,ngrids_ocean,options2);
    197                         delete options2;
    198                         xDelete<IssmDouble>(icemask);
    199 
    200                         /*Put +9999 for places where there is no ice!*/
    201                         for(int i=0;i<ngrids_ocean;i++) if(icemask_oceangrid[i]>0.) icebase_oceangrid[i]=+9999.;
    202                         xDelete<IssmDouble>(icemask_oceangrid);
    203 
    204                         /*Send and receive data*/
    205                         if(my_rank==0){
    206                                 ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
    207                                 ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
    208                                 if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge");
    209                                 oceanmelt = xNew<IssmDouble>(ngrids_ocean);
    210                                 ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status);
    211                                 ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
    212                         }
    213                         ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    214                         if(my_rank!=0) oceanmelt=xNew<IssmDouble>(ngrids_ocean);
    215                         ISSM_MPI_Bcast(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    216 
    217                         /*Interp melt onto ice grid*/
    218                         InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean,
    219                                                 oceanmelt,ngrids_ocean,1,
    220                                                 lon_ice,lat_ice,ngrids_ice,NULL);
    221 
    222                         for(int i=0;i<ngrids_ice;i++) melt_mesh[i]=-melt_mesh[i]/rho_ice; //heat flux provided by ocean is in kg/m^2/s
    223                         InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,VertexSIdEnum);
    224 
    225                         /*Delete*/
    226                         xDelete<int>(index_ice);
    227                         xDelete<int>(index_ocean);
    228                         xDelete<IssmDouble>(lat_ice);
    229                         xDelete<IssmDouble>(lon_ice);
    230                         xDelete<IssmDouble>(x_ice);
    231                         xDelete<IssmDouble>(y_ice);
    232                         xDelete<IssmDouble>(melt_mesh);
    233                         xDelete<IssmDouble>(oceangridx);
    234                         xDelete<IssmDouble>(oceangridy);
    235                         xDelete<IssmDouble>(oceanmelt);
    236                         xDelete<IssmDouble>(icebase_oceangrid);
    237 
    238                 #else
    239                 _error_("not supported");
    240                 #endif
    241                 }
    242                 /*}}}*/
     128        #if defined(_HAVE_OCEAN_ )
     129        if(isoceancoupling) OceanExchangeDatax(femmodel,false);
     130        #endif
    243131
    244132                if(isthermal && domaintype==Domain3DEnum){

  • issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.cpp

    r23714 r23717  
    99#include "./OceanExchangeDatax.h"
    1010
    11 void OceanExchangeDatax(FemModel* femmodel){
     11void OceanExchangeDatax(FemModel* femmodel, bool init_stage){
    1212
    1313        #ifndef _HAVE_AD_
    14         if(VerboseSolution()) _printf0_("   ocean coupling: initialization \n");
     14        if(VerboseSolution()) _printf0_("   ocean coupling: exchanging information\n");
    1515        int my_rank;
    1616        ISSM_MPI_Comm tomitgcmcomm;
    … …  
    2323
    2424        int oceangridnxsize,oceangridnysize,ngrids_ocean,nels_ocean;
    25         IssmDouble  oceantime,coupling_time,time;
     25        IssmDouble  oceantime,coupling_time,time,yts;
     26        IssmDouble rho_ice;
     27        IssmDouble *oceanmelt         = NULL;
    2628        IssmDouble *oceangridx;
    2729        IssmDouble *oceangridy;
    … …  
    3436        IssmDouble* icebase           = NULL;
    3537        IssmDouble* icemask           = NULL;
     38        IssmDouble* melt_mesh         = NULL;
    3639        int*        index_ice         = NULL;
    3740        int*        index_ocean       = NULL;
    … …  
    4144        /*Recover fixed parameters and store them*/
    4245        femmodel->parameters->FindParam(&coupling_time,TimesteppingCouplingTimeEnum);
    43         if(my_rank==0){
    44                 ISSM_MPI_Send(&coupling_time,1,ISSM_MPI_DOUBLE,0,10001000,tomitgcmcomm);
    45                 ISSM_MPI_Recv(&oceangridnxsize,1,ISSM_MPI_INT,0,10001003,tomitgcmcomm,&status);
    46                 ISSM_MPI_Recv(&oceangridnysize,1,ISSM_MPI_INT,0,10001004,tomitgcmcomm,&status);
     46
     47        /*Exchange or recover mesh and inputs needed*/
     48        if(init_stage==true){
     49                if(my_rank==0){
     50                        ISSM_MPI_Send(&coupling_time,1,ISSM_MPI_DOUBLE,0,10001000,tomitgcmcomm);
     51                        ISSM_MPI_Recv(&oceangridnxsize,1,ISSM_MPI_INT,0,10001003,tomitgcmcomm,&status);
     52                        ISSM_MPI_Recv(&oceangridnysize,1,ISSM_MPI_INT,0,10001004,tomitgcmcomm,&status);
     53                }
     54                ngrids_ocean=oceangridnxsize*oceangridnysize;
     55                ISSM_MPI_Bcast(&oceangridnxsize,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     56                ISSM_MPI_Bcast(&oceangridnysize,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     57                ISSM_MPI_Bcast(&ngrids_ocean,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     58                ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     59                femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum);
     60                femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum);
     61                if(my_rank==0){
     62                        oceangridx = xNew<IssmDouble>(ngrids_ocean);
     63                        ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status);
     64                        oceangridy = xNew<IssmDouble>(ngrids_ocean);
     65                        ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status);
     66
     67                        /*Exchange varying parameters for the initialization*/
     68                        ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
     69                        ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
     70                }
     71                if(my_rank!=0){
     72                        oceangridx=xNew<IssmDouble>(ngrids_ocean);
     73                        oceangridy=xNew<IssmDouble>(ngrids_ocean);
     74                }
     75                ISSM_MPI_Bcast(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     76                ISSM_MPI_Bcast(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     77                femmodel->parameters->SetParam(oceangridx,ngrids_ocean,OceanGridXEnum);
     78                femmodel->parameters->SetParam(oceangridy,ngrids_ocean,OceanGridYEnum);
    4779        }
    48         ngrids_ocean=oceangridnxsize*oceangridnysize;
    49         ISSM_MPI_Bcast(&oceangridnxsize,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    50         ISSM_MPI_Bcast(&oceangridnysize,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    51         ISSM_MPI_Bcast(&ngrids_ocean,1,ISSM_MPI_INT,0,IssmComm::GetComm());
    52         ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    53         femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum);
    54         femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum);
    55         if(my_rank==0){
    56                 oceangridx = xNew<IssmDouble>(ngrids_ocean);
    57                 ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status);
    58                 oceangridy = xNew<IssmDouble>(ngrids_ocean);
    59                 ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status);
    60 
    61                 /*Exchange varying parameters for the initialization*/
    62                 ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
    63                 ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
     80        else{
     81                femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum);
     82                femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum);
    6483        }
    65         if(my_rank!=0){
    66                 oceangridx=xNew<IssmDouble>(ngrids_ocean);
    67                 oceangridy=xNew<IssmDouble>(ngrids_ocean);
    68         }
    69         ISSM_MPI_Bcast(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    70         ISSM_MPI_Bcast(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
    71         femmodel->parameters->SetParam(oceangridx,ngrids_ocean,OceanGridXEnum);
    72         femmodel->parameters->SetParam(oceangridy,ngrids_ocean,OceanGridYEnum);
    7384
    7485        /*Interpolate ice base and mask onto ocean grid*/
    … …  
    8798        options->AddOption(odouble);
    8899        InterpFromMeshToMesh2dx(&icebase_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
    89                                         icebase,ngrids_ice,1,
    90                                         oceangridx,oceangridy,ngrids_ocean,options);
     100                                        icebase,ngrids_ice,1,oceangridx,oceangridy,ngrids_ocean,options);
    91101        delete options;
    92102        xDelete<IssmDouble>(icebase);
    … …  
    111121        xDelete<IssmDouble>(icemask_oceangrid);
    112122
    113         if(my_rank==0){
    114                 ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
     123        if(init_stage==true){ //just send icebase
     124                if(my_rank==0){
     125                        ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
     126                }
     127        }
     128        else{ //send and receive exchanged data
     129                femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
     130                femmodel->parameters->FindParam(&yts,ConstantsYtsEnum);
     131                if(my_rank==0){
     132                        ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
     133                        ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
     134                        if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge");
     135                        oceanmelt = xNew<IssmDouble>(ngrids_ocean);
     136                        ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status);
     137                        ISSM_MPI_Send(icebase_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
     138                }
     139                ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     140                if(my_rank!=0) oceanmelt=xNew<IssmDouble>(ngrids_ocean);
     141                ISSM_MPI_Bcast(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     142
     143                /*Interp melt onto ice grid*/
     144                InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean,
     145                                        oceanmelt,ngrids_ocean,1,
     146                                        lon_ice,lat_ice,ngrids_ice,NULL);
     147
     148                for(int i=0;i<ngrids_ice;i++) melt_mesh[i]=-melt_mesh[i]/rho_ice; //heat flux provided by ocean is in kg/m^2/s
     149                InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,VertexSIdEnum);
    115150        }
    116151
    … …  
    125160        xDelete<IssmDouble>(oceangridx);
    126161        xDelete<IssmDouble>(oceangridy);
     162        xDelete<IssmDouble>(melt_mesh);
     163        xDelete<IssmDouble>(oceanmelt);
    127164        #else
    128165        _error_("not supported");

  • issm/trunk-jpl/src/c/modules/OceanExchangeDatax/OceanExchangeDatax.h

    r23714 r23717  
    99
    1010/* local prototypes: */
    11 void OceanExchangeDatax(FemModel* femmodel);
     11void OceanExchangeDatax(FemModel* femmodel, bool init_stage);
    1212#endif  /* _OCEANEXCHANGEDATAX_H */
Note: See TracChangeset for help on using the changeset viewer.