Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 26900) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/export_netCDF.py (revision 26901) @@ -1,4 +1,5 @@ from netCDF4 import Dataset import numpy as np +import numpy.ma as ma import time import collections @@ -30,5 +31,4 @@ datasize = np.squeeze(self.sizes) maxsize = [] - try: dims = np.arange(np.shape(datasize)[1]) @@ -37,7 +37,5 @@ except IndexError: maxsize.append(np.nanmax(datasize[:])) - findim = np.insert(maxsize, 0, rows) - #first check if all steps are the same size SameSize = np.sum(np.abs(datasize - datasize[0])) == 0 @@ -50,4 +48,5 @@ outdat = np.nan * np.ones(findim) for step in range(rows): + #slicer is the data slice in the final array slicer = [slice(0, d) for d in datasize[step, :]] slicer = np.insert(slicer, 0, step) @@ -55,6 +54,7 @@ outdat[tuple(slicer)] = np.reshape(self.data[startpoint:startpoint + curlen], newshape=(datasize[step, :])) startpoint += curlen + #outmasked = ma.masked_array(outdat, mask=np.where(np.isnan(outdat), 1, 0)) return outdat - #as much scalars as stpes (or less) so just one value per step + #as much scalars as steps (or less) so just one value per step else: return np.squeeze(np.asarray(self.data)) Index: /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/restable.m =================================================================== --- /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/restable.m (revision 26900) +++ /issm/trunk-jpl/src/m/contrib/defleurian/netCDF/restable.m (revision 26901) @@ -15,5 +15,6 @@ %we save the size of the current step for further treatment self.sizes=[self.sizes;fliplr(size(stepvar))]; - flatdata=reshape(stepvar, 1, []); + %we need to transpose to follow the indexing + flatdata=reshape(stepvar', 1, []); self.data=[self.data,flatdata]; end @@ -32,5 +33,5 @@ findim=[maxsize, rows]; %first check if all steps are the same size - SameSize = sum(self.sizes - self.sizes(1))==0; + SameSize = sum(self.sizes - self.sizes(1, :))==0; if SameSize, %same size for all steps, just reshape @@ -40,8 +41,9 @@ startpoint=1; datadim=ndims(self.data); - outdat=NaN*ones(findim); + outdat=nan(findim); for r=1:rows curlen = prod(self.sizes(r, :)); - outdat(1:self.sizes(r,1), 1:self.sizes(r,2), r) = reshape(self.data(startpoint:startpoint+curlen-1),self.sizes(r,:)); + outdat(1:self.sizes(r,1), 1:self.sizes(r,2), r) = reshape(self.data(startpoint:startpoint+ ... + curlen-1),self.sizes(r,:)); startpoint = startpoint+curlen; end Index: /issm/trunk-jpl/src/m/io/loadvars.py =================================================================== --- /issm/trunk-jpl/src/m/io/loadvars.py (revision 26900) +++ /issm/trunk-jpl/src/m/io/loadvars.py (revision 26901) @@ -9,4 +9,5 @@ from netCDF4 import Dataset, chartostring import numpy as np +import numpy.ma as ma from importlib import import_module from model import * @@ -31,5 +32,5 @@ filename = '' nvdict = {} - debug = True # print messages if true + verbose = 0 # 0 for silent 5 for chatty if len(args) >= 1 and isinstance(args[0], str): @@ -90,5 +91,5 @@ for mod in dict.keys(classtype): #==== First we create the model structure {{{ - if debug: + if verbose > 0: print(' ==== Now treating classtype {}'.format(mod)) if mod not in classtree.keys(): @@ -97,5 +98,5 @@ # this points to a subclass (results.TransientSolution for example) curclass = NCFile.groups[classtree[mod][0]].groups[classtree[mod][1]] - if debug: + if verbose > 0: print(" ==> {} is of class {}".format(mod, classtype[mod])) if classtype[mod][0] == 'results.solutionstep': #Treating results {{{ @@ -148,9 +149,9 @@ #}}} else: - if debug: + if verbose > 0: print(" Using the default for md.{}.{}, is that right??".format(classtree[mod][0], classtree[mod][1])) try: modulename = split(r'\.', classtype[mod][0])[0] - if debug: + if verbose > 0: print(" trying to import {} from {}".format(classtype[mod][0], modulename)) nvdict['md'].__dict__[classtree[mod][0]].__dict__[classtree[mod][1]] = getattr(classtype[mod][1], modulename)() @@ -168,5 +169,5 @@ nvdict['md'].__dict__[mod] = getattr(classtype[mod][1], modulename)() Tree = nvdict['md'].__dict__[classtree[mod][0]] - if debug: + if verbose > 0: print(" for {} Tree is a {} with len {}".format(mod, Tree.__class__.__name__, len(curclass.groups))) # }}} @@ -190,5 +191,5 @@ varval = listclass.variables[str(var)] vardim = varval.ndim - if debug: + if verbose > 0: print(" ==> treating var {} of dimension {}".format(var, vardim)) #There is a special treatment for results to account for its specific structure @@ -226,18 +227,21 @@ timelist = np.arange(0, len(NCFile.dimensions['Time'])) for t in timelist: - if debug: + if verbose > 5: print("filing step {} for {}".format(t, var)) if vardim == 0: Tree[t].__dict__[str(var)] = varval[:].data elif vardim == 1: - Tree[t].__dict__[str(var)] = varval[t].data + stepval = ma.masked_array(varval[t].data, mask=np.where(np.isnan(varval[t]), 1, 0)) + Tree[t].__dict__[str(var)] = ma.compressed(stepval) elif vardim == 2: - Tree[t].__dict__[str(var)] = varval[t, :].data + stepval = ma.masked_array(varval[t, :].data, mask=np.where(np.isnan(varval[t, :]), 1, 0)) + Tree[t].__dict__[str(var)] = ma.compressed(stepval) elif vardim == 3: - Tree[t].__dict__[str(var)] = varval[t, :, :].data + stepval = ma.masked_array(varval[t, :, :].data, mask=np.where(np.isnan(varval[t, :, :]), 1, 0)) + Tree[t].__dict__[str(var)] = ma.compressed(stepval).reshape((stepval.count(0)[0], stepval.count(1)[0])) else: print('table dimension greater than 3 not implemented yet') else: - if debug: + if verbose > 0: print("filing step {} for {}".format(groupindex, var)) Tree[groupindex].__dict__[str(var)] = varval[:].data @@ -254,5 +258,5 @@ elif vardim == 1: #that is a vector - if debug: + if verbose > 0: print(" for variable {} type is {}".format(str(var), varval.dtype)) if varval.dtype == str: @@ -279,5 +283,5 @@ elif vardim == 2: #dealling with dict - if debug: + if verbose > 0: print(" for variable {} type is {}".format(str(var), varval.dtype)) if varval.dtype == str: #that is for dictionaries @@ -311,5 +315,5 @@ #==== And with atribute {{{ for attr in listclass.ncattrs(): - if debug: + if verbose > 0: print(" ==> treating attribute {}".format(attr)) if attr != 'classtype': #classtype is for treatment, don't get it back @@ -319,5 +323,5 @@ attribute = attr if type(Tree) == list: - if debug: + if verbose > 0: print(" printing with index 0") if listtype == 'dict': Index: /issm/trunk-jpl/test/NightlyRun/runme.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/runme.py (revision 26900) +++ /issm/trunk-jpl/test/NightlyRun/runme.py (revision 26901) @@ -18,5 +18,5 @@ try: from arch import archread -except: # ISSM_DIR is not on path +except: # ISSM_DIR is not on path import devpath @@ -133,4 +133,18 @@ test_ids = test_ids.difference(exclude_ids) # }}} + if procedure == 'runFromNC': + #That is a bamg test + test_ids = test_ids.difference([119, 514]) + # that is smbGEMB format is weird for the test + test_ids = test_ids.difference([243, 244, 252, 253]) + #those are amr runs where the index is missing from fieldnames + test_ids = test_ids.difference([462, 463, 464, 465]) + #test247 solves for thermal and transient which makes it complex to check + test_ids = test_ids.difference([247]) + #test 902 is running two models with different stepping + test_ids = test_ids.difference([902]) + #I have a size issue in 517 needs investigation + test_ids = test_ids.difference([517]) + #Process Ids according to benchmarks {{{ if benchmark == 'nightly': @@ -200,34 +214,114 @@ if 'Solution' in key: solvetype = split('Solution', key)[0] + + #we save the results, scrap them and solve. + loaded_res = mdl.results mdl.results = [] mdl = solve(mdl, solvetype) + + #we loop on the field_names from the nghtly test for k, fieldname in enumerate(Tmod.field_names): try: + #first look for indexing if search(r'\d+$', fieldname): index = int(search(r'\d+$', fieldname).group()) - 1 fieldname = fieldname[:search(r'\d+$', fieldname).start()] - else: + elif 'FirstStep' in fieldname: index = 0 - #Get field from nc run + fieldname = fieldname[:search('FirstStep', fieldname).start()] + elif 'SecondStep' in fieldname: + index = 1 + fieldname = fieldname[:search('SecondStep', fieldname).start()] + elif 'ThirdStep' in fieldname: + index = 2 + fieldname = fieldname[:search('ThirdStep', fieldname).start()] + else: + index = 0 + + #Then check if the key exists in the loaded results + try: + reskeys = mdl.results.__dict__[solvetype + 'Solution'][index].__dict__.keys() + except TypeError: + # most probably a steady state so no subscripting + reskeys = mdl.results.__dict__[solvetype + 'Solution'].__dict__.keys() + if fieldname not in reskeys: + sufixes = ["P1bubble", "P1bubbleCondensed", "LliboutryDuval", "CuffeyTemperate", "SSA", "HO", "FS", "P1xP", "P2xP", + 'MINI', 'MINIcondensed', 'TaylorHood', 'XTaylorHood', 'LATaylorHood', 'CrouzeixRaviart', 'LACrouzeixRaviart'] + namedifs = {'Misfits': 'J', + 'D': 'DamageDbar', + 'F': 'DamageF', + 'MaterialsRheologyB': 'MaterialsRheologyBbar', + 'SedimentWaterHead': 'SedimentHead', + 'EplWaterHead': 'EplHead', + 'SedimentWaterHeadSubstep': 'SedimentHeadSubstep', + 'EplWaterHeadSubstep': 'EplHeadSubstep', + 'Volume': 'IceVolume', + 'Bed': 'Base', + 'SMB': 'SmbMassBalance'} + + if fieldname in namedifs.keys(): + #Some fields are not consistent + fieldname = namedifs[fieldname] + elif any([suf in fieldname for suf in sufixes]): + #some test have loops that mess up with naming + try: + sufix = sufixes[np.squeeze(np.where([suf in fieldname for suf in sufixes]))] + except TypeError: + #probably severalmatches, we take the last one which should be the good one (Needs to be controled in the list above) + sufix = sufixes[np.squeeze(np.where([suf in fieldname for suf in sufixes]))[-1]] + fieldname = fieldname[:search(sufix, fieldname).start()] + elif fieldname.endswith("P") and index == 1: + #we are looking for P2 but 2 as been considered as an index and so shifted by -1 + fieldname = fieldname[:-1] + else: + # could be that the index selected above is part of the name + fieldname = fieldname + str(index + 1) try: field = mdl.results.__dict__[solvetype + 'Solution'][index].__dict__[fieldname] + loaded_field = loaded_res.__dict__[solvetype + 'Solution'][index].__dict__[fieldname] + except TypeError: + # most probably a steady state so no subscripting + try: + field = mdl.results.__dict__[solvetype + 'Solution'].__dict__[fieldname] + loaded_field = loaded_res.__dict__[solvetype + 'Solution'].__dict__[fieldname] + except KeyError: + print("WARNING: {}{} does not exist and checking will be skipped".format(fieldname, index + 1)) + continue except KeyError: - print("WARNING: {} does not exist and checking will be skipped".format(fieldname)) + print("WARNING: {}{} does not exist and checking will be skipped".format(fieldname, index + 1)) continue - #Get reference from std run + ref = Tmod.field_values[k] #Get tolerance tolerance = Tmod.field_tolerances[k] + #compute differences for the results computed from the nc file error_diff = np.amax(np.abs(ref - field), axis=0) / (np.amax(np.abs(ref), axis=0) + float_info.epsilon) if not np.isscalar(error_diff): error_diff = error_diff[0] + #compute the differences for the results of the nc file + 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) + if not np.isscalar(load_diff): + load_diff = load_diff[0] + #disp test result - if (np.any(error_diff > tolerance) or np.isnan(error_diff)): - print(('ERROR difference: {:7.2g} > {:7.2g} test id: {} test name: {} field: {}'.format(error_diff, tolerance, id, id_string, fieldname))) + if (np.any(error_diff > tolerance) or np.isnan(error_diff)) and (np.any(load_diff > tolerance) or np.isnan(load_diff)): + if abs(error_diff - load_diff) < tolerance: + 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))) + else: + 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))) + errorcount += 1 + erroredtest_list.append(id) + elif (np.any(error_diff > tolerance) or np.isnan(error_diff)): + print(('ERROR difference: {:7.2g} > {:7.2g} test id: {} test name: {} field: {}{}'.format(error_diff, tolerance, id, id_string, fieldname, index + 1))) errorcount += 1 erroredtest_list.append(id) - else: - print(('SUCCESS difference: {:7.2g} < {:7.2g} test id: {} test name: {} field: {}'.format(error_diff, tolerance, id, id_string, fieldname))) + elif (np.any(load_diff > tolerance) or np.isnan(load_diff)): + print(('SAVEERROR difference: {:7.2g} > {:7.2g} test id: {} test name: {} saved result : {}{}'.format(load_diff, tolerance, id, id_string, fieldname, index + 1))) + errorcount += 1 + erroredtest_list.append(id) + else: + print(('SUCCESS difference: {:7.2g} < {:7.2g} test id: {} test name: {} field: {}{}'.format(error_diff, tolerance, id, id_string, fieldname, index + 1))) + #disp only if errors for the results except Exception as message: