Changeset 22749

Timestamp:

05/07/18 14:29:32 (7 years ago)

Author:

tpelle

Message:

CHG: Update Pico floating ice melt rate parameterization

Location:

issm/trunk-jpl

Files:

• src/c/classes/Elements/Tria.cpp (modified) (4 diffs)
• src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp (modified) (5 diffs)
• src/c/modules/ModelProcessorx/CreateParameters.cpp (modified) (1 diff)
• src/m/classes/basalforcingspico.m (modified) (1 diff)
• test/Archives/Archive470.arch (modified) ( previous)
• test/NightlyRun/test470.m (modified) (1 diff)

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

@@ -2200 +2200 @@
         }
 }
+/*}}}*/
 void       Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){/*{{{*/
 
@@ -2893 +2894 @@
 
         inputs->GetInputValue(&basin_id,BasalforcingsPicoBasinIdEnum);
-        IssmDouble boxid_max=reCast<IssmDouble>(max_boxid_basin_list[basin_id]);
+        IssmDouble boxid_max=reCast<IssmDouble>(max_boxid_basin_list[basin_id])+1.;
 
         Input* dist_gl_input=inputs->GetInput(DistanceToGroundinglineEnum); _assert_(dist_gl_input);
@@ -2936 +2937 @@
         IssmDouble gamma_T, overturning_coeff, thickness;
         IssmDouble pressure, T_star,p_coeff, q_coeff;
-        IssmDouble* boxareas  = NULL;
 
         /*Get variables*/
@@ -2958 +2958 @@
         this->parameters->FindParam(&overturning_coeff,BasalforcingsPicoOverturningCoeffEnum);
         this->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
-        this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
         this->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
         Input* thickness_input=this->GetInput(ThicknessEnum); _assert_(thickness_input);
         _assert_(basinid<=num_basins);
+
+   IssmDouble* boxareas = xNew<IssmDouble>(num_basins*maxbox);
+        this->parameters->FindParam(&boxareas,&M,BasalforcingsPicoBoxAreaEnum);
 
         IssmDouble area_boxi        = boxareas[basinid*maxbox+boxid];

issm/trunk-jpl/src/c/modules/FloatingiceMeltingRatePicox/FloatingiceMeltingRatePicox.cpp

@@ -32 +32 @@
 
         int         numvertices,num_basins,maxbox,basinid;
-        IssmDouble* dmax_basin=NULL;
+        IssmDouble  dist_max;
+        IssmDouble* dmax_basin_cpu=NULL;
         IssmDouble* distances=NULL;
 
         femmodel->parameters->FindParam(&num_basins,BasalforcingsPicoNumBasinsEnum);
         femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
-        dmax_basin=xNew<IssmDouble>(num_basins);
+        dmax_basin_cpu=xNew<IssmDouble>(num_basins);
 
         femmodel->elements->InputDuplicate(MaskGroundediceLevelsetEnum,DistanceToGroundinglineEnum);
@@ -46 +47 @@
 
         /*find maximum distance to grounding line per domain and per basin*/
-        IssmDouble dist_max=-1.;
-        for(int i=0;i<num_basins;i++){dmax_basin[i]=-1;}
+        IssmDouble maxdist_cpu=-1.;
+        for(int i=0;i<num_basins;i++){dmax_basin_cpu[i]=-1;}
         for(int i=0;i<femmodel->elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
@@ -56 +57 @@
                 element->inputs->GetInputValue(&basinid,BasalforcingsPicoBasinIdEnum);
                 for(int k=0; k<numvertices; k++){
-                        if(fabs(distances[k])>dist_max){dist_max=fabs(distances[k]);}
-                        if(fabs(distances[k])>dmax_basin[basinid]){dmax_basin[basinid]=fabs(distances[k]);}
+                        if(fabs(distances[k])>maxdist_cpu){maxdist_cpu=fabs(distances[k]);}
+                        if(fabs(distances[k])>dmax_basin_cpu[basinid]){dmax_basin_cpu[basinid]=fabs(distances[k]);}
                 }
                 xDelete<IssmDouble>(distances);
         }
+
+        /*Synchronize across cpus*/
+        IssmDouble* dmax_basin=xNew<IssmDouble>(num_basins);
+        ISSM_MPI_Allreduce((void*)&maxdist_cpu,(void*)&dist_max,1,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
+        ISSM_MPI_Allreduce(dmax_basin_cpu,dmax_basin,num_basins,ISSM_MPI_DOUBLE,ISSM_MPI_MAX,IssmComm::GetComm());
 
         /*Define maximum number of boxes per basin*/
@@ -97 +103 @@
         
         IssmDouble* boxareas=xNewZeroInit<IssmDouble>(num_basins*maxbox);
+
         for(int i=0;i<femmodel->elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
@@ -108 +115 @@
         IssmDouble* sumareas =xNew<IssmDouble>(num_basins*maxbox);
         ISSM_MPI_Allreduce(boxareas,sumareas,num_basins*maxbox,ISSM_MPI_DOUBLE,ISSM_MPI_SUM,IssmComm::GetComm());
-
+        if(sumareas[0]==0){_error_("No elements in box 0, basal meltrates will be 0. Consider decreasing md.basalforcings.maxboxcount or refining your mesh!");}
+        
         /*Update parameters to keep track of the new areas in future calculations*/
         femmodel->parameters->AddObject(new DoubleVecParam(BasalforcingsPicoBoxAreaEnum,sumareas,maxbox*num_basins));

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

@@ -234 +234 @@
                                 parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.gamma_T",BasalforcingsPicoGammaTEnum));
                                 iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_temperature");
-                                _assert_(M>1 && N>1); 
+                                _assert_(M>=1 && N>=1); 
                                 parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceantemperatureEnum,transparam,&transparam[N*(M-1)],interp,N,M));
                                 xDelete<IssmDouble>(transparam);
                                 iomodel->FetchData(&transparam,&M,&N,"md.basalforcings.farocean_salinity");
-                                _assert_(M>1 && N>1); 
+                                _assert_(M>=1 && N>=1); 
                                 parameters->AddObject(new TransientArrayParam(BasalforcingsPicoFarOceansalinityEnum,transparam,&transparam[N*(M-1)],interp,N,M));
                                 xDelete<IssmDouble>(transparam);

issm/trunk-jpl/src/m/classes/basalforcingspico.m

@@ -68 +68 @@
                                 md = checkfield(md,'fieldname','basalforcings.overturning_coeff','numel',1,'NaN',1,'Inf',1,'>',0);
                                 md = checkfield(md,'fieldname','basalforcings.gamma_T','numel',1,'NaN',1,'Inf',1,'>',0);
-                                md = checkfield(md,'fieldname','basalforcings.farocean_temperature','NaN',1,'Inf',1,'>',0);
-                                md = checkfield(md,'fieldname','basalforcings.farocean_salinity','NaN',1,'Inf',1,'>',0);
+                                md = checkfield(md,'fieldname','basalforcings.farocean_temperature','NaN',1,'Inf',1,'>',0,'size',[md.basalforcings.num_basins+1 NaN]);
+                                md = checkfield(md,'fieldname','basalforcings.farocean_salinity','NaN',1,'Inf',1,'>',0,'size',[md.basalforcings.num_basins+1 NaN]);
                                 md = checkfield(md,'fieldname','basalforcings.geothermalflux','NaN',1,'Inf',1,'>=',0,'timeseries',1);
                                 md = checkfield(md,'fieldname','basalforcings.groundedice_melting_rate','NaN',1,'Inf',1,'timeseries',1);

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

@@ -46 +46 @@
            'Bed2','Surface2','Thickness2','Floatingice2','Vx2','Vy2','Pressure2','FloatingiceMeltingrate2',...
            'Bed3','Surface3','Thickness3','Floatingice3','Vx3','Vy3','Pressure3','FloatingiceMeltingrate3'};
-field_tolerances={1,1,1,1,1,1,1,1,...
-           1,1,1,1,1,1,1,1,...
-           1,1,1,1,1,1,1,1};
+field_tolerances={2e-11,5e-12,2e-11,1e-11,5e-10,1e-08,1e-13,1e-13,...
+           3e-11,3e-11,9e-10,7e-11,1e-09,5e-08,1e-10,1e-13,...
+           1e-10,3e-11,1e-10,7e-11,1e-09,5e-08,1e-10,1e-13};
 field_values={...
            (md.results.TransientSolution(1).Base),...