Changeset 24579

Timestamp:

02/25/20 14:07:05 (5 years ago)

Author:

Mathieu Morlighem

Message:

CHG: making PDD inputs a dataset instead, as we know its size is 12

Location:

issm/trunk-jpl

Files:

• src/c/analyses/SmbAnalysis.cpp (modified) (1 diff)
• src/c/classes/Elements/Element.cpp (modified) (2 diffs)
• src/c/classes/Elements/Tria.cpp (modified) (2 diffs)
• src/c/classes/Inputs2/DatasetInput2.cpp (modified) (1 diff)
• src/c/classes/IoModel.cpp (modified) (2 diffs)
• src/c/classes/IoModel.h (modified) (1 diff)
• test/NightlyRun/test245.m (modified) (1 diff)
• test/NightlyRun/test245.py (modified) (1 diff)

issm/trunk-jpl/src/c/analyses/SmbAnalysis.cpp

@@ -96 +96 @@
                         iomodel->FindConstant(&isfirnwarming,"md.smb.isfirnwarming");
                         iomodel->FetchDataToInput(inputs2,elements,"md.smb.smb_corr",SmbSmbCorrEnum);
-                        iomodel->FetchDataToInput(inputs2,elements,"md.smb.precipitation",SmbPrecipitationEnum);
-                        iomodel->FetchDataToInput(inputs2,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
                         iomodel->FetchDataToInput(inputs2,elements,"md.smb.precipitation_anomaly",SmbPrecipitationsAnomalyEnum);
                         iomodel->FetchDataToInput(inputs2,elements,"md.smb.temperature_anomaly",SmbTemperaturesAnomalyEnum);
+                        iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.monthlytemperatures",SmbMonthlytemperaturesEnum);
+                        iomodel->FetchDataToDatasetInput(inputs2,elements,"md.smb.precipitation",SmbPrecipitationEnum);
                         break;
                 case SMBd18opddEnum:

issm/trunk-jpl/src/c/classes/Elements/Element.cpp

@@ -3030 +3030 @@
         IssmDouble mu           = MU_0*(1000.0*86400.0)*(rho_ice/rho_water);   // (d*deg C)/(mm WE) --> (s*deg C)/(m IE)
 
+        /*Get inputs*/
+        DatasetInput2* dinput =this->GetDatasetInput2(SmbMonthlytemperaturesEnum); _assert_(dinput);
+        DatasetInput2* dinput2=this->GetDatasetInput2(SmbPrecipitationEnum);       _assert_(dinput2);
 
         /*loop over vertices: */
@@ -3035 +3038 @@
         for(int month=0;month<12;month++){
 
-                /*Recover monthly temperatures and precipitation*/
-                Input2* dinput =this->GetInput2(SmbMonthlytemperaturesEnum,reCast<IssmDouble>(month+1)/12.*yts); _assert_(dinput);
-                Input2* dinput2=this->GetInput2(SmbPrecipitationEnum,reCast<IssmDouble>(month+1)/12.*yts);       _assert_(dinput2);
-
                 for(int iv=0;iv<NUM_VERTICES;iv++){
                         gauss->GaussVertex(iv);
-                        dinput->GetInputValue(&monthlytemperatures[iv*12+month],gauss);
+                        dinput->GetInputValue(&monthlytemperatures[iv*12+month],gauss,month);
                         monthlytemperatures[iv*12+month]=monthlytemperatures[iv*12+month]-273.15; // conversion from Kelvin to celcius for PDD module
-                        dinput2->GetInputValue(&monthlyprec[iv*12+month],gauss);
+                        dinput2->GetInputValue(&monthlyprec[iv*12+month],gauss,month);
                         monthlyprec[iv*12+month]=monthlyprec[iv*12+month]*yts;
                 }

issm/trunk-jpl/src/c/classes/Elements/Tria.cpp

@@ -233 +233 @@
 
         /*Intermediaries*/
+        int        vertexsids[NUMVERTICES];
+        int        vertexlids[NUMVERTICES];
         IssmDouble nodeinputs[NUMVERTICES];
-        int        vertexlids[NUMVERTICES];
 
         /*Some sanity checks*/
@@ -244 +245 @@
         _assert_(iomodel->elements);
         for(int i=0;i<NUMVERTICES;i++){
-                int vertexid =reCast<int>(iomodel->elements[NUMVERTICES*this->Sid()+i]); //ids for vertices are in the elements array from Matlab
-
-                vertexlids[i]=iomodel->my_vertices_lids[vertexid-1];
-                nodeinputs[i]=array[vertexid*N+i];
+                vertexsids[i] = reCast<int>(iomodel->elements[NUMVERTICES*this->Sid()+i])-1;
+                vertexlids[i] = iomodel->my_vertices_lids[vertexsids[i]];
         }
 
         /*Create inputs and add to DataSetInput*/
         for(int i=0;i<num_inputs;i++){
-                this->inputs2->SetTriaDatasetInput(input_enum,individual_enums[i],P1Enum,NUMVERTICES,vertexlids,nodeinputs); break;
+                for(int j=0;j<NUMVERTICES;j++) nodeinputs[j]=array[vertexsids[j]*N+i];
+                inputs2->SetTriaDatasetInput(input_enum,individual_enums[i],P1Enum,NUMVERTICES,vertexlids,nodeinputs);
         }
 }

issm/trunk-jpl/src/c/classes/Inputs2/DatasetInput2.cpp

@@ -207 +207 @@
         if(index == -1){
                 this->Echo();
-                _error_("Could not find input "<<id<<" ("<<EnumToStringx(id)<<") in DatasetInput");
+                _error_("Could not find input "<<id<<" ("<<EnumToStringx(id)<<"?) in DatasetInput");
         }
 

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

@@ -1663 +1663 @@
         int         integer;
         IssmDouble  scalar;
-        char       *string           = NULL;
         IssmDouble *doublearray = NULL;
         int         M,N;
@@ -1723 +1722 @@
         /*Free ressources*/
         xDelete<IssmDouble>(doublearray);
-        xDelete<char>(string);
+}
+/*}}}*/
+void  IoModel::FetchDataToDatasetInput(Inputs2* inputs2,Elements* elements,const char* vector_name,int input_enum){/*{{{*/
+
+        /*First, look whether it is not already loaded in this->data*/
+        vector<IoData*>::iterator iter;
+        for(iter=data.begin();iter<data.end();iter++){
+                IoData* iodata=*iter;
+                if(strcmp(iodata->name,vector_name)==0){
+                        for(int i=0;i<elements->Size();i++){
+                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                                _error_("to be implemented...");
+                                //element->InputCreate(iodata->data,inputs2,this,iodata->M,iodata->N,iodata->layout,input_enum,iodata->code);//we need i to index into elements.
+                        }
+                        return;
+                }
+        }
+
+        /*intermediary: */
+        int         code,vector_layout;
+        IssmDouble *doublearray = NULL;
+        int         M,N;
+
+        /*First of, find the record for the name, and get code  of data type: */
+        this->SetFilePointerToData(&code,&vector_layout,vector_name);
+
+        switch(code){
+                case 1: //boolean constant
+                        _error_("not implemented yet");
+                        break;
+                case 2: //integer constant
+                        _error_("not implemented yet");
+                        break;
+                case 3: //IssmDouble constant
+                        _error_("not implemented yet");
+                        break;
+                case 5: //boolean vector
+                        _error_("not implemented yet");
+                        break;
+                case 6: //int vector
+                        _error_("not implemented yet");
+                        break;
+                case 7: //IssmDouble vector
+                          {
+                        this->FetchData(&doublearray,&M,&N,vector_name);
+                        if(!doublearray) _error_("\""<<vector_name<<"\" not found in binary file");
+
+                        int* ids = xNew<int>(N);
+                        for(int i=0;i<N;i++) ids[i] = i;
+
+                        for(int i=0;i<elements->Size();i++){
+                                Element* element=xDynamicCast<Element*>(elements->GetObjectByOffset(i));
+                                element->DatasetInputCreate(doublearray,M,N,ids,N,inputs2,this,input_enum);
+                        }
+                        xDelete<int>(ids);
+                          }
+                        break;
+                default:
+                        _error_("data code " << code << " not supported yet (detected while processing \""<<vector_name<<"\")");
+                        break;
+        }
+        /*Free ressources*/
+        xDelete<IssmDouble>(doublearray);
 }
 /*}}}*/

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

