Changeset 22778

Timestamp:

05/15/18 13:42:52 (7 years ago)

Author:

seroussi

Message:

NEW: starting to interpolate data between ice and ocean grids

Location:

issm/trunk-jpl/src/c

Files:

• classes/Vertex.cpp (modified) (2 diffs)
• classes/Vertex.h (modified) (1 diff)
• classes/Vertices.cpp (modified) (1 diff)
• classes/Vertices.h (modified) (1 diff)
• cores/transient_core.cpp (modified) (5 diffs)
• modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp (modified) (4 diffs)

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

@@ -31 +31 @@
         this->z            = iomodel->Data("md.mesh.z")[i];
         this->domaintype     = iomodel->domaintype;
+
+        this->latitute     = iomodel->Data("md.mesh.lat")[i];
+        this->longitude    = iomodel->Data("md.mesh.long")[i];
 
         switch(iomodel->domaintype){
@@ -195 +198 @@
 /*}}}*/
 int        Vertex::Sid(void){ return sid; }/*{{{*/
-/*}}}*/
-void       Vertex::ToXYZ(Matrix<IssmDouble>* matrix){/*{{{*/
-
-        IssmDouble xyz[3];
-        int        indices[3];
-
-        if (this->clone==true) return;
-
-        xyz[0]=x;
-        xyz[1]=y; 
-        xyz[2]=z;
-        indices[0]=0;
-        indices[1]=1; 
-        indices[2]=2;
-
-        matrix->SetValues(1,&sid,3,&indices[0],&xyz[0],INS_VAL);
-}
 /*}}}*/
 void       Vertex::UpdateClonePids(int* alltruepids){/*{{{*/

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

@@ -62 +62 @@
                 void       ShowTruePids(int* borderpids);
                 int        Sid(void); 
-                void       ToXYZ(Matrix<IssmDouble>* matrix);
                 void       UpdateClonePids(int* allborderpids);
                 void       UpdatePosition(Vector<IssmDouble>* vx,Vector<IssmDouble>* vy,Vector<IssmDouble>* vz,Parameters* parameters,IssmDouble* thickness,IssmDouble* bed);

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

@@ -180 +180 @@
 }
 /*}}}*/
-IssmDouble* Vertices::ToXYZ(void){/*{{{*/
-
-        /*intermediary: */
-        int i;
-        int my_rank;
-        int num_vertices;
+void Vertices::LatLonList(IssmDouble** plat,IssmDouble** plon){/*{{{*/
 
         /*output: */
-        Matrix<IssmDouble>* xyz = NULL;
         IssmDouble* xyz_serial=NULL;
 
         /*recover my_rank:*/
-        my_rank=IssmComm::GetRank();
+        int my_rank=IssmComm::GetRank();
 
         /*First, figure out number of vertices: */
-        num_vertices=this->NumberOfVertices();
-
-        /*Now, allocate matrix to hold all the vertices x,y and z values: */
-        xyz= new Matrix<IssmDouble>(num_vertices,3);
+        int num_vertices=this->NumberOfVertices();
+
+        /*Now, allocate vectors*/
+        Vector<IssmDouble>* lat = new Vector<IssmDouble>(num_vertices);
+        Vector<IssmDouble>* lon = new Vector<IssmDouble>(num_vertices);
 
         /*Go through vertices, and for each vertex, object, report it cpu: */
-        for(i=0;i<this->Size();i++){
-
-                /*let vertex fill matrix: */
-                Vertex* vertex=xDynamicCast<Vertex*>(this->GetObjectByOffset(i));
-                vertex->ToXYZ(xyz);
+        for(int i=0;i<this->Size();i++){
+                Vertex* vertex=xDynamicCast<Vertex*>(this->GetObjectByOffset(i));
+                lat->SetValue(vertex->sid,vertex->GetLatitude() ,INS_VAL);
+                lon->SetValue(vertex->sid,vertex->GetLongitude(),INS_VAL);
         }
 
         /*Assemble:*/
-        xyz->Assemble();
+        lat->Assemble();
+        lon->Assemble();
 
         /*gather on cpu 0: */
-        xyz_serial=xyz->ToSerial();
-
+        IssmDouble* lat_serial=lat->ToMPISerial();
+        IssmDouble* lon_serial=lon->ToMPISerial();
+
+        printf("lat_serial %g %g %g \n",lat_serial[0],lat_serial[2], lat_serial[799]);
         /*free ressources: */
-        delete xyz;
-        if(my_rank!=0)delete xyz_serial;
-
-        /*return matrix: */
-        return xyz_serial;
-}
-/*}}}*/
+        *plat = lat_serial;
+        *plon = lon_serial;
+        //delete lat;
+        //delete lon;
+}
+/*}}}*/

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

@@ -25 +25 @@
                 int   NumberOfVertices(void);
                 void  Ranks(int* ranks);
-                IssmDouble* ToXYZ(void);
+                void LatLonList(IssmDouble** lat,IssmDouble** lon);
 };
 

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

@@ -93 +93 @@
                 tomitgcmcomm=parcom->GetParameterValue();
                 if(my_rank==0){
+                        /*Recover fixed parameters and store them*/
                         ISSM_MPI_Send(&coupling_time,1,ISSM_MPI_DOUBLE,0,10001000,tomitgcmcomm);
                         ISSM_MPI_Recv(&oceangridnxsize,1,ISSM_MPI_INT,0,10001003,tomitgcmcomm,&status);
@@ -102 +103 @@
                         oceangridy = xNew<IssmDouble>(oceangridnxsize*oceangridnysize);
                         ISSM_MPI_Recv(oceangridy,oceangridnxsize*oceangridnysize,ISSM_MPI_DOUBLE,0,10001006,tomitgcmcomm,&status);
-
+                        femmodel->parameters->SetParam(oceangridx,oceangridnxsize*oceangridnysize,OceanGridXEnum);
+                        femmodel->parameters->SetParam(oceangridy,oceangridnxsize*oceangridnysize,OceanGridYEnum);
+
+                        /*Exchange varying parameters for the initialization*/
                         ISSM_MPI_Send(&time,1,ISSM_MPI_DOUBLE,0,10001001,tomitgcmcomm);
                         ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
@@ -162 +166 @@
                 femmodel->parameters->SetParam(save_results,SaveResultsEnum);
 
-                if(isoceancoupling){ {{{
+                if(isoceancoupling){ /*{{{*/
                         if(VerboseSolution()) _printf0_("   ocean coupling: exchanging information\n");
                         int my_rank;
-                        int oceangridnxsize,oceangridnysize;
+                        int oceangridnxsize,oceangridnysize,ngrids_ocean;
                         IssmDouble *oceanmelt;
                         IssmDouble *icebase;
+                        IssmDouble *oceangridx;
+                        IssmDouble *oceangridy;
                         IssmDouble oceantime;
                         ISSM_MPI_Comm tomitgcmcomm;
@@ -180 +186 @@
                                 ISSM_MPI_Recv(&oceantime,1,ISSM_MPI_DOUBLE,0,10001002,tomitgcmcomm,&status);
                                 if((oceantime - time > 0.1*yts) & (oceantime - time < -0.1*yts)) _error_("Ocean and ice time are starting to diverge");
+
+                                /*Recover inputs needed*/
                                 femmodel->parameters->FindParam(&oceangridnxsize,OceanGridNxEnum);
                                 femmodel->parameters->FindParam(&oceangridnysize,OceanGridNyEnum);
-                                oceanmelt = xNew<IssmDouble>(oceangridnxsize*oceangridnysize);
-                                ISSM_MPI_Recv(oceanmelt,oceangridnxsize*oceangridnysize,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status);
-                                icebase = xNew<IssmDouble>(oceangridnxsize*oceangridnysize);
+                                femmodel->parameters->FindParam(&oceangridx,&ngrids_ocean,OceanGridXEnum);
+                                femmodel->parameters->FindParam(&oceangridy,&ngrids_ocean,OceanGridYEnum);
+
+                                /*Exchange information*/
+                                oceanmelt = xNew<IssmDouble>(ngrids_ocean);
+                                ISSM_MPI_Recv(oceanmelt,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001007,tomitgcmcomm,&status);
+                                icebase = xNew<IssmDouble>(ngrids_ocean);
                                 for(int i=0;i<oceangridnxsize;i++){
                                         for(int j=0;j<oceangridnysize;j++){
@@ -190 +202 @@
                                         }
                                 }
-                                ISSM_MPI_Send(icebase,oceangridnxsize*oceangridnysize,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
+                                ISSM_MPI_Send(icebase,ngrids_ocean,ISSM_MPI_DOUBLE,0,10001008,tomitgcmcomm);
+
+                                /*Interpolate melt onto mesh*/
+                                IssmDouble* x_ice = NULL;
+                                IssmDouble* y_ice = NULL;
+                                IssmDouble* lat_ice = NULL;
+                                IssmDouble* lon_ice = NULL;
+                                IssmDouble* z_ice = NULL;
+                                IssmDouble* melt_mesh = NULL;
+                                int*        index_ice= NULL;
+                                int*        index_ocean = NULL;
+                                int         nels_ocean;
+                                int         ngrids_ice=femmodel->vertices->NumberOfVertices();
+                                lat_ice= xNew<IssmDouble>(ngrids_ice);
+                                lon_ice= xNew<IssmDouble>(ngrids_ice);
+                                melt_mesh= xNew<IssmDouble>(ngrids_ice);
+                                femmodel->GetMesh(femmodel->vertices,femmodel->elements,&x_ice,&y_ice,&z_ice,&index_ice);
+                                BamgTriangulatex(&index_ocean,&nels_ocean,oceangridx,oceangridy,ngrids_ocean);
+                                femmodel->vertices->LatLonList(&lat_ice,&lon_ice);
+                                InterpFromMeshToMesh2dx(&melt_mesh,index_ocean,oceangridx,oceangridy,ngrids_ocean,nels_ocean,
+                                                        oceanmelt,ngrids_ocean,1,
+                                                        lat_ice,lon_ice,ngrids_ice,NULL);           /*FIXME: may be wrong: y=lat x=lon*/
+                                InputUpdateFromVectorx(femmodel,melt_mesh,BasalforcingsFloatingiceMeltingRateEnum,VertexSIdEnum);
+                                xDelete<IssmDouble>(lat_ice);
+                                xDelete<IssmDouble>(lon_ice);
+                                xDelete<IssmDouble>(x_ice);
+                                xDelete<IssmDouble>(y_ice);
+                                xDelete<IssmDouble>(z_ice);
+                                xDelete<IssmDouble>(melt_mesh);
+                                xDelete<int>(index_ice);
+                                xDelete<IssmDouble>(oceangridx);
+                                xDelete<IssmDouble>(oceangridy);
                                 xDelete<IssmDouble>(oceanmelt);
                                 xDelete<IssmDouble>(icebase);
                         }
                 }
-                }}}
+                /*}}}*/
 
                 if(isthermal && domaintype==Domain3DEnum){ 

issm/trunk-jpl/src/c/modules/ModelProcessorx/CreateElementsVerticesAndMaterials.cpp

@@ -16 +16 @@
         bool dakota_analysis;
         bool adolc_analysis;
+        bool isoceancoupling;
         int nnat,dummy;
         int* nature=NULL;
@@ -24 +25 @@
         iomodel->FindConstant(&materials_type,"md.materials.type");
         iomodel->FindConstant(&adolc_analysis,"md.autodiff.isautodiff");
+        iomodel->FindConstant(&isoceancoupling,"md.transient.isoceancoupling");
 
         /*Did we already create the elements? : */
@@ -232 +234 @@
         if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->FetchData(3,"md.mesh.lat","md.mesh.long","md.mesh.r");
         else iomodel->FetchDataToInput(elements,"md.mesh.scale_factor",MeshScaleFactorEnum,1.);
+        if (isoceancoupling) iomodel->FetchData(2,"md.mesh.lat","md.mesh.long");
         
         CreateNumberNodeToElementConnectivity(iomodel,solution_type);
@@ -242 +245 @@
         iomodel->DeleteData(6,"md.mesh.x","md.mesh.y","md.mesh.z","md.geometry.base","md.geometry.thickness","md.mask.ice_levelset");
         if (iomodel->domaintype == Domain3DsurfaceEnum) iomodel->DeleteData(3,"md.mesh.lat","md.mesh.long","md.mesh.r");
+        if (isoceancoupling) iomodel->DeleteData(2,"md.mesh.lat","md.mesh.long");
 }