Changeset 24173

Timestamp:

Sep 27, 2019, 4:44:53 PM (7 years ago)

Author:

schlegel

Message:

CHG: added v and e partitions

Location:

issm/trunk-jpl/src/m

Files:

• classes/qmu.m (modified) (3 diffs)
• classes/qmu.py (modified) (4 diffs)
• partition/partitioner.m (modified) (7 diffs)
• partition/partitioner.py (modified) (4 diffs)
• qmu/importancefactors.m (modified) (1 diff)
• qmu/importancefactors.py (modified) (1 diff)

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

@@ -66 +66 @@
                                 end
                         end
-                        if ~isempty(md.qmu.vpartition),
-                                if numel(md.qmu.vpartition)~=md.mesh.numberofvertices 
+                        if ~isempty(md.qmu.vpartition) & ~any(isnan(md.qmu.vpartition)),
+                                if numel(md.qmu.vpartition)~=md.mesh.numberofvertices  
                                         md = checkmessage(md,['user supplied vertex partition for qmu analysis should have size md.mesh.numberofvertices x 1']);
                                 end
@@ -77 +77 @@
                                 end
                         end
-                        if ~isempty(md.qmu.epartition),
+                        if ~isempty(md.qmu.epartition) & ~any(isnan(md.qmu.epartition)),
                                 if numel(md.qmu.epartition)~=md.mesh.numberofelements, 
                                         md = checkmessage(md,['user supplied element partition for qmu analysis should have size md.mesh.numberofelements x 1']);
@@ -88 +88 @@
                                 end
                         end
-                        if isempty(md.qmu.vpartition) & isempty(md.qmu.epartition),
-                                        md = checkmessage(md,['for qmu analysis, at least one of element and partitioning vectors need to be supplied!']);
-                        end
+                        if isempty(md.qmu.vpartition) | any(isnan(md.qmu.vpartition)) | isempty(md.qmu.epartition) | any(isnan(md.qmu.epartition)),
+                                md = checkmessage(md,['for qmu analysis, both an element and partitioning vectors need to be supplied with no nan values! One can be defaulted to all zeros.']);
+                        end
+
                 end % }}}
                 function disp(self) % {{{

issm/trunk-jpl/src/m/classes/qmu.py

@@ -25 +25 @@
                 self.params                      = OrderedStruct()
                 self.results                     = OrderedDict()
-                self.partition                   = float('NaN')
+                self.vpartition                  = float('NaN')
+                self.epartition                  = float('NaN')
                 self.numberofpartitions          = 0
                 self.numberofresponses           = 0
@@ -107 +108 @@
                                 s+="            %-*s:    [%ix%i]    '%s'\n" % (maxlen+1,fname,a,b,type(getattr(result,fname)))
 
-                s+="%s\n" % fielddisplay(self,'partition','user provided mesh partitioning, defaults to metis if not specified') 
+                s+="%s\n" % fielddisplay(self,'vpartition','user provided mesh partitioning (vertex based)') 
+                s+="%s\n" % fielddisplay(self,'epartition','user provided mesh partitioning (element based)')
                 s+="%s\n" % fielddisplay(self,'numberofpartitions','number of partitions for semi-discrete qmu') 
                 s+="%s\n" % fielddisplay(self,'variabledescriptors','')
@@ -152 +154 @@
                                 md.checkmessage('in parallel library mode, the requirement is for md.cluster.np = md.qmu.params.processors_per_evaluation * number_of_slaves, where number_of_slaves will automatically be determined by Dakota. Modify md.cluster.np accordingly')
                 
-                if np.size(md.qmu.partition) > 0:
-                        if (np.size(md.qmu.partition)!=md.mesh.numberofvertices and np.size(md.qmu.partition) != md.mesh.numberofelements) or np.size(md.qmu.partition,1) != 1:
-                                md.checkmessage("user supplied partition for qmu analysis should have size (md.mesh.numberofvertices x 1) or (md.mesh.numberofelements x 1)")
-                        if not min(md.qmu.partition.flatten())==0:
-                                md.checkmessage("partition vector not indexed from 0 on")
-                        if max(md.qmu.partition.flatten())>=md.qmu.numberofpartitions:
-                                md.checkmessage("for qmu analysis, partitioning vector cannot go over npart, number of partition areas")
+                if np.size(md.qmu.vpartition) > 0:
+                        if np.size(md.qmu.vpartition)!=md.mesh.numberofvertices:
+                                md.checkmessage("user supplied vertex partition for qmu analysis should have size (md.mesh.numberofvertices x 1)")
+                        if not min(md.qmu.vpartition.flatten())==0:
+                                md.checkmessage("vertex partition vector not indexed from 0 on")
+                        if max(md.qmu.vpartition.flatten())>=md.qmu.numberofpartitions:
+                                md.checkmessage("for qmu analysis, vertex partitioning vector cannot go over npart, number of partition areas")
+
+                if np.size(md.qmu.epartition) > 0:
+                        if np.size(md.qmu.epartition) != md.mesh.numberofelements:
+                                md.checkmessage("user supplied element partition for qmu analysis should have size (md.mesh.numberofelements x 1)")
+                        if not min(md.qmu.epartition.flatten())==0:
+                                md.checkmessage("elememtn partition vector not indexed from 0 on")
+                        if max(md.qmu.epartition.flatten())>=md.qmu.numberofpartitions:
+                                md.checkmessage("for qmu analysis, element partitioning vector cannot go over npart, number of partition areas")
+
+                if np.size(md.qmu.vpartition) == 0 or np.any(np.isnan(md.qmu.vpartition)) or np.size(md.qmu.epartition) == 0 or np.any(np.isnan(md.qmu.epartition)):
+                        md.checkmessage("for qmu analysis, both an element and partitioning vectors need to be supplied with no nan values! One can be defaulted to all zeros.")
 
                 return md
@@ -167 +180 @@
                         WriteData(fid,prefix,'data',False,'name','md.qmu.mass_flux_segments_present','format','Boolean');
                         return
-                WriteData(fid,prefix,'object',self,'fieldname','partition','format','DoubleMat','mattype',2)
+                WriteData(fid,prefix,'object',self,'fieldname','vpartition','format','DoubleMat','mattype',2)
+                WriteData(fid,prefix,'object',self,'fieldname','epartition','format','DoubleMat','mattype',2)
                 WriteData(fid,prefix,'object',self,'fieldname','numberofpartitions','format','Integer')
                 WriteData(fid,prefix,'object',self,'fieldname','numberofresponses','format','Integer')

issm/trunk-jpl/src/m/partition/partitioner.m

@@ -9 +9 @@
 %   section:  1 by defaults(1=bisection, 2=quadrisection, 3=octasection)
 %   recomputeadjacency:  'on' by default (set to 'off' to compute existing one)
+%   type: 'node' or 'element' partition vector (default to 'node') 
 %   Output: md.qmu.partition recover the partition vector
 %   
@@ -24 +25 @@
 options=addfielddefault(options,'section',1);
 options=addfielddefault(options,'recomputeadjacency','on');
+options=addfielddefault(options,'type','node');
 
 %get package: 
@@ -29 +31 @@
 npart=getfieldvalue(options,'npart');
 recomputeadjacency=getfieldvalue(options,'recomputeadjacency');
+vectortype=getfieldvalue(options,'type');
 
 if(dimension(md.mesh)==3),
@@ -53 +56 @@
 if strcmpi(package,'chaco'),
 
-        %  default method (from chaco.m)
-        method=[1 1 0 0 1 1 50 0 .001 7654321]';
-        method(1)=3;    %  global method (3=inertial (geometric))
-        method(3)=0;    %  vertex weights (0=off, 1=on)
+        if strcmpi(vectortype,'element')
+                error(['partitioner error message: package ' package ' does not allow element partitions.']);
+        else
 
-        %specify bisection
-        method(6)=getfieldvalue(options,'section');%  ndims (1=bisection, 2=quadrisection, 3=octasection)
+                %  default method (from chaco.m)
+                method=[1 1 0 0 1 1 50 0 .001 7654321]';
+                method(1)=3;    %  global method (3=inertial (geometric))
+                method(3)=0;    %  vertex weights (0=off, 1=on)
 
-        %are we using weights? 
-        if strcmpi(getfieldvalue(options,'weighting'),'on'),
-                weights=floor(md.qmu.vertex_weight/min(md.qmu.vertex_weight));
-                method(3)=1;
-        else 
-                weights=[];
-        end
+                %specify bisection
+                method(6)=getfieldvalue(options,'section');%  ndims (1=bisection, 2=quadrisection, 3=octasection)
 
-        %  partition into nparts
-        if isa(md.mesh,'mesh2d'),
-                part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,zeros(md.mesh.numberofvertices,1),method,npart,[])'+1; %index partitions from 1 up. like metis.
-        else
-                part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,md.mesh.z,method,npart,[])'+1; %index partitions from 1 up. like metis.
+                %are we using weights? 
+                if strcmpi(getfieldvalue(options,'weighting'),'on'),
+                        weights=floor(md.qmu.vertex_weight/min(md.qmu.vertex_weight));
+                        method(3)=1;
+                else 
+                        weights=[];
+                end
+
+                %  partition into nparts
+                if isa(md.mesh,'mesh2d'),
+                        part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,zeros(md.mesh.numberofvertices,1),method,npart,[])'+1; %index partitions from 1 up. like metis.
+                else
+                        part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,md.mesh.z,method,npart,[])'+1; %index partitions from 1 up. like metis.
+                end
+
         end
 
 elseif strcmpi(package,'scotch'),
 
-        %are we using weights? 
-        if strcmpi(getfieldvalue(options,'weighting'),'on'),
-                weights=floor(md.qmu.vertex_weight/min(md.qmu.vertex_weight));
+        if strcmpi(vectortype,'element')
+                error(['partitioner error message: package ' package ' does not allow element partitions.']);
         else
-                weights=[];
+                %are we using weights? 
+                if strcmpi(getfieldvalue(options,'weighting'),'on'),
+                        weights=floor(md.qmu.vertex_weight/min(md.qmu.vertex_weight));
+                else
+                        weights=[];
+                end
+                maptab=Scotch(md.qmu.adjacency,[],weights,[],'cmplt',[npart]);
+
+                part=maptab(:,2)+1;%index partitions from 1 up. like metis.
         end
-        maptab=Scotch(md.qmu.adjacency,[],weights,[],'cmplt',[npart]);
-
-        part=maptab(:,2)+1;%index partitions from 1 up. like metis.
 
 elseif strcmpi(package,'linear'),
 
-        if npart==md.mesh.numberofelements | md.qmu.numberofpartitions==md.mesh.numberofelements
+        if strcmpi(vectortype,'element')
                 part=1:1:md.mesh.numberofelements;
                 disp('Linear partitioner requesting partitions on elements');
@@ -99 +112 @@
 elseif strcmpi(package,'metis'),
 
-        [element_partitioning,part]=MeshPartition(md,md.qmu.numberofpartitions);
+        if strcmpi(vectortype,'element')
+                error(['partitioner error message: package ' package ' does not allow element partitions.']);
+        else
+                [element_partitioning,part]=MeshPartition(md,md.qmu.numberofpartitions);
+        end
 
 else
@@ -112 +129 @@
         md3d.qmu.adjacency=md.qmu.adjacency;
         md=md3d;