@@ -141 +141 @@
                 void        FetchDataToInput(Inputs2* inputs2,Elements* elements,const char* vector_name,int input_enum);
                 void        FetchDataToInput(Inputs2* inputs2,Elements* elements,const char* vector_name,int input_enum,IssmDouble default_value);
+                void        FetchDataToDatasetInput(Inputs2* inputs2,Elements* elements,const char* vector_name,int input_enum);
                 void        FetchIndependent(const char* dependent_name);
                 void        FetchMultipleData(char***   pstringarray,int* pnumstrings,const char* data_name);

issm/trunk-jpl/test/NightlyRun/test245.m

@@ -18 +18 @@
 for imonth=0:11
     md.smb.monthlytemperatures(1:md.mesh.numberofvertices,imonth+1)=md.materials.meltingpoint+temp_ma_present+(temp_mj_present-temp_ma_present)*sin(double(imonth+1-4)*pi/6.0);
-    md.smb.monthlytemperatures(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
     md.smb.precipitation(1:md.mesh.numberofvertices,imonth+1)=precipitation;
-    md.smb.precipitation(md.mesh.numberofvertices+1,imonth+1)=((imonth+1)/12);
 end
 

issm/trunk-jpl/test/NightlyRun/test245.py

@@ -30 +30 @@
 for imonth in range(12):
     md.smb.monthlytemperatures[0:md.mesh.numberofvertices, imonth] = md.materials.meltingpoint + temp_ma_present + (temp_mj_present - temp_ma_present) * np.sin((imonth + 1. - 4.) * np.pi / 6.0)
-    md.smb.monthlytemperatures[md.mesh.numberofvertices, imonth] = ((imonth + 1) / 12.)
     md.smb.precipitation[0:md.mesh.numberofvertices, imonth] = precipitation
-    md.smb.precipitation[md.mesh.numberofvertices, imonth] = ((imonth + 1) / 12.)
 
 # time steps and resolution