Changeset 23067


Ignore:
Timestamp:
08/07/18 10:29:35 (7 years ago)
Author:
seroussi
Message:

FIX: fixed coupling for multiple ice cpus

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r22993 r23067  
    8989                tomitgcmcomm=parcom->GetParameterValue();
    9090
     91                int oceangridnxsize,oceangridnysize,ngrids_ocean,nels_ocean;
     92                IssmDouble  oceantime,coupling_time;
     93                IssmDouble *oceangridx;
     94                IssmDouble *oceangridy;
     95                IssmDouble *base_oceangrid;
     96                IssmDouble* x_ice = NULL;
     97                IssmDouble* y_ice = NULL;
     98                IssmDouble* lat_ice = NULL;
     99                IssmDouble* lon_ice = NULL;
     100                IssmDouble* z_ice = NULL;
     101                IssmDouble* icebase= NULL;
     102                int*        index_ice= NULL;
     103                int*        index_ocean = NULL;
     104                int         ngrids_ice=femmodel->vertices->NumberOfVertices();
     105                int         nels_ice=femmodel->elements->NumberOfElements();
     106
     107                /*Recover fixed parameters and store them*/
     108                femmodel->parameters->FindParam(&coupling_time,TimesteppingCouplingTimeEnum);
    91109                if(my_rank==0){
    92                         int oceangridnxsize,oceangridnysize,ngrids_ocean,nels_ocean;
    93                         IssmDouble  oceantime,coupling_time;
    94                         IssmDouble *oceangridx;
    95                         IssmDouble *oceangridy;
    96                         IssmDouble *base_oceangrid;
    97                         IssmDouble* x_ice = NULL;
    98                         IssmDouble* y_ice = NULL;
    99                         IssmDouble* lat_ice = NULL;
    100                         IssmDouble* lon_ice = NULL;
    101                         IssmDouble* z_ice = NULL;
    102                         IssmDouble* icebase= NULL;
    103                         int*        index_ice= NULL;
    104                         int*        index_ocean = NULL;
    105                         int         ngrids_ice=femmodel->vertices->NumberOfVertices();
    106                         int         nels_ice=femmodel->elements->NumberOfElements();
    107 
    108                         /*Recover fixed parameters and store them*/
    109                         femmodel->parameters->FindParam(&coupling_time,TimesteppingCouplingTimeEnum);
    110110                        ISSM_MPI_Send(&coupling_time,1,ISSM_MPI_DOUBLE,0,10001000,tomitgcmcomm);
    111111                        ISSM_MPI_Recv(&oceangridnxsize,1,ISSM_MPI_INT,0,10001003,tomitgcmcomm,&status);
    112112                        ISSM_MPI_Recv(&oceangridnysize,1,ISSM_MPI_INT,0,10001004,tomitgcmcomm,&status);
    113                         femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum);
    114                         femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum);
    115                         ngrids_ocean=oceangridnxsize*oceangridnysize;
     113                }
     114                ngrids_ocean=oceangridnxsize*oceangridnysize;
     115                ISSM_MPI_Bcast(&oceangridnxsize,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     116                ISSM_MPI_Bcast(&oceangridnysize,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     117                ISSM_MPI_Bcast(&ngrids_ocean,1,ISSM_MPI_INT,0,IssmComm::GetComm());
     118                ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     119                femmodel->parameters->SetParam(oceangridnxsize,OceanGridNxEnum);
     120                femmodel->parameters->SetParam(oceangridnysize,OceanGridNyEnum);
     121                if(my_rank==0){
    116122                        oceangridx = xNew<IssmDouble>(ngrids_ocean);
    117123                        ISSM_MPI_Recv(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001005,tomitgcmcomm,&status);
    118124                        oceangridy = xNew<IssmDouble>(ngrids_ocean);
    119125                        ISSM_MPI_Recv(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status);
    120                         femmodel->parameters->SetParam(oceangridx,ngrids_ocean,OceanGridXEnum);
    121                         femmodel->parameters->SetParam(oceangridy,ngrids_ocean,OceanGridYEnum);
    122126
    123127                        /*Exchange varying parameters for the initialization*/
    124128                        ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
    125129                        ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
    126 
    127                         /*Interpolate ice base onto ocean grid*/
    128                         femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice);
    129                         BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean);
    130                         femmodel->vertices->LatLonList(&lat_ice,&lon_ice);
    131                         GetVectorFromInputsx(&icebase,femmodel,BaseEnum,VertexSIdEnum);
    132                         Options* options = new Options();
    133                         GenericOption<double> *odouble = new GenericOption<double>();
    134                         const char* name = "default";
    135                         odouble->name =xNew<char>(strlen(name)+1);
    136                         memcpy(odouble->name,name,(strlen(name)+1)*sizeof(char));
    137                         odouble->value=+9999.;
    138                         odouble->numel=1;
    139                         odouble->ndims=1;
    140                         odouble->size=NULL;
    141                         options->AddOption(odouble);
    142                         InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
    143                                                         icebase,ngrids_ice,1,
    144                                                         oceangridx,oceangridy,ngrids_ocean,options);
    145                         delete options;
    146 
     130                }
     131                if(my_rank!=0){
     132                        oceangridx=xNew<IssmDouble>(ngrids_ocean);
     133                        oceangridy=xNew<IssmDouble>(ngrids_ocean);
     134                }
     135                ISSM_MPI_Bcast(oceangridx,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     136                ISSM_MPI_Bcast(oceangridy,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     137                femmodel->parameters->SetParam(oceangridx,ngrids_ocean,OceanGridXEnum);
     138                femmodel->parameters->SetParam(oceangridy,ngrids_ocean,OceanGridYEnum);
     139
     140                /*Interpolate ice base onto ocean grid*/
     141                femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice);
     142                BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean);
     143                femmodel->vertices->LatLonList(&lat_ice,&lon_ice);
     144                GetVectorFromInputsx(&icebase,femmodel,BaseEnum,VertexSIdEnum);
     145                Options* options = new Options();
     146                GenericOption<double> *odouble = new GenericOption<double>();
     147                const char* name = "default";
     148                odouble->name =xNew<char>(strlen(name)+1);
     149                memcpy(odouble->name,name,(strlen(name)+1)*sizeof(char));
     150                odouble->value=+9999.;
     151                odouble->numel=1;
     152                odouble->ndims=1;
     153                odouble->size=NULL;
     154                options->AddOption(odouble);
     155                InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
     156                                                icebase,ngrids_ice,1,
     157                                                oceangridx,oceangridy,ngrids_ocean,options);
     158                delete options;
     159
     160                if(my_rank==0){
    147161                        ISSM_MPI_Send(base_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
    148 
    149                         /*Delete*/
    150                         xDelete<int>(index_ice);
    151                         xDelete<int>(index_ocean);
    152                         xDelete<IssmDouble>(lat_ice);
    153                         xDelete<IssmDouble>(lon_ice);
    154                         xDelete<IssmDouble>(x_ice);
    155                         xDelete<IssmDouble>(y_ice);
    156                         xDelete<IssmDouble>(z_ice);
    157                         xDelete<IssmDouble>(icebase);
    158                         xDelete<IssmDouble>(base_oceangrid);
    159                         xDelete<IssmDouble>(oceangridx);
    160                         xDelete<IssmDouble>(oceangridy);
    161                 }
    162                 #else
    163                 _error_("not supported");
    164                 #endif
     162                }
     163
     164                /*Delete*/
     165                xDelete<int>(index_ice);
     166                xDelete<int>(index_ocean);
     167                xDelete<IssmDouble>(lat_ice);
     168                xDelete<IssmDouble>(lon_ice);
     169                xDelete<IssmDouble>(x_ice);
     170                xDelete<IssmDouble>(y_ice);
     171                xDelete<IssmDouble>(z_ice);
     172                xDelete<IssmDouble>(icebase);
     173                xDelete<IssmDouble>(base_oceangrid);
     174                xDelete<IssmDouble>(oceangridx);
     175                xDelete<IssmDouble>(oceangridy);
    165176        }
     177        #else
     178        _error_("not supported");
     179        #endif
    166180        /*}}}*/
    167181
    … …  
    220234                        if(!parcom)_error_("TransferForcing error message: could not find ToMITgcmCommEnum communicator");
    221235                        tomitgcmcomm=parcom->GetParameterValue();
     236                        int ngrids_ocean, nels_ocean;
     237                        IssmDouble oceantime;
     238                        IssmDouble rho_ice;
     239                        IssmDouble *oceanmelt;
     240                        IssmDouble *base_oceangrid;
     241                        IssmDouble *oceangridx;
     242                        IssmDouble *oceangridy;
     243                        IssmDouble* x_ice = NULL;
     244                        IssmDouble* y_ice = NULL;
     245                        IssmDouble* lat_ice = NULL;
     246                        IssmDouble* lon_ice = NULL;
     247                        IssmDouble* z_ice = NULL;
     248                        IssmDouble* icebase= NULL;
     249                        IssmDouble* melt_mesh = NULL;
     250                        int*        index_ice= NULL;
     251                        int*        index_ocean = NULL;
     252                        int         ngrids_ice=femmodel->vertices->NumberOfVertices();
     253                        int         nels_ice=femmodel->elements->NumberOfElements();
     254
     255                        /*Recover mesh and inputs needed*/
     256                        femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
     257                        femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice);
     258                        femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum);
     259                        femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum);
     260                        BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean);
     261
     262                        femmodel->vertices->LatLonList(&lat_ice,&lon_ice);
     263
     264                        /*Interpolate ice base onto ocean grid*/
     265                        GetVectorFromInputsx(&icebase,femmodel,BaseEnum,VertexSIdEnum);
     266                        InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
     267                                                icebase,ngrids_ice,1,
     268                                                oceangridx,oceangridy,ngrids_ocean,NULL);
     269
     270                        /*Send and receive data*/
    222271                        if(my_rank==0){
    223                                 int ngrids_ocean, nels_ocean;
    224                                 IssmDouble oceantime;
    225                                 IssmDouble rho_ice;
    226                                 IssmDouble *oceanmelt;
    227                                 IssmDouble *base_oceangrid;
    228                                 IssmDouble *oceangridx;
    229                                 IssmDouble *oceangridy;
    230                                 IssmDouble* x_ice = NULL;
    231                                 IssmDouble* y_ice = NULL;
    232                                 IssmDouble* lat_ice = NULL;
    233                                 IssmDouble* lon_ice = NULL;
    234                                 IssmDouble* z_ice = NULL;
    235                                 IssmDouble* icebase= NULL;
    236                                 IssmDouble* melt_mesh = NULL;
    237                                 int*        index_ice= NULL;
    238                                 int*        index_ocean = NULL;
    239                                 int         ngrids_ice=femmodel->vertices->NumberOfVertices();
    240                                 int         nels_ice=femmodel->elements->NumberOfElements();
    241 
    242                                 /*Recover mesh and inputs needed*/
    243                                 femmodel->parameters->FindParam(&rho_ice,MaterialsRhoIceEnum);
    244                                 femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice);
    245                                 femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum);
    246                                 femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum);
    247                                 BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean);
    248 
    249                                 femmodel->vertices->LatLonList(&lat_ice,&lon_ice);
    250 
    251                                 /*Interpolate ice base onto ocean grid*/
    252                                 GetVectorFromInputsx(&icebase,femmodel,BaseEnum,VertexSIdEnum);
    253                                 InterpFromMeshToMesh2dx(&base_oceangrid,index_ice,lon_ice,lat_ice,ngrids_ice,nels_ice,
    254                                                         icebase,ngrids_ice,1,
    255                                                         oceangridx,oceangridy,ngrids_ocean,NULL);
    256 
    257                                 /*Send and receive data*/
    258272                                ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
    259                                 oceanmelt = xNew<IssmDouble>(ngrids_ocean);
    260273                                ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
    261274                                if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge");
     275                                oceanmelt = xNew<IssmDouble>(ngrids_ocean);
    262276                                ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status);
    263                                 //for(int i=0;i<ngrids_ice;i++) base_oceangrid[i]=0;
    264277                                ISSM_MPI_Send(base_oceangrid,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
    265 
    266                                 /*Interp melt onto ice grid*/
    267                                 InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean,
    268                                                         oceanmelt,ngrids_ocean,1,
    269                                                         lon_ice,lat_ice,ngrids_ice,NULL);
    270 
    271                                 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
    272                                 InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,VertexSIdEnum);
    273 
    274                                 /*Delete*/
    275                                 xDelete<int>(index_ice);
    276                                 xDelete<int>(index_ocean);
    277                                 xDelete<IssmDouble>(lat_ice);
    278                                 xDelete<IssmDouble>(lon_ice);
    279                                 xDelete<IssmDouble>(x_ice);
    280                                 xDelete<IssmDouble>(y_ice);
    281                                 xDelete<IssmDouble>(z_ice);
    282                                 xDelete<IssmDouble>(melt_mesh);
    283                                 xDelete<IssmDouble>(oceangridx);
    284                                 xDelete<IssmDouble>(oceangridy);
    285                                 xDelete<IssmDouble>(oceanmelt);
    286                                 xDelete<IssmDouble>(icebase);
    287                                 xDelete<IssmDouble>(base_oceangrid);
    288278                        }
    289                         #else
    290                         _error_("not supported");
    291                         #endif
    292                 }
     279                        ISSM_MPI_Bcast(&oceantime,1,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     280                        if(my_rank!=0) oceanmelt=xNew<IssmDouble>(ngrids_ocean);
     281                        ISSM_MPI_Bcast(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,IssmComm::GetComm());
     282
     283                        /*Interp melt onto ice grid*/
     284                        InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean,
     285                                                oceanmelt,ngrids_ocean,1,
     286                                                lon_ice,lat_ice,ngrids_ice,NULL);
     287
     288                        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
     289                        InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,VertexSIdEnum);
     290
     291                        /*Delete*/
     292                        xDelete<int>(index_ice);
     293                        xDelete<int>(index_ocean);
     294                        xDelete<IssmDouble>(lat_ice);
     295                        xDelete<IssmDouble>(lon_ice);
     296                        xDelete<IssmDouble>(x_ice);
     297                        xDelete<IssmDouble>(y_ice);
     298                        xDelete<IssmDouble>(z_ice);
     299                        xDelete<IssmDouble>(melt_mesh);
     300                        xDelete<IssmDouble>(oceangridx);
     301                        xDelete<IssmDouble>(oceangridy);
     302                        xDelete<IssmDouble>(oceanmelt);
     303                        xDelete<IssmDouble>(icebase);
     304                        xDelete<IssmDouble>(base_oceangrid);
     305                }
     306                #else
     307                _error_("not supported");
     308                #endif
    293309                /*}}}*/
    294310
Note: See TracChangeset for help on using the changeset viewer.