-        part=project3d(md,'vector',part','type','node');
+        if strcmpi(vectortype,'element')
+                part=project3d(md,'vector',part','type','element');
+        else
+                part=project3d(md,'vector',part','type','node');
+        end
 end
 
@@ -118 +139 @@
         part=part';
 end
-md.qmu.partition=part;
+
+if strcmpi(vectortype,'element')
+        md.qmu.epartition=part;
+        if isempty(md.qmu.vpartition) | (prod(size(md.qmu.vpartition))==1 & isnan(md.qmu.vpartition))
+                md.qmu.vpartition=zeros(md.mesh.numberofvertices,1);
+        end
+else
+        md.qmu.vpartition=part;
+        if isempty(md.qmu.epartition) | (prod(size(md.qmu.epartition))==1 & isnan(md.qmu.epartition))
+                md.qmu.epartition=zeros(md.mesh.numberofelements,1);
+        end
+end

issm/trunk-jpl/src/m/partition/partitioner.py

@@ -21 +21 @@
    section:  1 by defaults(1=bisection, 2=quadrisection, 3=octasection)
    recomputeadjacency:  'on' by default (set to 'off' to compute existing one)
+   type: 'node' or 'element' partition vector (default to 'node')
    Output: md.qmu.partition recover the partition vector
    
@@ -36 +37 @@
         options.addfielddefault('section',1)
         options.addfielddefault('recomputeadjacency','on')
+        options.addfielddefault('type','node')
 
         #get package: 
@@ -41 +43 @@
         npart=options.getfieldvalue('npart')
         recomputeadjacency=options.getfieldvalue('recomputeadjacency')
+        vectortype=options.getfieldvalue('type')
+
+        if m.strcmpi(vectortype,'element') and not m.strcmpi(package,'linear'):
+                raise RuntimeError('partitioner error message: package '+str(package)+' does not allow element partitions.')
 
         if(md.mesh.dimension()==3):
@@ -121 +127 @@
                 md3d.qmu.adjacency=md.qmu.adjacency
                 md=md3d
-                part=project3d(md,'vector',np.squeeze(part),'type','node')
+                if m.strcmpi(vectortype,'element'):
+                    part=project3d(md,'vector',np.squeeze(part),'type','element')
+                else:
+                    part=project3d(md,'vector',np.squeeze(part),'type','node')
 
         part=part.reshape(-1,1)
-        md.qmu.partition=part
+
+        if m.strcmpi(vectortype,'element'):
+                md.qmu.epartition=part
+                if np.size(md.qmu.vpartition) == 0 or (np.size(md.qmu.vpartition) == 1 and np.isnan(md.qmu.vpartition)): 
+                    md.qmu.vpartition=np.zeros((md.mesh.numberofvertices,1))
+
+        else:
+                md.qmu.vpartition=part
+                if np.size(md.qmu.epartition) == 0 or (np.size(md.qmu.epartition) == 1 and np.isnan(md.qmu.epartition)):
+                    md.qmu.epartition=np.zeros((md.mesh.numberofelements,1))
 
         return md

issm/trunk-jpl/src/m/qmu/importancefactors.m

@@ -44 +44 @@
         factors=importancefactors;
         return;
+elseif count==max(md.qmu.epartition+1)
+        %distribute importance factor
+        factors=importancefactors(md.qmu.epartition'+1); %md.qmu.partition was created to index "c" style
 else
         %distribute importance factor
-        factors=importancefactors(md.qmu.partition'+1); %md.qmu.partition was created to index "c" style
+        factors=importancefactors(md.qmu.vpartition'+1); %md.qmu.partition was created to index "c" style
 end
 

issm/trunk-jpl/src/m/qmu/importancefactors.py

@@ -43 +43 @@
                 factors=importancefactors
                 return factors
+        elif count==np.max(md.qmu.epartition+1):
+                #distribute importance factor
+                factors=importancefactors[(md.qmu.epartition.conj().T).flatten().astype(int)]
+                #md.qmu.partition was created to index "c" style
         else:
                 #distribute importance factor
-                factors=importancefactors[(md.qmu.partition.conj().T).flatten().astype(int)]
+                factors=importancefactors[(md.qmu.vpartition.conj().T).flatten().astype(int)]
                 #md.qmu.partition was created to index "c" style