Index: /issm/trunk/src/m/solutions/jpl/ControlInitialization.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/ControlInitialization.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/ControlInitialization.m (revision 2548) @@ -0,0 +1,76 @@ +function [inputs models]=ControlInitialization(models,inputs); + +%recover models +m_dh=models.dh; +m_dv=models.dv; +m_ds=models.ds; +m_dhu=models.dhu; +m_sl=models.sl; + +%recover parameters common to all solutions +verbose=m_dhu.parameters.verbose; +dim=m_dhu.parameters.dim; +ishutter=m_dhu.parameters.ishutter; +ismacayealpattyn=m_dh.parameters.ismacayealpattyn; +isstokes=m_ds.parameters.isstokes; + +%If Pattyn or MacAyeal, just return, no initialization needed +if ~isstokes + models.active=m_dh; + return +end + +%1: Get slopes once for all + +%compute slopes +displaystring(verbose,'\n%s',['computing bed slope (x and y derivatives)...']); +slopex=diagnostic_core_linear(m_sl,inputs,SlopeComputeAnalysisEnum(),BedXAnalysisEnum()); +slopey=diagnostic_core_linear(m_sl,inputs,SlopeComputeAnalysisEnum(),BedYAnalysisEnum()); +slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,m_sl.parameters,slopex,'slopex',0); +slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,m_sl.parameters,slopey,'slopey',0); + +%Add to inputs +inputs=add(inputs,'bedslopex',slopex,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes); +inputs=add(inputs,'bedslopey',slopey,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes); + +%2: update spcs using MacAyeal or Pattyn + +%horizontal velocities +displaystring(verbose,'\n%s',['computing horizontal velocities...']); +u_g=diagnostic_core_nonlinear(m_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); +displaystring(verbose,'\n%s',['extruding horizontal velocities...']); +u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,u_g,'velocity',1); + +%vertical velocities +displaystring(verbose,'\n%s',['computing vertical velocities...']); +inputs=add(inputs,'velocity',u_g_horiz,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes); +u_g_vert=diagnostic_core_linear(m_dv,inputs,DiagnosticAnalysisEnum(),VertAnalysisEnum()); + +%create 3d u_g +displaystring(verbose,'\n%s',['combining horizontal and vertical velocities...']); +u_g=zeros(m_dh.parameters.numberofnodes*3,1); +u_g(dofsetgen([1,2],3,m_dh.parameters.numberofnodes*3))=u_g_horiz; +u_g(dofsetgen([3],3,m_dh.parameters.numberofnodes*3))=u_g_vert; + +% get pressure (reconditionned) and create 4d u_g +displaystring(verbose,'\n%s',['computing pressure according to Pattyn...']); +p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,inputs); +p_g=p_g/m_ds.parameters.stokesreconditioning; +u_g_stokes=zeros(m_ds.nodesets.gsize,1); +u_g_stokes(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g; +u_g_stokes(dofsetgen([4],4,m_ds.nodesets.gsize))=p_g; +inputs=add(inputs,'velocity',u_g_stokes,'doublevec',4,m_ds.parameters.numberofnodes); + +%update BC (spcs) +displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']); +m_ds.y_g=zeros(m_ds.nodesets.gsize,1); +m_ds.y_g(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g; +[m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets); + +%Compute Stokes velocities once to have a reasonably good ug in input +displaystring(verbose,'\n%s',['computing stokes velocities and pressure ...']); +u_g=diagnostic_core_nonlinear(m_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); +inputs=add(inputs,'velocity',u_g,'doublevec',4,m_ds.parameters.numberofnodes); + +%assign output +models.active=m_ds; Index: /issm/trunk/src/m/solutions/jpl/ControlOptions.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/ControlOptions.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/ControlOptions.m (revision 2548) @@ -0,0 +1,13 @@ +function options=ControlOptions(parameters); +%CONTROLOPTIONS - setup optimization options for fmincon, using workspace parameters. +% +% Usage: +% options=ControlOptions(parameters) +% +% See also: fmincon, optimset, control + +options=optimset('fminbnd'); +options=optimset(options,'Display','iter'); +options=optimset(options,'MaxIter',parameters.maxiter(1)); +options=optimset(options,'TolX',parameters.tolx); +options=optimset(options,'FunValCheck','on'); Index: /issm/trunk/src/m/solutions/jpl/CreateFemModel.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/CreateFemModel.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/CreateFemModel.m (revision 2548) @@ -0,0 +1,37 @@ +%----------------------------------------------------------------------- +% CreateFEMModel.m: from a matlab ISSM @model class object, create a finite element +% model structure, with "c" code datasets (elements, nodes, loads, constraints, +% parameters, partitioning vectors. +% +%----------------------------------------------------------------------- + +function m=CreateFEMModel(md) + + displaystring(md.verbose,'\n reading data from model %s...',md.name); + [m.elements,m.nodes,m.constraints,m.loads,m.materials,m.parameters]=ModelProcessor(md); + + displaystring(md.verbose,'%s',' generating degrees of freedom...'); + [m.nodes,m.part,m.tpart]=Dof(m.elements,m.nodes,m.parameters); + + displaystring(md.verbose,'%s',' generating single point constraints...'); + [m.nodes,m.yg]=SpcNodes(m.nodes,m.constraints); + + displaystring(md.verbose,'%s',' generating rigid body constraints...'); + [m.Rmg,m.nodes]=MpcNodes(m.nodes,m.constraints); + + displaystring(md.verbose,'%s',' generating node sets...'); + m.nodesets=BuildNodeSets(m.nodes); + + displaystring(md.verbose,'%s',' reducing single point constraints vector...'); + [m.ys m.ys0]=Reducevectorgtos(m.yg.vector,m.nodesets); + + displaystring(md.verbose,'%s',' normalizing rigid body constraints matrix...'); + m.Gmn = NormalizeConstraints(m.Rmg,m.nodesets); + + displaystring(md.verbose,'%s',' configuring element and loads...'); + [m.elements,m.loads,m.nodes,m.parameters] = ConfigureObjects( m.elements, m.loads, m.nodes, m.materials,m.parameters); + + displaystring(md.verbose,'%s',' processing parameters...'); + m.parameters= ProcessParams(m.parameters,m.part.vector); + +end Index: /issm/trunk/src/m/solutions/jpl/SpawnCore.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/SpawnCore.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/SpawnCore.m (revision 2548) @@ -0,0 +1,63 @@ +function responses=SpawnCore(models,inputs,variables,variabledescriptors,analysis_type,sub_analysis_type,counter); +%SPAWNCORE - for Qmu analysis, using Dakota. Spawn the core solution. +% +% Usage: +% responses=SpawnCore(models,inputs,variables,variabledescriptors) +% + +%recover parameters +verbose=models.dh.parameters.verbose; +parameters=models.dh.parameters; +responsedescriptors=models.dh.parameters.responsedescriptors; +npart=models.dh.parameters.qmu_npart; +part=models.dh.parameters.qmu_part+1; %C indexing +numberofnodes=models.dh.parameters.numberofnodes; + +disp([' qmu iteration: ',num2str(counter)]); + +%first update the inputs to the models using the variables provided to us by dakota. +count=1; +while count<=numel(variables) + + descriptor=variabledescriptors{count}; + if ~qmuisdistributed(descriptor), + inputs=add(inputs,descriptor,variables(count),'double'); + + count=count+1; + else + root=qmuroot(descriptor); + param=parameters.(root); + + %next npart values in variables are partition values for this param, collect them. + partition=variables(count:count+npart-1); + + %update parameter: + param=param.*partition(part); + + %add parameter to inputs + inputs=add(inputs,root,param,'doublevec',1,numberofnodes); + + %skip next npart iterations, they all deal with the same parameter descriptor + count=count+npart; + end +end + + +%now run the core solution +if analysis_type==DiagnosticAnalysisEnum(), + + results=diagnostic_core(models,inputs); + +else + error(['SpawnCore error message: could not find core solutoin for analysis type: ' analysis_type]); +end + +%process results +processedresults=processresults(models,results); + +%now process the results to get response function values +responses=zeros(numel(responsedescriptors),1); +for i=1:numel(responsedescriptors), + descriptor=responsedescriptors{i}; + responses(i)=qmuresponse(models,results,processedresults,descriptor); +end Index: /issm/trunk/src/m/solutions/jpl/control_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/control_core.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/control_core.m (revision 2548) @@ -0,0 +1,141 @@ +function results=control_core(models,inputs) +%CONTROL_CORE - compute the core inversion +% +% Usage: +% results=control_core(models,inputs); +% + +%Preprocess models +[inputs models]=ControlInitialization(models,inputs); + +%recover active model. +model=models.active; + +%recover parameters common to all solutions +verbose=model.parameters.verbose; +dim=model.parameters.dim; +isstokes=model.parameters.isstokes; + +%recover other parameters +eps_cm=model.parameters.eps_cm; +fit=model.parameters.fit; + +%initialize control parameters +param_g=model.parameters.param_g; + +%set optimization options. +options=ControlOptions(model.parameters); + +for n=1:model.parameters.nsteps, + + %set options + options=optimset(options,'MaxFunEvals',model.parameters.maxiter(n)); + + disp(sprintf('\n%s%s%s%s\n',[' control method step ' num2str(n) '/' num2str(model.parameters.nsteps)])); + + %In case we are running a steady state control method, compute new temperature field using new parameter distribution: + if model.parameters.control_steady; + steadystate_results=steadystate_core(models,inputs); t_g=results.t_g; + inputs=add(inputs,'temperature',t_g,'doublevec',1,model.parameters.numberofnodes); + end + + %update inputs with new fit + inputs=add(inputs,'fit',model.parameters.fit(n),'double'); + inputs=add(inputs,model.parameters.control_type,param_g,'doublevec',1,model.parameters.numberofnodes); + + %Update inputs in datasets + [model.elements,model.nodes,model.loads,model.materials,model.parameters]=UpdateFromInputs(model.elements,model.nodes,model.loads,model.materials,model.parameters,inputs); + + displaystring(verbose,'\n%s',[' computing gradJ...']); + results_grad=gradjcompute_core(models,inputs); + u_g=results_grad.u_g; c(n).grad_g=results_grad.grad_g; + if dim==3, + if isstokes, + inputs=add(inputs,'velocity',u_g,'doublevec',3,model.parameters.numberofnodes); + else + inputs=add(inputs,'velocity',u_g,'doublevec',3,model.parameters.numberofnodes); + end + else + inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.numberofnodes); + end + + displaystring(verbose,'\n%s',[' normalizing directions...']); + if n>=2, + c(n).grad_g=Orth(c(n).grad_g,c(n-1).grad_g); + else + c(n).grad_g=Orth(c(n).grad_g,{}); + end + + %visualize direction. + if model.parameters.plot + plot_direction; + end + + displaystring(verbose,'\n%s',[' optimizing along gradient direction...']); + [search_scalar c(n).J]=ControlOptimization('objectivefunctionC',0,1,options,models,inputs,param_g,c(n).grad_g,n,model.parameters); + + displaystring(verbose,'\n%s',[' updating parameter using optimized search scalar...']); + param_g=param_g+search_scalar*model.parameters.optscal(n)*c(n).grad_g; + + displaystring(verbose,'\n%s',[' constraining the new distribution...']); + param_g=ControlConstrain(param_g,model.parameters); + + disp([' value of misfit J after optimization #' num2str(n) ':' num2str(c(n).J)]); + + %Has convergence been reached? + convergence=0; + if ~isnan(eps_cm), + i=n-2; + %go through the previous misfits(starting from n-2) + while (i>=1), + if (fit(i)==fit(n)), + %convergence test only if we have the same misfits + if ((c(i).J-c(n).J)/c(n).J <= eps_cm), + %convergence if convergence criteria fullfilled + convergence=1; + displaystring(verbose,'\n%s%g%s%g\n',' Convergence criterion: dJ/J = ',(c(i).J-c(n).J)/c(n).J,'<',eps_cm); + else + displaystring(verbose,'\n%s%g%s%g\n',' Convergence criterion: dJ/J = ',(c(i).J-c(n).J)/c(n).J,'>',eps_cm); + end + break; + end + i=i-1; + end + end + %stop if convergence has been reached + if (convergence), break; end + +end + +%generate output +displaystring(verbose,'\n%s',[' preparing final velocity solution...']); + +%compute final velocity from diagnostic_core (horiz+vertical) +if model.parameters.control_steady; + inputs=add(inputs,model.parameters.control_type,param_g,'doublevec',1,model.parameters.numberofnodes); + steadystate_results=steadystate_core(models,inputs); t_g=results.t_g; + u_g=steadystate_results.u_g; + t_g=steadystate_results.t_g; + m_g=steadystate_results.m_g; +else + inputs=add(inputs,model.parameters.control_type,param_g,'doublevec',1,model.parameters.numberofnodes); + results_diag=diagnostic_core(models,inputs); + u_g=results_diag.u_g; +end + +%Recover misfit at each iteration of the control method +J=zeros(length(c),1); +for i=1:length(c), + J(i)=c(i).J; +end + +%build results +results.time=0; +results.step=1; +results.J=J; +results.param_g=param_g; +results.u_g=u_g; +if model.parameters.control_steady, + results.t_g=t_g; + results.m_g=m_g; +end Index: /issm/trunk/src/m/solutions/jpl/convergence.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/convergence.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/convergence.m (revision 2548) @@ -0,0 +1,81 @@ +function converged=convergence(K_ff,p_f,u_f,u_f_old,params) + +%Get convergence options +verbose=params.verbose; +yts=params.yts; +eps_res=params.eps_res; +eps_rel=params.eps_rel; +eps_abs=params.eps_abs; + +%initialization +converged=0; +displaystring(verbose,' '); + +%Display solver caracteristics +if (verbose>1), + solver_res=norm(K_ff*u_f-p_f,2)/norm(p_f,2); + displaystring(verbose>1,'%s%g',' condition number of stiffness matrix: ',condest(K_ff)); + displaystring(verbose>1,'%s%g',' solver residue: norm(KU-F)/norm(F)=',solver_res); +end + +%Force equilibrium (Mandatory) +res=norm(K_ff*u_f_old-p_f,2)/norm(p_f,2); +if isnan(res), + error('convergence error message: mechanical equilibrium convergence criterion is NaN!'); +end +if (res<=eps_res), + displaystring(verbose,'%-60s%g%s%g%s',' mechanical equilibrium convergence criterion',res*100,' < ',eps_res*100,' %'); + converged=1; +else + displaystring(verbose,'%-60s%g%s%g%s',' mechanical equilibrium convergence criterion',res*100,' > ',eps_res*100,' %'); + converged=0; +end + +%Relative criterion (optional) +if ((~isnan(eps_rel)) | (verbose>1)), + + %compute ndu/nu + duf=u_f-u_f_old; + ndu=norm(duf,2); + nu=norm(u_f_old,2); + if isnan(ndu), + error('convergence error message: relative convergence criterion is NaN!'); + end + + %print criterion + if ~isnan(eps_rel), + if (ndu/nu<=eps_rel), + displaystring(verbose,'%-60s%g%s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' < ',eps_rel*100,' %'); + else + displaystring(verbose,'%-60s%g%s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' > ',eps_rel*100,' %'); + converged=0; + end + else + displaystring(verbose,'%-60s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' %'); + end + +end + +%Absolute criterion (optional) +if ((~isnan(eps_abs)) | (verbose>1)), + + %compute max(du) + duf=u_f-u_f_old; + nduinf=norm(duf,inf)*yts; + if isnan(nduinf), + error('convergence error message: absolute convergence criterion is NaN!'); + end + + %print criterion + if ~isnan(eps_abs), + if (nduinf<=eps_abs), + displaystring(verbose,'%-60s%g%s%g%s',' absolute convergence criterion: max(du)',nduinf,' < ',eps_abs,' m/yr'); + else + displaystring(verbose,'%-60s%g%s%g%s',' absolute convergence criterion: max(du)',nduinf,' > ',eps_abs,' m/yr'); + converged=0; + end + else + displaystring(verbose,'%-60s%g%s',' absolute convergence criterion: max(du)',nduinf,' m/yr'); + end + +end Index: /issm/trunk/src/m/solutions/jpl/diagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/diagnostic.m (revision 2548) @@ -0,0 +1,62 @@ +function md=diagnostic(md); +%DIAGNOSTIC - compute the velocity field of a model +% +% Usage: +% md=diagnostic(md) +% + %timing + t1=clock; + + %Build all models requested for diagnostic simulation + models.analysis_type=DiagnosticAnalysisEnum; %needed for processresults + + displaystring(md.verbose,'%s',['reading diagnostic horiz model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HorizAnalysisEnum; models.dh=CreateFemModel(md); + + displaystring(md.verbose,'\n%s',['reading diagnostic vert model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=VertAnalysisEnum; models.dv=CreateFemModel(md); + + displaystring(md.verbose,'\n%s',['reading diagnostic stokes model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=StokesAnalysisEnum; models.ds=CreateFemModel(md); + + displaystring(md.verbose,'\n%s',['reading diagnostic hutter model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HutterAnalysisEnum; models.dhu=CreateFemModel(md); + + displaystring(md.verbose,'\n%s',['reading surface and bed slope computation model data']); + md.analysis_type=SlopeComputeAnalysisEnum; md.sub_analysis_type=NoneAnalysisEnum; models.sl=CreateFemModel(md); + + % figure out number of dof: just for information purposes. + md.dof=modelsize(models); + + %initialize inputs + inputs=inputlist; + inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes); + if md.control_analysis, + inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes); + end + + %compute solution + if ~models.dh.parameters.qmu_analysis, + if md.control_analysis, + %launch core of control solution. + results=control_core(models,inputs); + + %process results + if ~isstruct(md.results), md.results=struct(); end + md.results.diagnostic=processresults(models,results); + else, + %launch core of diagnostic solution. + results=diagnostic_core(models,inputs); + + %process results + if ~isstruct(md.results), md.results=struct(); end + md.results.diagnostic=processresults(models,results); + end + else + %launch dakota driver for diagnostic core solution + Qmu(models,inputs,models.dh.parameters); + end + + %stop timing + t2=clock; + displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic_core.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/diagnostic_core.m (revision 2548) @@ -0,0 +1,123 @@ +function results=diagnostic_core(models,inputs); +%DIAGNOSTIC_CORE - compute the core velocity field +% +% Usage: +% results=diagnostic_core(models,inputs); +% + +%recover models +m_dh=models.dh; +m_dv=models.dv; +m_ds=models.ds; +m_dhu=models.dhu; +m_sl=models.sl; + +%recover parameters common to all solutions +verbose=m_dhu.parameters.verbose; +dim=m_dh.parameters.dim; +ishutter=m_dhu.parameters.ishutter; +ismacayealpattyn=m_dh.parameters.ismacayealpattyn; +isstokes=m_ds.parameters.isstokes; +numrifts=m_dhu.parameters.numrifts; + +debug=1; + +if ishutter, + + displaystring(verbose,'\n%s',['computing surface slope (x and y derivatives)...']); + slopex=diagnostic_core_linear(m_sl,inputs,SlopeComputeAnalysisEnum(),SurfaceXAnalysisEnum()); + slopey=diagnostic_core_linear(m_sl,inputs,SlopeComputeAnalysisEnum(),SurfaceYAnalysisEnum()); + + if dim==3, + displaystring(verbose,'\n%s',['extruding slopes in 3d...']); + slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,m_sl.parameters,slopex,'slopex',0); + slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,m_sl.parameters,slopey,'slopey',0); + end + + displaystring(verbose,'\n%s',['computing slopes...']); + inputs=add(inputs,'surfaceslopex',slopex,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes); + inputs=add(inputs,'surfaceslopey',slopey,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes); + + displaystring(verbose,'\n%s',['computing hutter velocities...']); + u_g=diagnostic_core_linear(m_dhu,inputs,DiagnosticAnalysisEnum(),HutterAnalysisEnum()); + + displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']); + p_g=ComputePressure(m_dhu.elements,m_dhu.nodes,m_dhu.loads,m_dhu.materials,m_dhu.parameters,inputs); + + displaystring(verbose,'\n%s',['update boundary conditions for macyeal pattyn using hutter results...']); + if ismacayealpattyn, + m_dh.y_g=u_g; + [m_dh.ys m_dh.ys0]=Reducevectorgtos(m_dh.y_g,m_dh.nodesets); + end + +end + +if ismacayealpattyn, + + displaystring(verbose,'\n%s',['computing horizontal velocities...']); + [u_g m_dh.loads]=diagnostic_core_nonlinear(m_dh,inputs,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); + + displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']); + p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,inputs); +end + +if dim==3, + + displaystring(verbose,'\n%s',['extruding horizontal velocities...']); + u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,u_g,'velocity',1); + + displaystring(verbose,'\n%s',['computing vertical velocities...']); + inputs=add(inputs,'velocity',u_g_horiz,'doublevec',m_dh.parameters.numberofdofspernode,m_dh.parameters.numberofnodes); + u_g_vert=diagnostic_core_linear(m_dv,inputs,DiagnosticAnalysisEnum(),VertAnalysisEnum()); + + displaystring(verbose,'\n%s',['combining horizontal and vertical velocities...']); + u_g=zeros(m_dh.parameters.numberofnodes*3,1); + u_g(dofsetgen([1,2],3,m_dh.parameters.numberofnodes*3))=u_g_horiz; + u_g(dofsetgen([3],3,m_dh.parameters.numberofnodes*3))=u_g_vert; + + displaystring(verbose,'\n%s',['computing pressure according to Pattyn...']); + p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.loads,m_dh.materials,m_dh.parameters,inputs); + + if isstokes, + + %"recondition" pressure + p_g=p_g/m_ds.parameters.stokesreconditioning; + + displaystring(verbose,'\n%s',['computing bed slope (x and y derivatives)...']); + slopex=diagnostic_core_linear(m_sl,inputs,SlopeComputeAnalysisEnum(),BedXAnalysisEnum()); + slopey=diagnostic_core_linear(m_sl,inputs,SlopeComputeAnalysisEnum(),BedYAnalysisEnum()); + slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,m_sl.parameters,slopex,'slopex',0); + slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.loads,m_sl.materials,m_sl.parameters,slopey,'slopey',0); + + inputs=add(inputs,'bedslopex',slopex,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes); + inputs=add(inputs,'bedslopey',slopey,'doublevec',m_sl.parameters.numberofdofspernode,m_sl.parameters.numberofnodes); + + %recombine u_g and p_g: + u_g_stokes=zeros(m_ds.nodesets.gsize,1); + u_g_stokes(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g; + u_g_stokes(dofsetgen([4],4,m_ds.nodesets.gsize))=p_g; + + inputs=add(inputs,'velocity',u_g_stokes,'doublevec',4,m_ds.parameters.numberofnodes); + + displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']); + m_ds.y_g=zeros(m_ds.nodesets.gsize,1); + m_ds.y_g(dofsetgen([1,2,3],4,m_ds.nodesets.gsize))=u_g; + [m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets); + + displaystring(verbose,'\n%s',['computing stokes velocities and pressure ...']); + u_g=diagnostic_core_nonlinear(m_ds,inputs,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); + + %"decondition" pressure + p_g=u_g(4:4:end)*m_dh.parameters.stokesreconditioning; + end +end + +%load onto results +results.step=1; +results.time=0; +results.u_g=u_g; +results.p_g=p_g; + +if numrifts, + results.riftproperties=OutputRifts(m_dh.loads,m_dh.parameters); +end Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m (revision 2548) @@ -0,0 +1,30 @@ +function u_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type) +%DIAGNOSTIC_CORE_LINEAR - linear solution sequence +% +% Usage: +% u_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type) + + %stiffness and load generation only: + m.parameters.kflag=1; m.parameters.pflag=1; + + %Update inputs in datasets + [m.elements,m.nodes, m.loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes, m.loads,m.materials,m.parameters,inputs); + + %system matrices + [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + [K_gg, p_g,kmax]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + + %Reduce tangent matrix from g size to f size + [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); + displaystring(m.parameters.verbose>1,'%s%g',' condition number of stiffness matrix: ',condest(K_ff)); + + %Reduce load from g size to f size + [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets); + + %Solve + [u_f]=Solver(K_ff,p_f,[],m.parameters); + + %Merge back to g set + [u_g]= Mergesolutionfromftog( u_f, m.Gmn, m.ys, m.nodesets ); + +end %end function Index: /issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/diagnostic_core_nonlinear.m (revision 2548) @@ -0,0 +1,94 @@ +function [u_g varargout]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type) +%DIAGNOSTIC_CORE_NONLINEAR - core solution of diagnostic non-linear +% +% Usage: +% [u_g varargout]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type) + +% first off! We are going to modify the loads dataset. We need to shield the loads from those changes once we return; + loads=m.loads; + +% stiffness and load generation only: + m.parameters.kflag=1; m.parameters.pflag=1; + +% initialize solution vector + converged=0; count=1; + soln(count).u_g=get(inputs,'velocity',[1 1 (m.parameters.numberofdofspernode>=3) (m.parameters.numberofdofspernode>=3)]); %only keep vz if running with more than 3 dofs per node + soln(count).u_f= Reducevectorgtof(soln(count).u_g,m.nodesets); + + %add converged=0 flag first in inputs. + inputs=add(inputs,'converged',converged,'double'); + + displaystring(m.parameters.verbose,'\n%s',[' starting direct shooting method']); + while(~converged), + + %add velocity to inputs + inputs=add(inputs,'velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); + if count>1, + inputs=add(inputs,'old_velocity',soln(count-1).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); + else + inputs=add(inputs,'old_velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); + end + + %Update inputs in datasets + [m.elements,m.nodes, loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,m.parameters,inputs); + + %system matrices + [K_gg_nopenalty , p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + + %penalties + [K_gg , p_g, kmax]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + + %Reduce tangent matrix from g size to f size + [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); + + %Reduce load from g size to f size + [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets); + + %Increment index + count=count+1; + + %Solve + [soln(count).u_f]=Solver(K_ff,p_f,[],m.parameters); + + %Merge back to g set + [soln(count).u_g]= Mergesolutionfromftog( soln(count).u_f, m.Gmn, m.ys, m.nodesets ); + + %Deal with penalty loads + inputs=add(inputs,'velocity',soln(count).u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); + + %penalty constraints + [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + + displaystring(m.parameters.verbose,'%s%i',' number of unstable constraints: ',num_unstable_constraints); + + %Figure out if convergence have been reached + converged=convergence(K_ff,p_f,soln(count).u_f,soln(count-1).u_f,m.parameters); + + %add convergence status into inputs + inputs=add(inputs,'converged',converged,'double'); + + %rift convergence criterion + if ~constraints_converged, + converged=0; + end + + end + + %some cleanup + u_g=soln(count).u_g; + clear soln + + %more output might be needed, when running in control_core.m + nout=max(nargout,1)-1; + if nout==2, + inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); + m.parameters.kflag=1; m.parameters.pflag=0; + [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); + varargout(1)={K_ff}; + varargout(2)={K_fs}; + end + if nout==1, + varargout(1)={loads}; + end +end Index: /issm/trunk/src/m/solutions/jpl/dofsetgen.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/dofsetgen.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/dofsetgen.m (revision 2548) @@ -0,0 +1,19 @@ +function dofset=dofsetgen(doflist,numdofs,dofsize) +%DOFSETGEN - generate a dof set to index a solution vector +% +% Usage: +% dofset=dofsetgen(doflist,numdofs,dofsize) +% +% ex: +% dof1set=dofsetgen([1],3,10); will generate 1,4,7,10 +% dof2set=dofsetgen([2],4,10); will generate 2,6,10 +% dof12set=dofsetgen([1,2],5,15); will generate 1,2,6,7,11,12 + + +flags=zeros(dofsize,1); + +for i=1:length(doflist), + flags(doflist(i):numdofs:end)=1; +end + +dofset=find(flags); Index: /issm/trunk/src/m/solutions/jpl/gradjcompute_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/gradjcompute_core.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/gradjcompute_core.m (revision 2548) @@ -0,0 +1,53 @@ +function results=gradjcompute_core(models,inputs); + +%recover active models +m=models.active; + +%recover parameters +verbose=m.parameters.verbose; +dim=m.parameters.dim; +extrude_param=m.parameters.extrude_param; +ishutter=m.parameters.ishutter; +ismacayealpattyn=m.parameters.ismacayealpattyn; +isstokes=m.parameters.isstokes; +analysis_type=m.parameters.analysis_type; +sub_analysis_type=m.parameters.sub_analysis_type; + +%Recover solution for this stiffness and right hand side: +displaystring(verbose,'%s',' computing velocities...'); +[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type); +inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); + + +%Buid Du, difference between observed velocity and model velocity. +displaystring(verbose,'%s',' computing Du...'); +[Du_g]=Du(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + +%Reduce adjoint load from g-set to f-set +[Du_f] = Reduceloadfromgtof( Du_g, m.Gmn, K_fs0, m.ys0, m.nodesets); + +%Solve for adjoint vector: +displaystring(verbose,'%s',' computing adjoint state...'); +lambda_f=Solver(K_ff0,Du_f,[],m.parameters); + +%Merge back to g set +lambda_g= Mergesolutionfromftog( lambda_f, m.Gmn, m.ys0, m.nodesets ); +inputs=add(inputs,'adjoint',lambda_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); + +%Compute gradJ +grad_g=Gradj(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); +if (dim==3 & extrude_param), + displaystring(verbose,'%s',' extruding gradient...'); + grad_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,m.parameters,grad_g,'gradj',0); +end + +if m.parameters.control_steady, + %compute initial velocity from diagnostic_core (horiz+vertical) + displaystring(verbose,'%s',' compute 3d initial velocity...'); + results_diag=diagnostic_core(models,inputs); + u_g=results_diag.u_g; +end + +%Assign output +results.u_g=u_g; +results.grad_g=grad_g; Index: /issm/trunk/src/m/solutions/jpl/modelsize.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/modelsize.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/modelsize.m (revision 2548) @@ -0,0 +1,20 @@ +function dof=modelsize(models) +%DOF - return the maximum number of degrees of freedom of the model +% +% Usage: +% dof=modelsize(models) + +%initialize dof +dof=0; + +%get all models +modelsname=fieldnames(models); + +%get max dof +for i=1:length(modelsname); + if ~strcmpi(modelsname,'analysis_type'), + if ~isempty(models.(modelsname{i}).nodesets), + dof=max(dof,models.(modelsname{i}).nodesets.fsize); + end + end +end Index: /issm/trunk/src/m/solutions/jpl/objectivefunctionC.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/objectivefunctionC.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/objectivefunctionC.m (revision 2548) @@ -0,0 +1,25 @@ +function J =objectivefunctionC(search_scalar,models,inputs,p_g,grad_g,n,analysis_type,sub_analysis_type); + +%recover active model. +m=models.active; + +%recover some parameters +optscal=m.parameters.optscal(n); +fit=m.parameters.fit(n); +control_type=m.parameters.control_type; +analysis_type=m.parameters.analysis_type; + +%Update along gradient using scalar supplied by fmincon optimization routine +parameter=p_g+search_scalar*optscal*grad_g; + +%Plug parameter into inputs +inputs=add(inputs,m.parameters.control_type,parameter,'doublevec',1,m.parameters.numberofnodes); + +%Run diagnostic with updated parameters. +u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type); + +%add velocity to inputs. +inputs=add(inputs,'velocity',u_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); + +%Compute misfit for this velocity field. +J=Misfit(m.elements,m.nodes,m.loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); Index: /issm/trunk/src/m/solutions/jpl/plot_direction.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/plot_direction.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/plot_direction.m (revision 2548) @@ -0,0 +1,1 @@ +plotmodel(md,'data',c(n).grad_g(1:2:end),'title',['Normalized Direction for ' m.parameters.control_type],'figure',1,'colorbar#all','on'); pause(1); Index: /issm/trunk/src/m/solutions/jpl/plot_newdistribution.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/plot_newdistribution.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/plot_newdistribution.m (revision 2548) @@ -0,0 +1,1 @@ +plotmodel(md,'data',p_g(1:2:end),'title',['Distribution of ' m.parameters.control_type ' at iteration ' num2str(n)],'figure',2,'colorbar#all','on'); pause(1); Index: /issm/trunk/src/m/solutions/jpl/processresults.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/processresults.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/processresults.m (revision 2548) @@ -0,0 +1,136 @@ +function newresults=processresults(models,results) +%PROCESSRESULTS - process results from core solution +% +% The solution (vel,pressure,...) are on the g-set, +% use index to get the value on each node +% +% Usage: +% results=processresults(results,models,analysis_type) + +%initialize output +newresults=struct(); + +%get analysis_type: +analysis_type=models.analysis_type; + +%recover models first +if (analysis_type==DiagnosticAnalysisEnum | analysis_type==TransientAnalysisEnum | analysis_type==SteadystateAnalysisEnum()), + m_dh=models.dh; + m_ds=models.ds; + m_dhu=models.dhu; + + %some flags needed + dim=m_dhu.parameters.dim; + ishutter=m_dhu.parameters.ishutter; + ismacayealpattyn=m_dh.parameters.ismacayealpattyn; + isstokes=m_ds.parameters.isstokes; +end +if (analysis_type==ControlAnalysisEnum()), + m_dh=models.dh; + + %some flags needed + dim=m_dh.parameters.dim; + ishutter=m_dh.parameters.ishutter; + ismacayealpattyn=m_dh.parameters.ismacayealpattyn; + isstokes=m_dh.parameters.isstokes; +end +if (analysis_type==ThermalAnalysisEnum() | analysis_type==TransientAnalysisEnum() | analysis_type==SteadystateAnalysisEnum()), + m_m=models.m; +end + + +%go through results +for i=1:length(results), + + %get field names of results + resultsname=fieldnames(results(i)); + + %go through all field + for j=1:length(resultsname), + + if strcmpi(resultsname{j},'u_g'), + + u_g=results(i).u_g; + yts=m_dh.parameters.yts; + + if dim==2, + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).vx=u_g(1:2:end)*yts; + newresults(i).vy=u_g(2:2:end)*yts; + newresults(i).vel=sqrt(newresults(i).vx.^2+newresults(i).vy.^2); + + else + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + if isstokes, + newresults(i).vx=u_g(1:4:end)*yts; + newresults(i).vy=u_g(2:4:end)*yts; + newresults(i).vz=u_g(3:4:end)*yts; + else + newresults(i).vx=u_g(1:3:end)*yts; + newresults(i).vy=u_g(2:3:end)*yts; + newresults(i).vz=u_g(3:3:end)*yts; + end + newresults(i).vel=sqrt(newresults(i).vx.^2+newresults(i).vy.^2+newresults(i).vz.^2); + end + + elseif strcmpi(resultsname{j},'p_g'), + + p_g=results(i).p_g; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).pressure=p_g; + + elseif strcmpi(resultsname{j},'t_g'), + + t_g=results(i).t_g; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).temperature=t_g; + + elseif strcmpi(resultsname{j},'m_g'), + + m_g=results(i).m_g; + yts=m_m.parameters.yts; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).melting=m_g*yts; + + elseif strcmpi(resultsname{j},'h_g'), + + h_g=results(i).h_g; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).thickness=h_g; + + elseif strcmpi(resultsname{j},'s_g'), + + s_g=results(i).s_g; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).surface=s_g; + + elseif strcmpi(resultsname{j},'b_g'), + + b_g=results(i).b_g; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).bed=b_g; + + elseif strcmpi(resultsname{j},'param_g'), + + param_g=results(i).param_g; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).parameter=param_g; + + else + + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).(resultsname{j})=results(i).(resultsname{j}); + + end + end +end Index: /issm/trunk/src/m/solutions/jpl/prognostic.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/prognostic.m (revision 2548) @@ -0,0 +1,37 @@ +function md=prognostic(md) +%PROGNOSITC - prognostic solution sequence. +% +% Usage: +% md=prognostic(md) + + %timing + t1=clock; + + %Build all models requested for diagnostic simulation + models.analysis_type=PrognosticAnalysisEnum; %needed for processresults + + displaystring(md.verbose,'%s',['reading prognostic model data']); + md.analysis_type=PrognosticAnalysisEnum; models.p=CreateFemModel(md); + + % figure out number of dof: just for information purposes. + md.dof=modelsize(models); + + %initialize inputs + displaystring(md.verbose,'\n%s',['setup inputs...']); + inputs=inputlist; + inputs=add(inputs,'velocity',models.p.parameters.u_g,'doublevec',3,models.p.parameters.numberofnodes); + inputs=add(inputs,'thickness',models.p.parameters.h_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.numberofnodes); + inputs=add(inputs,'dt',models.p.parameters.dt*models.p.parameters.yts,'double'); + + displaystring(md.verbose,'\n%s',['call computational core:']); + results=prognostic_core(models,inputs,PrognosticAnalysisEnum(),NoneAnalysisEnum()); + + displaystring(md.verbose,'\n%s',['load results...']); + if ~isstruct(md.results), md.results=struct(); end + md.results.prognostic=processresults(models, results); + + %stop timing + t2=clock; + displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: /issm/trunk/src/m/solutions/jpl/prognostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic_core.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/prognostic_core.m (revision 2548) @@ -0,0 +1,24 @@ +function results=prognostic_core(models,inputs,analysis_type,sub_analysis_type) +%PROGNOSTIC_CORE - linear solution sequence +% +% Usage: +% h_g=prognostic_core(m,inputs,analysis_type,sub_analysis_type) + + %get FE model + m=models.p; + results.time=0; + results.step=1; + + displaystring(m.parameters.verbose,'\n%s',['depth averaging velocity...']); + %Take only the first two dofs of m.parameters.u_g + u_g=get(inputs,'velocity',[1 1 0 0]); + u_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,u_g,'velocity'); + inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.numberofnodes); + + displaystring(m.parameters.verbose,'\n%s',['call computational core:']); + results.h_g=diagnostic_core_linear(m,inputs,analysis_type,sub_analysis_type); + + displaystring(m.parameters.verbose,'\n%s',['extrude computed thickness on all layers:']); + results.h_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness',0); + +end %end function Index: /issm/trunk/src/m/solutions/jpl/steadystate.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/steadystate.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/steadystate.m (revision 2548) @@ -0,0 +1,76 @@ +function md=steadystate(md); +%STEADYSTATE - compute the velocity and temperature field of a model in steady state. +% +% Usage: +% md=steadystate(md) +% + + %timing + t1=clock; + + models.analysis_type=SteadystateAnalysisEnum; %needed for processresults + + %Build all models requested for diagnostic simulation + displaystring(md.verbose,'%s',['reading diagnostic horiz model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HorizAnalysisEnum; models.dh=CreateFemModel(md); + + displaystring(md.verbose,'\n%s',['reading diagnostic vert model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=VertAnalysisEnum; models.dv=CreateFemModel(md); + + displaystring(md.verbose,'\n%s',['reading diagnostic stokes model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=StokesAnalysisEnum; models.ds=CreateFemModel(md); + + displaystring(md.verbose,'\n%s',['reading diagnostic hutter model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HutterAnalysisEnum; models.dhu=CreateFemModel(md); + + displaystring(md.verbose,'\n%s',['reading surface and bed slope computation model data']); + md.analysis_type=SlopeComputeAnalysisEnum; md.sub_analysis_type=NoneAnalysisEnum; models.sl=CreateFemModel(md); + + %Build all models requested for thermal simulation + displaystring(md.verbose,'%s',['reading thermal model data']); + md.analysis_type=ThermalAnalysisEnum(); md.sub_analysis_type=NoneAnalysisEnum(); models.t=CreateFemModel(md); + + displaystring(md.verbose,'%s',['reading melting model data']); + md.analysis_type=MeltingAnalysisEnum(); md.sub_analysis_type=NoneAnalysisEnum(); models.m=CreateFemModel(md); + + % figure out number of dof: just for information purposes. + md.dof=modelsize(models); + + %initialize inputs + displaystring(md.verbose,'\n%s',['setup inputs...']); + inputs=inputlist; + inputs=add(inputs,'velocity',models.dh.parameters.u_g,'doublevec',3,models.dh.parameters.numberofnodes); + inputs=add(inputs,'pressure',models.t.parameters.p_g,'doublevec',1,models.t.parameters.numberofnodes); + inputs=add(inputs,'dt',models.t.parameters.dt*models.t.parameters.yts,'double'); + if md.control_analysis, + inputs=add(inputs,'velocity_obs',models.dh.parameters.u_g_obs,'doublevec',2,models.dh.parameters.numberofnodes); + end + + %compute solution + if ~models.dh.parameters.qmu_analysis, + if md.control_analysis, + %change control_steady to 1 + models.dh.parameters.control_steady=1; + models.ds.parameters.control_steady=1; + %launch core of control solution. + results=control_core(models,inputs); + + %process results + if ~isstruct(md.results), md.results=struct(); end + md.results.steadystate=processresults(models,results); + else, + %launch core of steadystate solution. + results=steadystate_core(models,inputs); + + %process results + if ~isstruct(md.results), md.results=struct(); end + md.results.steadystate=processresults(models,results); + end + else + %launch dakota driver for steadystate core solution + Qmu(models,inputs,models.dh.parameters); + end + + %stop timing + t2=clock; + displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: /issm/trunk/src/m/solutions/jpl/steadystate_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/steadystate_core.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/steadystate_core.m (revision 2548) @@ -0,0 +1,87 @@ +function results=steadystate_core(models,inputs); +%STEADYSTATE_CORE - compute the core temperature and velocity field at thermal steady state. +% +% Usage: +% results=steadystate_core(models,inputs); +% + +%recover models +m_dh=models.dh; +m_dv=models.dv; +m_ds=models.ds; +m_dhu=models.dhu; +m_sl=models.sl; +m_t=models.t; +m_m=models.m; + +%recover parameters common to all solutions +verbose=m_dhu.parameters.verbose; +dim=m_dhu.parameters.dim; +eps_rel=m_dhu.parameters.eps_rel; +ishutter=m_dhu.parameters.ishutter; +ismacayealpattyn=m_dh.parameters.ismacayealpattyn; +isstokes=m_ds.parameters.isstokes; + +%convergence +converged=0; + +step=1; + +%initialize: +t_g_old=0; +u_g_old=0; + +if isstokes, ndof=4; else ndof=3; end + +while(~converged), + + displaystring(verbose,'%s%i','computing temperature and velocity for step ',step); + + %first compute temperature at steady state. + if (step>1), + inputs=add(inputs,'velocity',results_diagnostic.u_g,'doublevec',ndof,m_t.parameters.numberofnodes); + end + results_thermal=thermal_core(models,inputs); + + %add temperature to inputs. + %Compute depth averaged temperature + temperature_average=FieldDepthAverage(m_t.elements,m_t.nodes,m_t.loads,m_t.materials,m_t.parameters,results_thermal.t_g,'temperature'); + inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,m_t.parameters.numberofnodes); + inputs=add(inputs,'temperature',results_thermal.t_g,'doublevec',1,m_t.parameters.numberofnodes); + + %now compute diagnostic velocity using the steady state temperature. + results_diagnostic=diagnostic_core(models,inputs); + + %convergence? + du_g=results_diagnostic.u_g-u_g_old; + ndu=norm(du_g,2); + nu=norm(u_g_old,2); + + dt_g=results_thermal.t_g-t_g_old; + ndt=norm(dt_g,2); + nt=norm(t_g_old,2); + + u_g_old=results_diagnostic.u_g; + t_g_old=results_thermal.t_g; + + displaystring(verbose,'%-60s%g\n %s%g\n %s%g%s',... + ' relative convergence criterion: velocity -> norm(du)/norm(u)= ',ndu/nu*100,' temperature -> norm(dt)/norm(t)=',ndt/nt*100,' eps_rel: ',eps_rel*100,' %'); + if (ndu/nu<=eps_rel) & (ndt/nt<=eps_rel), + converged=1; + else + converged=0; + end + + step=step+1; + if converged, + break; + end +end + +%save results from thermal and diagnostic +results.u_g=results_diagnostic.u_g; +results.p_g=results_diagnostic.p_g; +results.t_g=results_thermal.t_g; +results.m_g=results_thermal.m_g; +results.step=step; +results.time=0; Index: /issm/trunk/src/m/solutions/jpl/thermal.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/thermal.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/thermal.m (revision 2548) @@ -0,0 +1,44 @@ +function md=thermal(md) +%THERMAL - thermal solution sequence: steady state and transient +% +% Usage: +% md=thermal(md) + + %timing + t1=clock; + + models.analysis_type=ThermalAnalysisEnum(); %needed for processresults + + %Build all models requested for diagnostic simulation + displaystring(md.verbose,'%s',['reading thermal model data']); + md.analysis_type=ThermalAnalysisEnum(); models.t=CreateFemModel(md); + + displaystring(md.verbose,'%s',['reading melting model data']); + md.analysis_type=MeltingAnalysisEnum(); models.m=CreateFemModel(md); + + % figure out number of dof: just for information purposes. + md.dof=modelsize(models); + + %initialize inputs + displaystring(md.verbose,'\n%s',['setup inputs...']); + inputs=inputlist; + inputs=add(inputs,'velocity',models.t.parameters.u_g,'doublevec',3,models.t.parameters.numberofnodes); + inputs=add(inputs,'pressure',models.t.parameters.p_g,'doublevec',1,models.t.parameters.numberofnodes); + inputs=add(inputs,'dt',models.t.parameters.dt*models.t.parameters.yts,'double'); + + %compute solution + if ~models.t.parameters.qmu_analysis, + %launch core of diagnostic solution. + results=thermal_core(models,inputs); + + %process results + if ~isstruct(md.results), md.results=struct(); end + md.results.thermal=processresults(models,results); + else + %launch dakota driver for diagnostic core solution + Qmu(models,inputs,models.t.parameters); + end + + %stop timing + t2=clock; + displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: /issm/trunk/src/m/solutions/jpl/thermal_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/thermal_core.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/thermal_core.m (revision 2548) @@ -0,0 +1,54 @@ +function results=thermal_core(models,inputs) +%THERMAL_CORE - core of thermal solution +% +% Usage: +% solution=thermal_core(models,inputs) + +%get FE model +m_t=models.t; +m_m=models.m; + +if m_t.parameters.dt==0, + + results.time=0; + results.step=1; + + displaystring(m_t.parameters.verbose,'\n%s',['computing temperatures...']); + [results.t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,ThermalAnalysisEnum(),NoneAnalysisEnum()); + + displaystring(m_t.parameters.verbose,'\n%s',['computing melting...']); + inputs=add(inputs,'melting_offset',melting_offset,'double'); + inputs=add(inputs,'temperature',results.t_g,'doublevec',1,m_t.parameters.numberofnodes); + results.m_g=diagnostic_core_linear(m_m,inputs,MeltingAnalysisEnum(),NoneAnalysisEnum()); + +else + + %initialize temperature and melting + t_g=m_t.parameters.t_g; + m_g=m_m.parameters.m_g; + nsteps=m_t.parameters.ndt/m_t.parameters.dt; + + %initialize temperature and melting + results.step=1; + results.time=0; + results.t_g=t_g; + results.m_g=m_g; + + for n=1:nsteps, + + displaystring(m_t.parameters.verbose,'\n%s%i/%i\n','time step: ',n,nsteps); + results(n+1).step=n+1; + results(n+1).time=n*m_t.parameters.dt; + + displaystring(m_t.parameters.verbose,'\n%s',[' computing temperatures...']); + inputs=add(inputs,'temperature',results(n).t_g,'doublevec',1,m_t.parameters.numberofnodes); + [results(n+1).t_g m_t.loads melting_offset]=thermal_core_nonlinear(m_t,inputs,ThermalAnalysisEnum(),NoneAnalysisEnum()); + + displaystring(m_t.parameters.verbose,'\n%s',[' computing melting...']); + inputs=add(inputs,'temperature',results(n+1).t_g,'doublevec',1,m_t.parameters.numberofnodes); + inputs=add(inputs,'melting_offset',melting_offset,'double'); + results(n+1).m_g=diagnostic_core_linear(m_m,inputs,MeltingAnalysisEnum(),NoneAnalysisEnum()); + + end + +end Index: /issm/trunk/src/m/solutions/jpl/thermal_core_nonlinear.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/thermal_core_nonlinear.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/thermal_core_nonlinear.m (revision 2548) @@ -0,0 +1,76 @@ +function [t_g ,loads, melting_offset]=thermal_core_nonlinear(m,inputs,analysis_type,sub_analysis_type) +%THERMAL_CORE_NONLINEAR - core of thermal solution sequence. +% model is return together with temperature +% +% +% Usage: +% [t_g ,loads, melting_offset]=thermal_core_nonlinear(m,inputs,analysis_type,sub_analysis_type); +% +% + + count=1; + converged=0; + +% we are going to return the loads, make them a variable of this routine + loads=m.loads; + + %stiffness and load generation only: + m.parameters.kflag=1; m.parameters.pflag=1; + + displaystring(m.parameters.verbose,'\n%s',[' starting direct shooting method']); + + while(~converged), + + %Update inputs in datasets + [m.elements,m.nodes, loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,m.parameters,inputs); + + %system matrices + if ~m.parameters.lowmem + if count==1 + displaystring(m.parameters.verbose,'%s',[' system matrices']); + [K_gg_nopenalty, p_g_nopenalty]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + end + displaystring(m.parameters.verbose,'%s',[' penalty system matrices']); + [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg_nopenalty,p_g_nopenalty,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + else + displaystring(m.parameters.verbose,'%s',[' system matrices']); + [K_gg , p_g]=SystemMatrices(m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + displaystring(m.parameters.verbose,'%s',[' penalty system matrices']); + [K_gg , p_g, melting_offset]=PenaltySystemMatrices(K_gg,p_g,m.elements,m.nodes,loads,m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + end + + %Reduce tangent matrix from g size to f size + [K_ff, K_fs] = Reducematrixfromgtof( K_gg, m.Gmn, m.nodesets); + + %Reduce load from g size to f size + [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets); + + %Solve + displaystring(m.parameters.verbose,'%s%g',' condition number of stiffness matrix: ',condest(K_ff)); + [t_f]=Solver(K_ff,p_f,[],m.parameters); + + %Merge back to g set + displaystring(m.parameters.verbose,'%s',[' merging solution back to g set']); + [t_g]= Mergesolutionfromftog( t_f, m.Gmn, m.ys, m.nodesets ); + + %Update inputs in datasets + inputs=add(inputs,'temperature',t_g,'doublevec',m.parameters.numberofdofspernode,m.parameters.numberofnodes); + displaystring(m.parameters.verbose,'%s',[' update inputs']); + [m.elements,m.nodes, loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes, loads,m.materials,m.parameters,inputs); + + %penalty constraints + displaystring(m.parameters.verbose,'%s',[' penalty constraints']); + [loads,constraints_converged,num_unstable_constraints] =PenaltyConstraints( m.elements,m.nodes, loads, m.materials,m.parameters,inputs,analysis_type,sub_analysis_type); + + if ~converged, + displaystring(m.parameters.verbose,'%s%i',' #unstable constraints ',num_unstable_constraints); + + if num_unstable_constraints<=m.parameters.min_thermal_constraints, + converged=1; + end + end + + count=count+1; + end + +end Index: /issm/trunk/src/m/solutions/jpl/transient.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/transient.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/transient.m (revision 2548) @@ -0,0 +1,18 @@ +function md=transient(md); +%TRANSIENT - transient solution +% +% Usage: +% md=transient(md) + +%start timing +t1=clock; + +if strcmpi(md.type,'2d'), + md=transient2d(md); +else + md=transient3d(md); +end + +%stop timing +t2=clock; +disp(sprintf('\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds'])); Index: /issm/trunk/src/m/solutions/jpl/transient2d.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/transient2d.m (revision 2548) +++ /issm/trunk/src/m/solutions/jpl/transient2d.m (revision 2548) @@ -0,0 +1,116 @@ +function md=transient2d(md); +%TRANSIENT2D - 2d transient solution +% +% Usage: +% md=transient(md) +% +% See also: TRANSIENT3D, TRANSIENT + +%Build all models requested +models.analysis_type='transient'; %needed for processresults + +displaystring(md.verbose,'%s',['reading diagnostic horiz model data']); +md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=HorizAnalysisEnum(); models.dh=CreateFemModel(md); + +displaystring(md.verbose,'\n%s',['reading diagnostic vert model data']); +md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=VertAnalysisEnum(); models.dv=CreateFemModel(md); + +displaystring(md.verbose,'\n%s',['reading diagnostic stokes model data']); +md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=StokesAnalysisEnum(); models.ds=CreateFemModel(md); + +displaystring(md.verbose,'\n%s',['reading diagnostic hutter model data']); +md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=HutterAnalysisEnum(); models.dhu=CreateFemModel(md); + +displaystring(md.verbose,'\n%s',['reading surface and bed slope computation model data']); +md.analysis_type=SlopeComputeAnalysisEnum(); md.sub_analysis_type=NoneAnalysisEnum(); models.sl=CreateFemModel(md); + +displaystring(md.verbose,'%s',['reading prognostic model data']); +md.analysis_type=PrognosticAnalysisEnum(); models.p=CreateFemModel(md); + +%initialize solution +solution=struct('step',[],'time',[],'u_g',[],'p_g',[],'h_g',[],'s_g',[],'b_g',[]); +solution.step=1; +solution.time=0; +solution.u_g=models.dh.parameters.u_g(dofsetgen([1,2],3,length(models.dh.parameters.u_g))); +solution.p_g=[]; +solution.h_g=models.p.parameters.h_g; +solution.s_g=models.p.parameters.s_g; +solution.b_g=models.p.parameters.b_g; + +%initialize inputs +displaystring(md.verbose,'\n%s',['setup inputs...']); +inputs=inputlist; +inputs=add(inputs,'velocity',solution.u_g,'doublevec',2,models.p.parameters.numberofnodes); +inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.numberofnodes); +inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.numberofnodes); +inputs=add(inputs,'dt',models.p.parameters.dt*models.p.parameters.yts,'double'); %change time from yrs to sec + +% figure out number of dof: just for information purposes. +md.dof=modelsize(models); + +%first time step is given by model. +dt=models.p.parameters.dt; +finaltime=models.p.parameters.ndt; +time=dt; +yts=models.p.parameters.yts; +n=1; %counter + +while time