Index: /issm/trunk-jpl/src/m/modules/Chaco.py =================================================================== --- /issm/trunk-jpl/src/m/modules/Chaco.py (revision 23230) +++ /issm/trunk-jpl/src/m/modules/Chaco.py (revision 23231) @@ -5,5 +5,5 @@ Usage: - [assgn] = Chaco(A,vwgts,ewgts,x,y,z,options,nparts,goal); + assgn = Chaco(A,vwgts,ewgts,x,y,z,options,nparts,goal); A: Input adjacency matrix Index: /issm/trunk-jpl/src/m/modules/Scotch.py =================================================================== --- /issm/trunk-jpl/src/m/modules/Scotch.py (revision 23230) +++ /issm/trunk-jpl/src/m/modules/Scotch.py (revision 23231) @@ -2,5 +2,5 @@ def Scotch(*varargin): -'''SCOTCH - Scotch partitioner + '''SCOTCH - Scotch partitioner Usage: Index: /issm/trunk-jpl/src/m/partition/partitioner.py =================================================================== --- /issm/trunk-jpl/src/m/partition/partitioner.py (revision 23230) +++ /issm/trunk-jpl/src/m/partition/partitioner.py (revision 23231) @@ -3,8 +3,9 @@ import MatlabFuncs as m from adjacency import * -#from Chaco import * +from Chaco import * #from Scotch import * -#from MeshPartition import * +from MeshPartition import * from project3d import * +from mesh2d import * def partitioner(md,*varargin): @@ -60,27 +61,28 @@ if m.strcmpi(package,'chaco'): - raise RuntimeError('Chaco is not currently supported for this function') + #raise RuntimeError('Chaco is not currently supported for this function') # default method (from chaco.m) - #method=np.array([1,1,0,0,1,1,50,0,.001,7654321]).reshape(-1,1) - #method[0]=3 # global method (3=inertial (geometric)) - #method[2]=0 # vertex weights (0=off, 1=on) + method=np.array([1,1,0,0,1,1,50,0,.001,7654321]) + method[0]=3 # global method (3=inertial (geometric)) + method[2]=0 # vertex weights (0=off, 1=on) #specify bisection - #method[5]=options.getfieldvalue('section')# ndims (1=bisection, 2=quadrisection, 3=octasection) + method[5]=options.getfieldvalue('section')# ndims (1=bisection, 2=quadrisection, 3=octasection) #are we using weights? - #if m.strcmpi(options.getfieldvalue('weighting'),'on'): - #weights=np.floor(md.qmu.vertex_weight/min(md.qmu.vertex_weight)) - #method[2]=1 - #else: - #weights=[] + if m.strcmpi(options.getfieldvalue('weighting'),'on'): + weights=np.floor(md.qmu.vertex_weight/min(md.qmu.vertex_weight)) + method[2]=1 + else: + weights=[] + method = method.reshape(-1,1) # transpose to 1x10 instead of 10 + # partition into nparts - #if isinstance(md.mesh,mesh2d): - #part=Chaco(md.qmu.adjacency,weights,[],md.mesh.x, md.mesh.y,np.zeros((md.mesh.numberofvertices,)),method,npart,[]).T+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,[]).T+1 #index partitions from 1 up. like metis. - + if isinstance(md.mesh,mesh2d): + part=np.array(Chaco(md.qmu.adjacency,weights,np.array([]),md.mesh.x, md.mesh.y,np.zeros((md.mesh.numberofvertices,)),method,npart,np.array([]))).T+1 #index partitions from 1 up. like metis. + else: + part=np.array(Chaco(md.qmu.adjacency,weights,np.array([]),md.mesh.x, md.mesh.y,md.mesh.z,method,npart,np.array([]))).T+1 #index partitions from 1 up. like metis. elif m.strcmpi(package,'scotch'): Index: /issm/trunk-jpl/src/m/qmu/dakota_in_write.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_in_write.py (revision 23230) +++ /issm/trunk-jpl/src/m/qmu/dakota_in_write.py (revision 23231) @@ -62,5 +62,5 @@ filei2=fullfile(pathstr,name+ext) - print 'Opening Dakota input file \''+filei2 + '\'' + print 'Opening Dakota input file \''+filei2 + '\'.' try: with open(filei2,'w+') as fidi: Index: /issm/trunk-jpl/src/m/qmu/dakota_out_parse.py =================================================================== --- /issm/trunk-jpl/src/m/qmu/dakota_out_parse.py (revision 23230) +++ /issm/trunk-jpl/src/m/qmu/dakota_out_parse.py (revision 23231) @@ -241,4 +241,15 @@ dmax95 =prctile_issm(data,95,0) + # Note: the following line may cause the following warning + # (should not crash or invalidate results) when one of + # the inputs does not change with respect to the + # other/s causing an internal divide-by-zero error + + #/usr/local/lib/python2.7/dist-packages/numpy/lib/function_base.py:3163: + # RuntimeWarning: invalid value encountered in true_divide + # c /= stddev[:,None] + + # (and/or the same but with "c /= stddev[None, :]") + dcorrel=np.corrcoef(data.T) Index: /issm/trunk-jpl/src/m/qmu/process_qmu_response_data.m =================================================================== --- /issm/trunk-jpl/src/m/qmu/process_qmu_response_data.m (revision 23230) +++ /issm/trunk-jpl/src/m/qmu/process_qmu_response_data.m (revision 23231) @@ -43,8 +43,6 @@ %ok, process the domains named in qmu_mass_flux_profiles, to build a list of segments (MatArray) md.qmu.mass_flux_segments=cell(num_mass_flux,1); - md.qmu.mass_flux_segments for i=1:num_mass_flux, md.qmu.mass_flux_segments{i}=MeshProfileIntersection(md.mesh.elements,md.mesh.x,md.mesh.y,[md.qmu.mass_flux_profile_directory '/' md.qmu.mass_flux_profiles{i}]); end - md.qmu.mass_flux_segments end Index: /issm/trunk-jpl/src/wrappers/Chaco/Chaco.cpp =================================================================== --- /issm/trunk-jpl/src/wrappers/Chaco/Chaco.cpp (revision 23230) +++ /issm/trunk-jpl/src/wrappers/Chaco/Chaco.cpp (revision 23231) @@ -47,13 +47,13 @@ /*Fetch Data*/ FetchChacoData(&nvtxs,&adjacency,&start,&ewgts,A_IN,EWGTS_IN); - FetchData(&vwgts,&nterms,VWGTS_IN); - FetchData(&x,&nterms,X_IN); - FetchData(&y,&nterms,Y_IN); - FetchData(&z,&nterms,Z_IN); - FetchData(&in_options,&nterms,OPTNS_IN); + FetchData(&vwgts,&nterms,VWGTS_IN); + FetchData(&x,&nterms,X_IN); + FetchData(&y,&nterms,Y_IN); + FetchData(&z,&nterms,Z_IN); + FetchData(&in_options,&nterms,OPTNS_IN); for (i=0;i<(nterms<10?nterms:10);i++) options[i]=in_options[i]; //copy in_options into default options - FetchData(&npart,NPARTS_IN); + FetchData(&npart,NPARTS_IN); //int * nparts=xNew(1); nparts[0]=npart; //weird Chacox interface ain't it? - FetchData(&goal,&nterms,GOAL_IN); + FetchData(&goal,&nterms,GOAL_IN); /*Allocate output: */ Index: /issm/trunk-jpl/src/wrappers/python/io/FetchPythonData.cpp =================================================================== --- /issm/trunk-jpl/src/wrappers/python/io/FetchPythonData.cpp (revision 23230) +++ /issm/trunk-jpl/src/wrappers/python/io/FetchPythonData.cpp (revision 23231) @@ -32,5 +32,5 @@ } else if (PyInt_Check(py_float)) - dscalar=(double)PyInt_AsLong(py_float); + dscalar=(double)PyInt_AsLong(py_float); else if (PyBool_Check(py_float)) dscalar=(double)PyLong_AsLong(py_float); @@ -132,5 +132,6 @@ double* dmatrix=NULL; double* matrix=NULL; - int M,N; + int M=0; + int N=0; int ndim; npy_intp* dims=NULL; @@ -236,5 +237,6 @@ /*output: */ int* matrix=NULL; - int M,N; + int M=0; + int N=0; int ndim; npy_intp* dims=NULL; @@ -325,5 +327,6 @@ bool* bmatrix=NULL; bool* matrix=NULL; - int M,N; + int M=0; + int N=0; int ndim; npy_intp* dims=NULL; @@ -411,5 +414,5 @@ double* dvector=NULL; double* vector=NULL; - int M; + int M=0; int ndim; npy_intp* dims=NULL; @@ -499,5 +502,5 @@ /*output: */ float* vector=NULL; - int M; + int M=0; int ndim; npy_intp* dims=NULL; @@ -539,5 +542,5 @@ /*transform into int vector: */ vector=xNew(M); - for(i=0;i(M); - for(i=0;i(nvtxs+1); + pyGetJc(a,nvtxs,mwstart); + + int* mwadjacency = xNewZeroInit(nzmax); + pyGetIr(a,nvtxs,nzmax,mwadjacency); + + int* start = xNew(nvtxs+1); + for(int i=0;i(nzmax); + for(int i=0; i(nedges); + for(int i = 0; i2) or int(dver[0])>4): + md.qmu.method = dmeth_params_set(md.qmu.method,'rng','rnum2') + +#parameters +md.qmu.params.direct = True +md.qmu.params.analysis_components = '' +md.qmu.params.interval_type = 'forward' +md.qmu.params.tabular_graphics_data = True +md.qmu.isdakota = 1 + +if version >= 6: + md.qmu.params.analysis_driver = 'matlab' + md.qmu.params.evaluation_scheduling = 'master' + md.qmu.params.processors_per_evaluation = 2 +else: + md.qmu.params.analysis_driver = 'stressbalance' + md.qmu.params.evaluation_concurrency = 1 + +md.stressbalance.reltol = 10**-5 #tighten for qmu analyses +md.transient.requested_outputs = ['IceVolume'] + +#solve +md.verbose = verbose('000000000') # this line is recommended +md = solve(md,'Transient','overwrite','y') +md.qmu.results = md.results.dakota + +#Fields and tolerances to track changes +md.results.dakota.moments = [] +for i in range(8): + md.results.dakota.moments.append(md.results.dakota.dresp_out[i].mean) + +for i in range(8): + md.results.dakota.moments.append(md.results.dakota.dresp_out[i].stddev) + +field_names = ['moments'] +field_tolerances = [1e-11] +field_values = [md.results.dakota.moments] Index: /issm/trunk-jpl/test/NightlyRun/test235.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test235.py (revision 23231) +++ /issm/trunk-jpl/test/NightlyRun/test235.py (revision 23231) @@ -0,0 +1,97 @@ +#Test Name: SquareShelfTranForceNeg2dDakotaLocal +import numpy as np +import scipy.io as spio +from os import getcwd +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from SMBgemb import * +from partitioner import * +from dmeth_params_set import * + +md = triangle(model(),'../Exp/Square.exp',180000.) +md = setmask(md,'all','') +md = parameterize(md,'../Par/SquareShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) + +md.timestepping.time_step = 1 +md.settings.output_frequency = 1 +md.timestepping.final_time = 4 + +smb = np.ones((md.mesh.numberofvertices,))*3.6 +smb = np.array([smb,smb*-1]).T + +md.smb.mass_balance = smb +md.smb.mass_balance = np.concatenate((md.smb.mass_balance,[[1.5,3]])) +md.transient.isthermal = 0 + +#Dakota options + +#dakota version +version = IssmConfig('_DAKOTA_VERSION_') +version = float(version[0]) + +#partitioning +md.qmu.numberofpartitions = 20 +md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on') +md.qmu.partition = md.qmu.partition - 1 + +#variables +md.qmu.variables.surface_mass_balance = normal_uncertain.normal_uncertain('scaled_SmbMassBalance',1,100) + +#responses +md.qmu.responses.MaxVel = response_function.response_function('MaxVel',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.IceVolume = response_function.response_function('IceVolume',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux1 = response_function.response_function('indexed_MassFlux_1',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux2 = response_function.response_function('indexed_MassFlux_2',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux3 = response_function.response_function('indexed_MassFlux_3',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux4 = response_function.response_function('indexed_MassFlux_4',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux5 = response_function.response_function('indexed_MassFlux_5',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.massFlux6 = response_function.response_function('indexed_MassFlux_6',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) + +#mass flux profiles +md.qmu.mass_flux_profiles = ['../Exp/MassFlux1.exp','../Exp/MassFlux2.exp','../Exp/MassFlux3.exp','../Exp/MassFlux4.exp','../Exp/MassFlux5.exp','../Exp/MassFlux6.exp'] +md.qmu.mass_flux_profile_directory = getcwd() + +#method +md.qmu.method = dakota_method.dakota_method('nond_l') + +#parameters +md.qmu.params.direct = True +md.qmu.params.analysis_components = '' +md.qmu.params.interval_type = 'forward' +md.qmu.params.fd_gradient_step_size = '0.1' +md.qmu.isdakota = 1 + +if version >= 6: + md.qmu.params.analysis_driver = 'matlab' + md.qmu.params.evaluation_scheduling = 'master' + md.qmu.params.processors_per_evaluation = 2 +else: + md.qmu.params.analysis_driver = 'stressbalance' + md.qmu.params.evaluation_concurrency = 1 + +md.stressbalance.reltol = 10**-5 #tighten for qmu analyses +md.transient.requested_outputs = ['IceVolume'] + +#solve +md.verbose = verbose('000000000') # this line is recommended +md = solve(md,'Transient','overwrite','y') +md.qmu.results = md.results.dakota + +#Fields and tolerances to track changes +md.results.dakota.moments = [] +for i in range(8): + md.results.dakota.moments.append(md.results.dakota.dresp_out[i].mean) + +for i in range(8): + md.results.dakota.moments.append(md.results.dakota.dresp_out[i].stddev) + +field_names = ['moments'] +field_tolerances = [1e-11] +field_values = [md.results.dakota.moments] Index: /issm/trunk-jpl/test/NightlyRun/test413.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test413.py (revision 23231) +++ /issm/trunk-jpl/test/NightlyRun/test413.py (revision 23231) @@ -0,0 +1,66 @@ +#Test Name: SquareSheetShelfDiadSSA3dDakotaPart +import numpy as np +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from partitioner import * +from importancefactors import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) + +#Dakota options + +#dakota version +version = IssmConfig('_DAKOTA_VERSION_') +version = float(version[0]) + +#partitioning +md.qmu.numberofpartitions = 20 +md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on') +md.qmu.partition = md.qmu.partition-1 + +#variables +md.qmu.variables.rho_ice = normal_uncertain.normal_uncertain('MaterialsRhoIce',md.materials.rho_ice,0.01) +md.qmu.variables.drag_coefficient = normal_uncertain.normal_uncertain('scaled_FrictionCoefficient',1,0.01) + +#responses +md.qmu.responses.MaxVel = response_function.response_function('MaxVel',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) + +#method +md.qmu.method = dakota_method.dakota_method('nond_l') + +#parameters +md.qmu.params.direct = True +md.qmu.params.interval_type = 'forward' + +if version >= 6: + md.qmu.params.analysis_driver = 'matlab' + md.qmu.params.evaluation_scheduling = 'master' + md.qmu.params.processors_per_evaluation = 2 +else: + md.qmu.params.analysis_driver = 'stressbalance' + md.qmu.params.evaluation_concurrency = 1 + + +#imperative! +md.stressbalance.reltol = 10**-5 #tighten for qmu analyses +md.qmu.isdakota = 1 + +#solve +md.verbose = verbose('000000000') # this line is recommended +md = solve(md,'Stressbalance','overwrite','y') + +#Fields and tolerances to track changes +md.qmu.results = md.results.dakota +md.results.dakota.importancefactors = importancefactors(md,'scaled_FrictionCoefficient','MaxVel').T +field_names = ['importancefactors'] +field_tolerances = [1e-10] +field_values = [md.results.dakota.importancefactors] Index: /issm/trunk-jpl/test/NightlyRun/test414.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test414.py (revision 23231) +++ /issm/trunk-jpl/test/NightlyRun/test414.py (revision 23231) @@ -0,0 +1,92 @@ +#Test Name: SquareSheetShelfDiadSSA3dDakotaMassFlux +import numpy as np +from os import getcwd +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from partitioner import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) +md.materials.rho_ice = 10**7 #involved in the mass flux, make it easy +md.geometry.thickness[:] = 1 #make it easy +md.geometry.surface = md.geometry.base+md.geometry.thickness + +#constrain all velocities to 1 m/yr, in the y-direction +md.stressbalance.spcvx = np.zeros((md.mesh.numberofvertices,)) +md.stressbalance.spcvy = np.ones((md.mesh.numberofvertices,)) +md.stressbalance.spcvz = np.zeros((md.mesh.numberofvertices,)) + +#Dakota options + +#dakota version +version = IssmConfig('_DAKOTA_VERSION_') +version = float(version[0]) + +#partitioning +md.qmu.numberofpartitions = 20 +md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on') +md.qmu.partition = md.qmu.partition-1 + +#variables +md.qmu.variables.drag_coefficient = normal_uncertain.normal_uncertain('scaled_FrictionCoefficient',1,0.01) + +#responses +md.qmu.responses.MaxVel = response_function.response_function('MaxVel',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux1 = response_function.response_function('indexed_MassFlux_1',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux2 = response_function.response_function('indexed_MassFlux_2',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux3 = response_function.response_function('indexed_MassFlux_3',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux4 = response_function.response_function('indexed_MassFlux_4',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux5 = response_function.response_function('indexed_MassFlux_5',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux6 = response_function.response_function('indexed_MassFlux_6',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux7 = response_function.response_function('indexed_MassFlux_7',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) + +#mass flux profiles +md.qmu.mass_flux_profiles = ['../Exp/MassFlux1.exp','../Exp/MassFlux2.exp','../Exp/MassFlux3.exp','../Exp/MassFlux4.exp','../Exp/MassFlux5.exp','../Exp/MassFlux6.exp','../Exp/Square.exp'] +md.qmu.mass_flux_profile_directory = getcwd() + +#method +md.qmu.method = dakota_method.dakota_method('nond_l') + +#parameters +md.qmu.params.direct = True +md.qmu.params.interval_type = 'forward' +md.qmu.isdakota = 1 +md.stressbalance.reltol = 10**-5 #tighten for qmu analyses + +if version >= 6: + md.qmu.params.analysis_driver = 'matlab' + md.qmu.params.evaluation_scheduling = 'master' + md.qmu.params.processors_per_evaluation = 2 +else: + md.qmu.params.analysis_driver = 'stressbalance' + md.qmu.params.evaluation_concurrency = 1 + +#solve +md.verbose = verbose('000000000') # this line is recommended +md = solve(md,'Stressbalance','overwrite','y') +md.qmu.results = md.results.dakota + +#Fields and tolerances to track changes +#ok, mass flux of 3 profiles should be -3 Gt/yr -3 Gt/yr and the sum, which is -6 Gt/yr +#we recover those mass fluxes through the mean of the response. +#also, we recover the max velo, which should be 1m/yr. +#we put all that data in the moments, which we will use to test for success. +#also, check that the stddev are 0. +md.results.dakota.moments = [] +for i in range(8): + md.results.dakota.moments.append(md.results.dakota.dresp_out[i].mean) + +for i in range(8): + md.results.dakota.moments.append(md.results.dakota.dresp_out[i].stddev) + +field_names = ['moments'] +field_tolerances = [1e-11] +field_values = [md.results.dakota.moments] Index: /issm/trunk-jpl/test/NightlyRun/test417.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test417.py (revision 23231) +++ /issm/trunk-jpl/test/NightlyRun/test417.py (revision 23231) @@ -0,0 +1,106 @@ +#Test Name: SquareSheetShelfDiadSSA3dDakotaSamp +import numpy as np +from os import getcwd +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from partitioner import * +from dmeth_params_set import * + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',gethostname(),'np',3) +md.materials.rho_ice = 10**7 #involved in the mass flux, make it easy +md.geometry.thickness[:] = 1 #make it easy +md.geometry.surface = md.geometry.base+md.geometry.thickness + +#constrain all velocities to 1 m/yr, in the y-direction +md.stressbalance.spcvx[:] = 0 +md.stressbalance.spcvy[:] = 1 +md.stressbalance.spcvz[:] = 0 + +#Dakota options + +#dakota version +version = IssmConfig('_DAKOTA_VERSION_') +version = float(version[0]) + +#partitioning +md.qmu.numberofpartitions = 20 +md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on') +md.qmu.partition = md.qmu.partition-1 + +#variables +md.qmu.variables.drag_coefficient = normal_uncertain.normal_uncertain('scaled_FrictionCoefficient',1,0.01) + +#responses +md.qmu.responses.MaxVel = response_function.response_function('MaxVel',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux1 = response_function.response_function('indexed_MassFlux_1',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux2 = response_function.response_function('indexed_MassFlux_2',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux3 = response_function.response_function('indexed_MassFlux_3',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux4 = response_function.response_function('indexed_MassFlux_4',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux5 = response_function.response_function('indexed_MassFlux_5',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux6 = response_function.response_function('indexed_MassFlux_6',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.MassFlux7 = response_function.response_function('indexed_MassFlux_7',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) + +#mass flux profiles +md.qmu.mass_flux_profiles = ['../Exp/MassFlux1.exp','../Exp/MassFlux2.exp','../Exp/MassFlux3.exp','../Exp/MassFlux4.exp','../Exp/MassFlux5.exp','../Exp/MassFlux6.exp','../Exp/Square.exp'] +md.qmu.mass_flux_profile_directory = getcwd() + +#method +md.qmu.method = dakota_method.dakota_method('nond_samp') +md.qmu.method = dmeth_params_set(md.qmu.method,'seed',1234,'samples',20,'sample_type','lhs') + +#parameters +md.qmu.params.direct = True +md.qmu.params.analysis_components = '' +md.qmu.params.tabular_graphics_data = True; + +if version >= 6: + md.qmu.params.analysis_driver = 'matlab' + md.qmu.params.evaluation_scheduling = 'master' + md.qmu.params.processors_per_evaluation = 2 +else: + md.qmu.params.analysis_driver = 'stressbalance' + md.qmu.params.evaluation_concurrency = 1 + +#partitioning +md.qmu.numberofpartitions = 20 +md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on') +md.qmu.partition = md.qmu.partition-1 +md.qmu.isdakota = 1 + +md.stressbalance.reltol = 10**-5 #tighten for qmu analyses + +#solve +md.verbose = verbose('000000000') # this line is recommended + +# There may be a pair of numpy warnings in the function true_divide, +# this is normal and will not affect the results +# See src/m/qmu/dakota_out_parse.py, function "dak_tab_out" for details +md = solve(md,'Stressbalance','overwrite','y') + +#Fields and tolerances to track changes +md.qmu.results = md.results.dakota; + +#ok, mass flux of 3 profiles should be -3 Gt/yr -3 Gt/yr and the sum, which is -6 Gt/yr +#we recover those mass fluxes through the mean of the response. +#also, we recover the max velo, which should be 1m/yr. +#we put all that data in the montecarlo field, which we will use to test for success. +#also, check that the stddev are 0. +md.results.dakota.montecarlo = [] +for i in range(8): + md.results.dakota.montecarlo.append(md.results.dakota.dresp_out[i].mean) + +for i in range(8): + md.results.dakota.montecarlo.append(md.results.dakota.dresp_out[i].stddev) + +field_names = ['moments'] +field_tolerances = [1e-11] +field_values = [md.results.dakota.montecarlo] Index: /issm/trunk-jpl/test/NightlyRun/test444.py =================================================================== --- /issm/trunk-jpl/test/NightlyRun/test444.py (revision 23231) +++ /issm/trunk-jpl/test/NightlyRun/test444.py (revision 23231) @@ -0,0 +1,127 @@ +#Test Name: SquareShelfTranForceNeg2dDakotaLocal +import numpy as np +from os import getcwd +from model import * +from socket import gethostname +from triangle import * +from setmask import * +from parameterize import * +from setflowequation import * +from solve import * +from SMBgemb import * +from partitioner import * +from dmeth_params_set import * +from ContourToMesh import * +from regionaloutput import * + +#model not consistent: equality thickness=surface-base violated + +md = triangle(model(),'../Exp/Square.exp',150000.) +md = setmask(md,'../Exp/SquareShelf.exp','') +md = parameterize(md,'../Par/SquareSheetShelf.py') +md.geometry.bed = md.geometry.base.copy() +pos = np.where(md.mask.groundedice_levelset < 0) +md.geometry.bed[pos] = md.geometry.base[pos]-10 +md.friction.coefficient = 20. * np.ones((md.mesh.numberofvertices,)) +md.friction.p = np.ones((md.mesh.numberofelements,)) +md.friction.q = np.ones((md.mesh.numberofelements,)) +md.transient.isthermal = 0 +md.transient.isgroundingline = 1 +md.groundingline.migration = 'AggressiveMigration' + +md.settings.output_frequency = 3 +md = setflowequation(md,'SSA','all') +md.cluster = generic('name',oshostname(),'np',3) + +regionalmask = np.zeros((md.mesh.numberofvertices,)) +c_in = ContourToMesh(md.mesh.elements,md.mesh.x,md.mesh.y,'../Exp/SquareHalfRight.exp','node',1) +regionalmask[np.where(c_in)] = 1 +md.transient.requested_outputs = ['default','GroundedArea','FloatingArea','IceVolumeAboveFloatation','GroundedArea1','FloatingArea1','TotalFloatingBmb1','TotalGroundedBmb1','TotalSmb1', + 'IceMass1','IceVolume1','IceVolumeAboveFloatation1','IceVolumeAboveFloatation'] +md.outputdefinition.definitions.append(regionaloutput('name','GroundedArea1','outputnamestring','GroundedArea','mask',regionalmask, + 'definitionstring','Outputdefinition1')) +md.outputdefinition.definitions.append(regionaloutput('name','FloatingArea1','outputnamestring','FloatingArea','mask',regionalmask, + 'definitionstring','Outputdefinition2')) +md.outputdefinition.definitions.append(regionaloutput('name','TotalFloatingBmb1','outputnamestring','TotalFloatingBmb','mask',regionalmask, + 'definitionstring','Outputdefinition3')) +md.outputdefinition.definitions.append(regionaloutput('name','TotalGroundedBmb1','outputnamestring','TotalGroundedBmb','mask',regionalmask, + 'definitionstring','Outputdefinition4')) +md.outputdefinition.definitions.append(regionaloutput('name','IceMass1','outputnamestring','IceMass','mask',regionalmask, + 'definitionstring','Outputdefinition5')) +md.outputdefinition.definitions.append(regionaloutput('name','IceVolume1','outputnamestring','IceVolume','mask',regionalmask, + 'definitionstring','Outputdefinition6')) +md.outputdefinition.definitions.append(regionaloutput('name','IceVolumeAboveFloatation1','outputnamestring','IceVolumeAboveFloatation','mask',regionalmask, + 'definitionstring','Outputdefinition7')) +md.outputdefinition.definitions.append(regionaloutput('name','TotalSmb1','outputnamestring','TotalSmb','mask',regionalmask, + 'definitionstring','Outputdefinition8')) +md.outputdefinition.definitions.append(regionaloutput('name','TotalSmb2','outputnamestring','TotalSmb','mask',regionalmask, + 'definitionstring','Outputdefinition9')) + +md.extrude(3,1.) +md.collapse() + +#Dakota options + +#dakota version +version = IssmConfig('_DAKOTA_VERSION_') +version = float(version[0]) + +#variables +md.qmu.variables.drag_coefficient = normal_uncertain.normal_uncertain('scaled_BasalforcingsFloatingiceMeltingRate',1,0.1) + +#responses +md.qmu.responses.IceMass1 = response_function.response_function('Outputdefinition5',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.IceVolume1 = response_function.response_function('Outputdefinition6',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.IceVolumeAboveFloatation1 = response_function.response_function('Outputdefinition7',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.IceVolumeAboveFloatation = response_function.response_function('IceVolumeAboveFloatation',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.GroundedArea1 = response_function.response_function('Outputdefinition1',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.FloatingArea1 = response_function.response_function('Outputdefinition2',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.TotalFloatingBmb1 = response_function.response_function('Outputdefinition3',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.TotalGroundedBmb1 = response_function.response_function('Outputdefinition4',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.TotalSmb1 = response_function.response_function('Outputdefinition8',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.TotalSmb2 = response_function.response_function('Outputdefinition9',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) +md.qmu.responses.FloatingArea = response_function.response_function('FloatingArea',[],[0.0001,0.001,0.01,0.25,0.5,0.75,0.99,0.999,0.9999]) + +#method +md.qmu.method = dakota_method.dakota_method('nond_samp') +md.qmu.method = dmeth_params_set(md.qmu.method,'seed',1234,'samples',20,'sample_type','random') + +#parameters +md.qmu.params.direct = True +md.qmu.params.analysis_components = '' +md.qmu.params.tabular_graphics_data = True + +if version >= 6: + md.qmu.params.analysis_driver = 'matlab' + md.qmu.params.evaluation_scheduling = 'master' + md.qmu.params.processors_per_evaluation = 2 +else: + md.qmu.params.analysis_driver = 'stressbalance' + md.qmu.params.evaluation_concurrency = 1 + +#partitioning +md.qmu.numberofpartitions = 10 +md = partitioner(md,'package','chaco','npart',md.qmu.numberofpartitions,'weighting','on') +md.qmu.partition = md.qmu.partition-1 +md.qmu.isdakota = 1 + +md.stressbalance.reltol = 10**-5 #tighten for qmu analyses + +#solve +md.verbose = verbose('000000000') # this line is recommended +md = solve(md,'Transient','overwrite','y') + +#Fields and tolerances to track changes +md.qmu.results = md.results.dakota + +#we put all the mean and stdev data in the montecarlo field, which we will use to test for success. +md.results.dakota.montecarlo = [] +for i in range(11): + md.results.dakota.montecarlo.append(md.results.dakota.dresp_out[i].mean) + +for i in range(11): + md.results.dakota.montecarlo.append(md.results.dakota.dresp_out[i].stddev) + +field_names = ['montecarlo'] +field_tolerances = [1e-11] +field_values = [md.results.dakota.montecarlo]