Changeset 17422

Timestamp:

03/12/14 11:20:26 (11 years ago)

Author:

bdef

Message:

NEW: Adding new stuff to deal with moulin input

Location:

issm/trunk-jpl

Files:

• src/c/analyses/HydrologyDCInefficientAnalysis.cpp (modified) (4 diffs)
• src/c/shared/Enum/EnumDefinitions.h (modified) (1 diff)
• src/c/shared/Enum/EnumToStringx.cpp (modified) (1 diff)
• src/c/shared/Enum/StringToEnumx.cpp (modified) (5 diffs)
• src/m/classes/hydrologydc.m (modified) (5 diffs)
• src/m/classes/model.m (modified) (2 diffs)
• src/m/contrib/paraview/exportVTK.m (modified) (5 diffs)
• src/m/enum/EnumDefinitions.py (modified) (1 diff)
• test/NightlyRun/test3300.m (modified) (1 diff)
• test/NightlyRun/test332.m (modified) (1 diff)
• test/NightlyRun/test333.m (modified) (1 diff)
• test/NightlyRun/test334.m (modified) (1 diff)
• test/NightlyRun/test335.m (modified) (1 diff)

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

@@ -86 +86 @@
         iomodel->FetchDataToInput(elements,MeshElementonbedEnum);
         iomodel->FetchDataToInput(elements,BasalforcingsMeltingRateEnum);
+        iomodel->FetchDataToInput(elements,HydrologydcBasalMoulinInputEnum);
         iomodel->FetchDataToInput(elements,SedimentHeadEnum);
         iomodel->FetchDataToInput(elements,HydrologydcSedimentTransmitivityEnum);
@@ -278 +279 @@
         /*Intermediaries */
         bool       active_element,isefficientlayer;
-        IssmDouble dt,scalar,water_head;
+        IssmDouble dt,scalar;
+        IssmDouble moulin_load,water_head;
         IssmDouble water_load,transfer;
         IssmDouble Jdet;
@@ -304 +306 @@
         Input* thickness_input   = basalelement->GetInput(ThicknessEnum);
         Input* bed_input         = basalelement->GetInput(BedEnum);
-        Input* water_input       = basalelement->GetInput(BasalforcingsMeltingRateEnum); _assert_(water_input);
-        if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);          _assert_(old_wh_input);}
+        Input* water_input       = basalelement->GetInput(BasalforcingsMeltingRateEnum);    _assert_(water_input);
+        Input* moulin_input      = basalelement->GetInput(HydrologydcBasalMoulinInputEnum); _assert_(moulin_input);
+        if(dt!= 0.){old_wh_input = basalelement->GetInput(SedimentHeadOldEnum);             _assert_(old_wh_input);}
 
         IssmDouble sediment_storing    = SedimentStoring(basalelement);
@@ -324 +327 @@
                 /*Loading term*/
                 water_input->GetInputValue(&water_load,gauss);
+                moulin_input->GetInputValue(&moulin_load,gauss);
         
                 scalar = Jdet*gauss->weight*(water_load);
+                scalar = scalar + Jdet* moulin_load;
                 if(dt!=0.) scalar = scalar*dt;
                 for(int i=0;i<numnodes;i++){

issm/trunk-jpl/src/c/shared/Enum/EnumDefinitions.h

@@ -126 +126 @@
         HydrologydcPenaltyFactorEnum,
         HydrologydcPenaltyLockEnum,
+        HydrologydcBasalMoulinInputEnum,
         HydrologyLayerEnum,
         HydrologySedimentEnum,

issm/trunk-jpl/src/c/shared/Enum/EnumToStringx.cpp

@@ -134 +134 @@
                 case HydrologydcPenaltyFactorEnum : return "HydrologydcPenaltyFactor";
                 case HydrologydcPenaltyLockEnum : return "HydrologydcPenaltyLock";
+                case HydrologydcBasalMoulinInputEnum : return "HydrologydcBasalMoulinInput";
                 case HydrologyLayerEnum : return "HydrologyLayer";
                 case HydrologySedimentEnum : return "HydrologySediment";

issm/trunk-jpl/src/c/shared/Enum/StringToEnumx.cpp

@@ -134 +134 @@
               else if (strcmp(name,"HydrologydcPenaltyFactor")==0) return HydrologydcPenaltyFactorEnum;
               else if (strcmp(name,"HydrologydcPenaltyLock")==0) return HydrologydcPenaltyLockEnum;
+              else if (strcmp(name,"HydrologydcBasalMoulinInput")==0) return HydrologydcBasalMoulinInputEnum;
               else if (strcmp(name,"HydrologyLayer")==0) return HydrologyLayerEnum;
               else if (strcmp(name,"HydrologySediment")==0) return HydrologySedimentEnum;
-              else if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum;
          else stage=2;
    }
    if(stage==2){
-              if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
+              if (strcmp(name,"HydrologyEfficient")==0) return HydrologyEfficientEnum;
+              else if (strcmp(name,"HydrologySedimentKmax")==0) return HydrologySedimentKmaxEnum;
               else if (strcmp(name,"IndependentObject")==0) return IndependentObjectEnum;
               else if (strcmp(name,"InversionControlParameters")==0) return InversionControlParametersEnum;
@@ -259 +260 @@
               else if (strcmp(name,"MaxIterationConvergenceFlag")==0) return MaxIterationConvergenceFlagEnum;
               else if (strcmp(name,"SteadystateMaxiter")==0) return SteadystateMaxiterEnum;
-              else if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
          else stage=3;
    }
    if(stage==3){
-              if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
+              if (strcmp(name,"SteadystateNumRequestedOutputs")==0) return SteadystateNumRequestedOutputsEnum;
+              else if (strcmp(name,"SteadystateReltol")==0) return SteadystateReltolEnum;
               else if (strcmp(name,"SteadystateRequestedOutputs")==0) return SteadystateRequestedOutputsEnum;
               else if (strcmp(name,"Surface")==0) return SurfaceEnum;
@@ -382 +383 @@
               else if (strcmp(name,"Nodes")==0) return NodesEnum;
               else if (strcmp(name,"Contours")==0) return ContoursEnum;
-              else if (strcmp(name,"Parameters")==0) return ParametersEnum;
          else stage=4;
    }
    if(stage==4){
-              if (strcmp(name,"Vertices")==0) return VerticesEnum;
+              if (strcmp(name,"Parameters")==0) return ParametersEnum;
+              else if (strcmp(name,"Vertices")==0) return VerticesEnum;
               else if (strcmp(name,"Results")==0) return ResultsEnum;
               else if (strcmp(name,"GenericParam")==0) return GenericParamEnum;
@@ -505 +506 @@
               else if (strcmp(name,"VzSSA")==0) return VzSSAEnum;
               else if (strcmp(name,"VzHO")==0) return VzHOEnum;
-              else if (strcmp(name,"VzPicard")==0) return VzPicardEnum;
          else stage=5;
    }
    if(stage==5){
-              if (strcmp(name,"VzFS")==0) return VzFSEnum;
+              if (strcmp(name,"VzPicard")==0) return VzPicardEnum;
+              else if (strcmp(name,"VzFS")==0) return VzFSEnum;
               else if (strcmp(name,"VxMesh")==0) return VxMeshEnum;
               else if (strcmp(name,"VyMesh")==0) return VyMeshEnum;
@@ -628 +629 @@
               else if (strcmp(name,"ToolkitsOptionsAnalyses")==0) return ToolkitsOptionsAnalysesEnum;
               else if (strcmp(name,"ToolkitsOptionsStrings")==0) return ToolkitsOptionsStringsEnum;
-              else if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
          else stage=6;
    }
    if(stage==6){
-              if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
+              if (strcmp(name,"QmuErrName")==0) return QmuErrNameEnum;
+              else if (strcmp(name,"QmuInName")==0) return QmuInNameEnum;
               else if (strcmp(name,"QmuOutName")==0) return QmuOutNameEnum;
               else if (strcmp(name,"Regular")==0) return RegularEnum;

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

@@ -16 +16 @@
                 transfer_flag            = 0;
                 leakage_factor           = 0;
+                basal_moulin_input       = NaN;
 
                 spcsediment_head         = NaN;
@@ -42 +43 @@
                 end 
                 % }}}
+                function self = initialize(self,md) % {{{
+                        if isnan(self.basal_moulin_input),
+                                self.basal_moulin_input=zeros(md.mesh.numberofvertices,1);
+                                disp('      no hydrology.basal_moulin_input specified: values set as zero');
+                        end
+
+                end % }}}
                 % {{{ function obj = setdefaultparameters(obj) 
                 function obj = setdefaultparameters(obj) 
@@ -88 +96 @@
                                 md = checkfield(md,'fieldname','hydrology.leakage_factor','>',0,'numel',1);
             end
+                        md = checkfield(md,'fieldname','hydrology.basal_moulin_input','NaN',1,'forcing',1);
 
                         md = checkfield(md,'fieldname','hydrology.spcsediment_head','forcing',1);
@@ -129 +138 @@
                                 fielddisplay(obj,'leakage_factor','user defined leakage factor [m]');
             end
+                        fielddisplay(obj,'basal_moulin_input','Figure out what it is');
                         disp(sprintf('   - for the sediment layer'));
                         fielddisplay(obj,'spcsediment_head','sediment water head constraints (NaN means no constraint) [m above MSL]');
@@ -165 +175 @@
                                 WriteData(fid,'object',obj,'fieldname','leakage_factor','format','Double');
             end
+                        WriteData(fid,'object',obj,'fieldname','basal_moulin_input','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1)
 
                         WriteData(fid,'object',obj,'fieldname','spcsediment_head','format','DoubleMat','mattype',1,'forcinglength',md.mesh.numberofvertices+1);

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

@@ -727 +727 @@
                         md.stressbalance.referential=project3d(md,'vector',md.stressbalance.referential,'type','node');
                         md.stressbalance.loadingforce=project3d(md,'vector',md.stressbalance.loadingforce,'type','node');
+
+                        % Hydrologydc variables
                         if isa(md.hydrology,'hydrologydc');
                                 md.hydrology.spcsediment_head=project3d(md,'vector',md.hydrology.spcsediment_head,'type','node','layer',1);
@@ -732 +734 @@
                                 md.hydrology.mask_eplactive_node=project3d(md,'vector',md.hydrology.mask_eplactive_node,'type','node','layer',1);
                                 md.hydrology.sediment_transmitivity=project3d(md,'vector',md.hydrology.sediment_transmitivity,'type','node','layer',1);
+                                md.hydrology.basal_moulin_input=project3d(md,'vector',md.hydrology.basal_moulin_input,'type','node','layer',1);
             end
 

issm/trunk-jpl/src/m/contrib/paraview/exportVTK.m

@@ -60 +60 @@
         timestep=step;
 
-        FID = fopen(strcat(path,filesep,name,filesep,name,'.vtk',int2str(timestep),'.vtk'),'w+');
-        fprintf(FID,'# vtk DataFile Version 2.0 \n');
-        fprintf(FID,'Data for run %s \n',model.miscellaneous.name);
-        fprintf(FID,'ASCII \n');
-        fprintf(FID,'DATASET UNSTRUCTURED_GRID \n');
+        fid = fopen(strcat(path,filesep,name,filesep,name,'.vtk',int2str(timestep),'.vtk'),'w+');
+        fprintf(fid,'# vtk DataFile Version 2.0 \n');
+        fprintf(fid,'Data for run %s \n',model.miscellaneous.name);
+        fprintf(fid,'ASCII \n');
+        fprintf(fid,'DATASET UNSTRUCTURED_GRID \n');
         
-        fprintf(FID,'POINTS %d float\n',num_of_points);
+        fprintf(fid,'POINTS %d float\n',num_of_points);
         if(dim==3);
                 s='%f %f %f \n';
@@ -73 +73 @@
   end
         P=[points zeros(num_of_points,3-dim)];
-        fprintf(FID,s,P');
+        fprintf(fid,s,P');
         
-        fprintf(FID,'CELLS %d %d\n',num_of_elt,num_of_elt*(point_per_elt+1));
+        fprintf(fid,'CELLS %d %d\n',num_of_elt,num_of_elt*(point_per_elt+1));
         s='%d';
         for j=1:point_per_elt
@@ -81 +81 @@
   end
         s=cell2mat(horzcat(s,{'\n'}));
-        fprintf(FID,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]');
+        fprintf(fid,s,[(point_per_elt)*ones(num_of_elt,1) model.mesh.elements-1]');
         
-        fprintf(FID,'CELL_TYPES %d\n',num_of_elt);
+        fprintf(fid,'CELL_TYPES %d\n',num_of_elt);
         s='%d\n';
-        fprintf(FID,s,celltype*ones(num_of_elt,1));
-        fprintf(FID,'POINT_DATA %s \n',num2str(num_of_points));
+        fprintf(fid,s,celltype*ones(num_of_elt,1));
+        fprintf(fid,'POINT_DATA %s \n',num2str(num_of_points));
 
         %loop over the different solution structures
@@ -112 +112 @@
                                         smallval=(abs(sol_struct{j}(timestep).(fieldnames{k}))<1.0e-20);
                                         sol_struct{j}(timestep).(fieldnames{k})(smallval)=0.0;
-                                        fprintf(FID,'SCALARS %s float 1 \n',fieldnames{k});
-                                        fprintf(FID,'LOOKUP_TABLE default\n');
+                                        fprintf(fid,'SCALARS %s float 1 \n',fieldnames{k});
+                                        fprintf(fid,'LOOKUP_TABLE default\n');
                                         s='%e\n';
-                                        fprintf(FID,s,sol_struct{j}(timestep).(fieldnames{k}));
+                                        fprintf(fid,s,sol_struct{j}(timestep).(fieldnames{k}));
                     end         
             end 
@@ -134 +134 @@
                                 smallval=(abs(res_struct.(fieldnames{k}))<1.0e-20);
                                 res_struct.(fieldnames{k})(smallval)=0.0;
-                                fprintf(FID,'SCALARS %s float 1 \n',fieldnames{k});
-                                fprintf(FID,'LOOKUP_TABLE default\n');
+                                fprintf(fid,'SCALARS %s float 1 \n',fieldnames{k});
+                                fprintf(fid,'LOOKUP_TABLE default\n');
                                 s='%e\n';
-                                fprintf(FID,s,res_struct.(fieldnames{k}));
+                                fprintf(fid,s,res_struct.(fieldnames{k}));
             end         
                 end 
         end
-        fclose(FID);
+        fclose(fid);
 end

issm/trunk-jpl/src/m/enum/EnumDefinitions.py

@@ -126 +126 @@
 def HydrologydcPenaltyFactorEnum(): return StringToEnum("HydrologydcPenaltyFactor")[0]
 def HydrologydcPenaltyLockEnum(): return StringToEnum("HydrologydcPenaltyLock")[0]
+def HydrologydcBasalMoulinInputEnum(): return StringToEnum("HydrologydcBasalMoulinInput")[0]
 def HydrologyLayerEnum(): return StringToEnum("HydrologyLayer")[0]
 def HydrologySedimentEnum(): return StringToEnum("HydrologySediment")[0]

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

@@ -5 +5 @@
 md.cluster=generic('name',oshostname(),'np',1);
 md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
 md.hydrology.isefficientlayer=1;
 md.hydrology.sedimentlimit_flag=1;

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

@@ -5 +5 @@
 md.cluster=generic('name',oshostname(),'np',1);
 md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
 md.hydrology.isefficientlayer=0;
 md.hydrology.sedimentlimit_flag=1;

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

@@ -5 +5 @@
 md.cluster=generic('name',oshostname(),'np',1);
 md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
 md.hydrology.isefficientlayer=1;
 md.hydrology.sedimentlimit_flag=1;

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

@@ -5 +5 @@
 md.cluster=generic('name',oshostname(),'np',1);
 md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
 md.hydrology.isefficientlayer=0;
 md.hydrology.sedimentlimit_flag=1;

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

@@ -5 +5 @@
 md.cluster=generic('name',oshostname(),'np',1);
 md.hydrology=(hydrologydc);
+md.hydrology=initialize(md.hydrology,md);
 md.hydrology.isefficientlayer=1;
 md.hydrology.sedimentlimit_flag=1;