source: issm/oecreview/Archive/21724-22754/ISSM-22748-22749.diff@ 22755

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

@@ -67 +67 @@
                                 md = checkfield(md,'fieldname','basalforcings.maxboxcount','numel',1,'NaN',1,'Inf',1,'>',0);
                                 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);
 

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

@@ -31 +31 @@
 void UpdateBoxIdsPico(FemModel* femmodel){/*{{{*/
 
         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);
         femmodel->DistanceToFieldValue(MaskGroundediceLevelsetEnum,0.,DistanceToGroundinglineEnum);
@@ -45 +46 @@
         femmodel->DistanceToFieldValue(MaskIceLevelsetEnum,0.,DistanceToCalvingfrontEnum);
 
         /*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));
                 if(!element->IsIceInElement() || !element->IsFloating()) continue;
@@ -55 +56 @@
                 element->GetInputListOnVertices(&distances[0],DistanceToGroundinglineEnum);
                 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*/
         int* nd=xNew<int>(num_basins);
         for(int i=0; i<num_basins;i++){
@@ -96 +102 @@
         femmodel->parameters->FindParam(&maxbox,BasalforcingsPicoMaxboxcountEnum);
         
         IssmDouble* boxareas=xNewZeroInit<IssmDouble>(num_basins*maxbox);
+
         for(int i=0;i<femmodel->elements->Size();i++){
                 Element* element=xDynamicCast<Element*>(femmodel->elements->GetObjectByOffset(i));
                 if(!element->IsIceInElement() || !element->IsFloating()) continue;
@@ -107 +114 @@
         /*Synchronize across cpus*/
         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));
 

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

@@ -233 +233 @@
                                 parameters->AddObject(iomodel->CopyConstantObject("md.basalforcings.overturning_coeff",BasalforcingsPicoOverturningCoeffEnum));
                                 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);
                         break;

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

@@ -2199 +2199 @@
                 xDelete<char*>(cost_functions);
         }
 }
+/*}}}*/
 void       Tria::InputUpdateFromSolutionOneDof(IssmDouble* solution,int enum_type){/*{{{*/
 
         /*Intermediary*/
@@ -2892 +2893 @@
         IssmDouble dist_gl,dist_cf;
 
         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);
         Input* dist_cf_input=inputs->GetInput(DistanceToCalvingfrontEnum);  _assert_(dist_cf_input);
@@ -2935 +2936 @@
         int        basinid, maxbox, num_basins, numnodes, M;
         IssmDouble gamma_T, overturning_coeff, thickness;
         IssmDouble pressure, T_star,p_coeff, q_coeff;
-        IssmDouble* boxareas  = NULL;
 
         /*Get variables*/
         IssmDouble rhoi       = this->GetMaterialParameter(MaterialsRhoIceEnum);
@@ -2957 +2957 @@
         this->parameters->FindParam(&gamma_T,BasalforcingsPicoGammaTEnum);
         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];
         IssmDouble g1               = area_boxi*gamma_T;
 

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

@@ -45 +45 @@
 field_names     ={'Bed1','Surface1','Thickness1','Floatingice1','Vx1','Vy1','Pressure1','FloatingiceMeltingrate1',...
            '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),...
            (md.results.TransientSolution(1).Surface),...