![(please configure the [header_logo] section in trac.ini)](/trac/issm/chrome/common/trac_banner.png)
Changeset 26901
- Timestamp:
- 02/25/22 00:11:37 (3 years ago)
- Location:
- issm/trunk-jpl
- Files:
-
- 4 edited
-
src/m/contrib/defleurian/netCDF/export_netCDF.py (modified) (5 diffs)
-
src/m/contrib/defleurian/netCDF/restable.m (modified) (3 diffs)
-
src/m/io/loadvars.py (modified) (12 diffs)
-
test/NightlyRun/runme.py (modified) (3 diffs)
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py
r26871 r26901 1 1 from netCDF4 import Dataset 2 2 import numpy as np 3 import numpy.ma as ma 3 4 import time 4 5 import collections … … 30 31 datasize = np.squeeze(self.sizes) 31 32 maxsize = [] 32 33 33 try: 34 34 dims = np.arange(np.shape(datasize)[1]) … … 37 37 except IndexError: 38 38 maxsize.append(np.nanmax(datasize[:])) 39 40 39 findim = np.insert(maxsize, 0, rows) 41 42 40 #first check if all steps are the same size 43 41 SameSize = np.sum(np.abs(datasize - datasize[0])) == 0 … … 50 48 outdat = np.nan * np.ones(findim) 51 49 for step in range(rows): 50 #slicer is the data slice in the final array 52 51 slicer = [slice(0, d) for d in datasize[step, :]] 53 52 slicer = np.insert(slicer, 0, step) … … 55 54 outdat[tuple(slicer)] = np.reshape(self.data[startpoint:startpoint + curlen], newshape=(datasize[step, :])) 56 55 startpoint += curlen 56 #outmasked = ma.masked_array(outdat, mask=np.where(np.isnan(outdat), 1, 0)) 57 57 return outdat 58 #as much scalars as st pes (or less) so just one value per step58 #as much scalars as steps (or less) so just one value per step 59 59 else: 60 60 return np.squeeze(np.asarray(self.data)) -
issm/trunk-jpl/src/m/contrib/defleurian/netCDF/restable.m
r26871 r26901 15 15 %we save the size of the current step for further treatment 16 16 self.sizes=[self.sizes;fliplr(size(stepvar))]; 17 flatdata=reshape(stepvar, 1, []); 17 %we need to transpose to follow the indexing 18 flatdata=reshape(stepvar', 1, []); 18 19 self.data=[self.data,flatdata]; 19 20 end … … 32 33 findim=[maxsize, rows]; 33 34 %first check if all steps are the same size 34 SameSize = sum(self.sizes - self.sizes(1 ))==0;35 SameSize = sum(self.sizes - self.sizes(1, :))==0; 35 36 if SameSize, 36 37 %same size for all steps, just reshape … … 40 41 startpoint=1; 41 42 datadim=ndims(self.data); 42 outdat= NaN*ones(findim);43 outdat=nan(findim); 43 44 for r=1:rows 44 45 curlen = prod(self.sizes(r, :)); 45 outdat(1:self.sizes(r,1), 1:self.sizes(r,2), r) = reshape(self.data(startpoint:startpoint+curlen-1),self.sizes(r,:)); 46 outdat(1:self.sizes(r,1), 1:self.sizes(r,2), r) = reshape(self.data(startpoint:startpoint+ ... 47 curlen-1),self.sizes(r,:)); 46 48 startpoint = startpoint+curlen; 47 49 end -
issm/trunk-jpl/src/m/io/loadvars.py
r26872 r26901 9 9 from netCDF4 import Dataset, chartostring 10 10 import numpy as np 11 import numpy.ma as ma 11 12 from importlib import import_module 12 13 from model import * … … 31 32 filename = '' 32 33 nvdict = {} 33 debug = True # print messages if true34 verbose = 0 # 0 for silent 5 for chatty 34 35 35 36 if len(args) >= 1 and isinstance(args[0], str): … … 90 91 for mod in dict.keys(classtype): 91 92 #==== First we create the model structure {{{ 92 if debug:93 if verbose > 0: 93 94 print(' ==== Now treating classtype {}'.format(mod)) 94 95 if mod not in classtree.keys(): … … 97 98 # this points to a subclass (results.TransientSolution for example) 98 99 curclass = NCFile.groups[classtree[mod][0]].groups[classtree[mod][1]] 99 if debug:100 if verbose > 0: 100 101 print(" ==> {} is of class {}".format(mod, classtype[mod])) 101 102 if classtype[mod][0] == 'results.solutionstep': #Treating results {{{ … … 148 149 #}}} 149 150 else: 150 if debug:151 if verbose > 0: 151 152 print(" Using the default for md.{}.{}, is that right??".format(classtree[mod][0], classtree[mod][1])) 152 153 try: 153 154 modulename = split(r'\.', classtype[mod][0])[0] 154 if debug:155 if verbose > 0: 155 156 print(" trying to import {} from {}".format(classtype[mod][0], modulename)) 156 157 nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = getattr(classtype[mod][1], modulename)() … … 168 169 nvdict['md'].__dict__[mod] = getattr(classtype[mod][1], modulename)() 169 170 Tree = nvdict['md'].__dict__[classtree[mod][0]] 170 if debug:171 if verbose > 0: 171 172 print(" for {} Tree is a {} with len {}".format(mod, Tree.__class__.__name__, len(curclass.groups))) 172 173 # }}} … … 190 191 varval = listclass.variables[str(var)] 191 192 vardim = varval.ndim 192 if debug:193 if verbose > 0: 193 194 print(" ==> treating var {} of dimension {}".format(var, vardim)) 194 195 #There is a special treatment for results to account for its specific structure … … 226 227 timelist = np.arange(0, len(NCFile.dimensions['Time'])) 227 228 for t in timelist: 228 if debug:229 if verbose > 5: 229 230 print("filing step {} for {}".format(t, var)) 230 231 if vardim == 0: 231 232 Tree[t].__dict__[str(var)] = varval[:].data 232 233 elif vardim == 1: 233 Tree[t].__dict__[str(var)] = varval[t].data 234 stepval = ma.masked_array(varval[t].data, mask=np.where(np.isnan(varval[t]), 1, 0)) 235 Tree[t].__dict__[str(var)] = ma.compressed(stepval) 234 236 elif vardim == 2: 235 Tree[t].__dict__[str(var)] = varval[t, :].data 237 stepval = ma.masked_array(varval[t, :].data, mask=np.where(np.isnan(varval[t, :]), 1, 0)) 238 Tree[t].__dict__[str(var)] = ma.compressed(stepval) 236 239 elif vardim == 3: 237 Tree[t].__dict__[str(var)] = varval[t, :, :].data 240 stepval = ma.masked_array(varval[t, :, :].data, mask=np.where(np.isnan(varval[t, :, :]), 1, 0)) 241 Tree[t].__dict__[str(var)] = ma.compressed(stepval).reshape((stepval.count(0)[0], stepval.count(1)[0])) 238 242 else: 239 243 print('table dimension greater than 3 not implemented yet') 240 244 else: 241 if debug:245 if verbose > 0: 242 246 print("filing step {} for {}".format(groupindex, var)) 243 247 Tree[groupindex].__dict__[str(var)] = varval[:].data … … 254 258 255 259 elif vardim == 1: #that is a vector 256 if debug:260 if verbose > 0: 257 261 print(" for variable {} type is {}".format(str(var), varval.dtype)) 258 262 if varval.dtype == str: … … 279 283 elif vardim == 2: 280 284 #dealling with dict 281 if debug:285 if verbose > 0: 282 286 print(" for variable {} type is {}".format(str(var), varval.dtype)) 283 287 if varval.dtype == str: #that is for dictionaries … … 311 315 #==== And with atribute {{{ 312 316 for attr in listclass.ncattrs(): 313 if debug:317 if verbose > 0: 314 318 print(" ==> treating attribute {}".format(attr)) 315 319 if attr != 'classtype': #classtype is for treatment, don't get it back … … 319 323 attribute = attr 320 324 if type(Tree) == list: 321 if debug:325 if verbose > 0: 322 326 print(" printing with index 0") 323 327 if listtype == 'dict': -
issm/trunk-jpl/test/NightlyRun/runme.py
r26857 r26901 18 18 try: 19 19 from arch import archread 20 except: # ISSM_DIR is not on path20 except: # ISSM_DIR is not on path 21 21 import devpath 22 22 … … 133 133 test_ids = test_ids.difference(exclude_ids) 134 134 # }}} 135 if procedure == 'runFromNC': 136 #That is a bamg test 137 test_ids = test_ids.difference([119, 514]) 138 # that is smbGEMB format is weird for the test 139 test_ids = test_ids.difference([243, 244, 252, 253]) 140 #those are amr runs where the index is missing from fieldnames 141 test_ids = test_ids.difference([462, 463, 464, 465]) 142 #test247 solves for thermal and transient which makes it complex to check 143 test_ids = test_ids.difference([247]) 144 #test 902 is running two models with different stepping 145 test_ids = test_ids.difference([902]) 146 #I have a size issue in 517 needs investigation 147 test_ids = test_ids.difference([517]) 148 135 149 #Process Ids according to benchmarks {{{ 136 150 if benchmark == 'nightly': … … 200 214 if 'Solution' in key: 201 215 solvetype = split('Solution', key)[0] 216 217 #we save the results, scrap them and solve. 218 loaded_res = mdl.results 202 219 mdl.results = [] 203 220 mdl = solve(mdl, solvetype) 221 222 #we loop on the field_names from the nghtly test 204 223 for k, fieldname in enumerate(Tmod.field_names): 205 224 try: 225 #first look for indexing 206 226 if search(r'\d+$', fieldname): 207 227 index = int(search(r'\d+$', fieldname).group()) - 1 208 228 fieldname = fieldname[:search(r'\d+$', fieldname).start()] 209 el se:229 elif 'FirstStep' in fieldname: 210 230 index = 0 211 #Get field from nc run 231 fieldname = fieldname[:search('FirstStep', fieldname).start()] 232 elif 'SecondStep' in fieldname: 233 index = 1 234 fieldname = fieldname[:search('SecondStep', fieldname).start()] 235 elif 'ThirdStep' in fieldname: 236 index = 2 237 fieldname = fieldname[:search('ThirdStep', fieldname).start()] 238 else: 239 index = 0 240 241 #Then check if the key exists in the loaded results 242 try: 243 reskeys = mdl.results.__dict__[solvetype + 'Solution'][index].__dict__.keys() 244 except TypeError: 245 # most probably a steady state so no subscripting 246 reskeys = mdl.results.__dict__[solvetype + 'Solution'].__dict__.keys() 247 if fieldname not in reskeys: 248 sufixes = ["P1bubble", "P1bubbleCondensed", "LliboutryDuval", "CuffeyTemperate", "SSA", "HO", "FS", "P1xP", "P2xP", 249 'MINI', 'MINIcondensed', 'TaylorHood', 'XTaylorHood', 'LATaylorHood', 'CrouzeixRaviart', 'LACrouzeixRaviart'] 250 namedifs = {'Misfits': 'J', 251 'D': 'DamageDbar', 252 'F': 'DamageF', 253 'MaterialsRheologyB': 'MaterialsRheologyBbar', 254 'SedimentWaterHead': 'SedimentHead', 255 'EplWaterHead': 'EplHead', 256 'SedimentWaterHeadSubstep': 'SedimentHeadSubstep', 257 'EplWaterHeadSubstep': 'EplHeadSubstep', 258 'Volume': 'IceVolume', 259 'Bed': 'Base', 260 'SMB': 'SmbMassBalance'} 261 262 if fieldname in namedifs.keys(): 263 #Some fields are not consistent 264 fieldname = namedifs[fieldname] 265 elif any([suf in fieldname for suf in sufixes]): 266 #some test have loops that mess up with naming 267 try: 268 sufix = sufixes[np.squeeze(np.where([suf in fieldname for suf in sufixes]))] 269 except TypeError: 270 #probably severalmatches, we take the last one which should be the good one (Needs to be controled in the list above) 271 sufix = sufixes[np.squeeze(np.where([suf in fieldname for suf in sufixes]))[-1]] 272 fieldname = fieldname[:search(sufix, fieldname).start()] 273 elif fieldname.endswith("P") and index == 1: 274 #we are looking for P2 but 2 as been considered as an index and so shifted by -1 275 fieldname = fieldname[:-1] 276 else: 277 # could be that the index selected above is part of the name 278 fieldname = fieldname + str(index + 1) 212 279 try: 213 280 field = mdl.results.__dict__[solvetype + 'Solution'][index].__dict__[fieldname] 281 loaded_field = loaded_res.__dict__[solvetype + 'Solution'][index].__dict__[fieldname] 282 except TypeError: 283 # most probably a steady state so no subscripting 284 try: 285 field = mdl.results.__dict__[solvetype + 'Solution'].__dict__[fieldname] 286 loaded_field = loaded_res.__dict__[solvetype + 'Solution'].__dict__[fieldname] 287 except KeyError: 288 print("WARNING: {}{} does not exist and checking will be skipped".format(fieldname, index + 1)) 289 continue 214 290 except KeyError: 215 print("WARNING: {} does not exist and checking will be skipped".format(fieldname))291 print("WARNING: {}{} does not exist and checking will be skipped".format(fieldname, index + 1)) 216 292 continue 217 #Get reference from std run 293 218 294 ref = Tmod.field_values[k] 219 295 #Get tolerance 220 296 tolerance = Tmod.field_tolerances[k] 297 #compute differences for the results computed from the nc file 221 298 error_diff = np.amax(np.abs(ref - field), axis=0) / (np.amax(np.abs(ref), axis=0) + float_info.epsilon) 222 299 if not np.isscalar(error_diff): 223 300 error_diff = error_diff[0] 224 301 302 #compute the differences for the results of the nc file 303 load_diff = np.amax(np.abs(np.squeeze(ref) - loaded_field), axis=0) / (np.amax(np.abs(np.squeeze(ref)), axis=0) + float_info.epsilon) 304 if not np.isscalar(load_diff): 305 load_diff = load_diff[0] 306 225 307 #disp test result 226 if (np.any(error_diff > tolerance) or np.isnan(error_diff)): 227 print(('ERROR difference: {:7.2g} > {:7.2g} test id: {} test name: {} field: {}'.format(error_diff, tolerance, id, id_string, fieldname))) 308 if (np.any(error_diff > tolerance) or np.isnan(error_diff)) and (np.any(load_diff > tolerance) or np.isnan(load_diff)): 309 if abs(error_diff - load_diff) < tolerance: 310 print(('WARNING difference: {:7.2g} > {:7.2g} test id: {} field: {}{} differs from computation but equal to saved results'.format(error_diff, tolerance, id, fieldname, index + 1))) 311 else: 312 print(('ERROR difference: {:7.2g} > {:7.2g} test id: {} field: {}{} is false in both loaded and computed results'.format(error_diff, tolerance, id, fieldname, index + 1))) 313 errorcount += 1 314 erroredtest_list.append(id) 315 elif (np.any(error_diff > tolerance) or np.isnan(error_diff)): 316 print(('ERROR difference: {:7.2g} > {:7.2g} test id: {} test name: {} field: {}{}'.format(error_diff, tolerance, id, id_string, fieldname, index + 1))) 228 317 errorcount += 1 229 318 erroredtest_list.append(id) 230 else: 231 print(('SUCCESS difference: {:7.2g} < {:7.2g} test id: {} test name: {} field: {}'.format(error_diff, tolerance, id, id_string, fieldname))) 319 elif (np.any(load_diff > tolerance) or np.isnan(load_diff)): 320 print(('SAVEERROR difference: {:7.2g} > {:7.2g} test id: {} test name: {} saved result : {}{}'.format(load_diff, tolerance, id, id_string, fieldname, index + 1))) 321 errorcount += 1 322 erroredtest_list.append(id) 323 else: 324 print(('SUCCESS difference: {:7.2g} < {:7.2g} test id: {} test name: {} field: {}{}'.format(error_diff, tolerance, id, id_string, fieldname, index + 1))) 325 #disp only if errors for the results 232 326 233 327 except Exception as message:
Note:
See TracChangeset
for help on using the changeset viewer.
