Index: /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m (revision 1859) +++ /issm/trunk/src/m/solutions/cielo/objectivefunctionC.m (revision 1860) @@ -8,5 +8,5 @@ fit=m.parameters.fit(n); control_type=m.parameters.control_type; -thermalstatic=m.parameters.thermalstatic; +steadystate=m.parameters.steadystate; %Update along gradient using scalar supplied by fmincon optimization routine @@ -17,10 +17,10 @@ %Run diagnostic with updated parameters. -if ~thermalstatic, +if ~steadystate, %do a simple diagnostic, with the current temperature profile, do not look for steady state. u_g=diagnostic_core_nonlinear(m,inputs,analysis_type,sub_analysis_type); else %do a full thermal mechanical steady state converged computation, much slower! - results=thermalstatic_core(models,inputs); u_g=results.u_g; + results=steadystate_core(models,inputs); u_g=results.u_g; %u_g ships with 3 or 4 dofs, we only want the horizontal ones! Index: /issm/trunk/src/m/solutions/cielo/steadystate.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/steadystate.m (revision 1860) +++ /issm/trunk/src/m/solutions/cielo/steadystate.m (revision 1860) @@ -0,0 +1,61 @@ +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.debug,'%s',['reading diagnostic horiz model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HorizAnalysisEnum; models.dh=CreateFemModel(md); + + displaystring(md.debug,'\n%s',['reading diagnostic vert model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=VertAnalysisEnum; models.dv=CreateFemModel(md); + + displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=StokesAnalysisEnum; models.ds=CreateFemModel(md); + + displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']); + md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HutterAnalysisEnum; models.dhu=CreateFemModel(md); + + displaystring(md.debug,'\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.debug,'%s',['reading thermal model data']); + md.analysis_type=ThermalAnalysisEnum(); md.sub_analysis_type=SteadyAnalysisEnum(); models.t=CreateFemModel(md); + + displaystring(md.debug,'%s',['reading melting model data']); + md.analysis_type=MeltingAnalysisEnum(); md.sub_analysis_type=SteadyAnalysisEnum(); models.m=CreateFemModel(md); + + % figure out number of dof: just for information purposes. + md.dof=modelsize(models); + + %initialize inputs + displaystring(md.debug,'\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'); + + %compute solution + if ~models.dh.parameters.qmu_analysis, + %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); + else + %launch dakota driver for steadystate core solution + Qmu(models,inputs,models.dh.parameters); + end + + %stop timing + t2=clock; + displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: /issm/trunk/src/m/solutions/cielo/steadystate_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/steadystate_core.m (revision 1860) +++ /issm/trunk/src/m/solutions/cielo/steadystate_core.m (revision 1860) @@ -0,0 +1,89 @@ +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 +debug=m_dhu.parameters.debug; +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(debug,'%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,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(debug,'%-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 + if step>1, + converged=1; + end + 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: sm/trunk/src/m/solutions/cielo/thermalstatic.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/thermalstatic.m (revision 1859) +++ (revision ) @@ -1,61 +1,0 @@ -function md=thermalstatic(md); -%THERMALSTATIC - compute the velocity and temperature field of a model in steady state. -% -% Usage: -% md=thermalstatic(md) -% - - %timing - t1=clock; - - models.analysis_type=ThermalstaticAnalysisEnum; %needed for processresults - - %Build all models requested for diagnostic simulation - displaystring(md.debug,'%s',['reading diagnostic horiz model data']); - md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HorizAnalysisEnum; models.dh=CreateFemModel(md); - - displaystring(md.debug,'\n%s',['reading diagnostic vert model data']); - md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=VertAnalysisEnum; models.dv=CreateFemModel(md); - - displaystring(md.debug,'\n%s',['reading diagnostic stokes model data']); - md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=StokesAnalysisEnum; models.ds=CreateFemModel(md); - - displaystring(md.debug,'\n%s',['reading diagnostic hutter model data']); - md.analysis_type=DiagnosticAnalysisEnum; md.sub_analysis_type=HutterAnalysisEnum; models.dhu=CreateFemModel(md); - - displaystring(md.debug,'\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.debug,'%s',['reading thermal model data']); - md.analysis_type=ThermalAnalysisEnum(); md.sub_analysis_type=SteadyAnalysisEnum(); models.t=CreateFemModel(md); - - displaystring(md.debug,'%s',['reading melting model data']); - md.analysis_type=MeltingAnalysisEnum(); md.sub_analysis_type=SteadyAnalysisEnum(); models.m=CreateFemModel(md); - - % figure out number of dof: just for information purposes. - md.dof=modelsize(models); - - %initialize inputs - displaystring(md.debug,'\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'); - - %compute solution - if ~models.dh.parameters.qmu_analysis, - %launch core of thermalstatic solution. - results=thermalstatic_core(models,inputs); - - %process results - if ~isstruct(md.results), md.results=struct(); end - md.results.thermalstatic=processresults(models,results); - else - %launch dakota driver for thermalstatic core solution - Qmu(models,inputs,models.dh.parameters); - end - - %stop timing - t2=clock; - displaystring(md.debug,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/cielo/thermalstatic_core.m =================================================================== --- /issm/trunk/src/m/solutions/cielo/thermalstatic_core.m (revision 1859) +++ (revision ) @@ -1,87 +1,0 @@ -function results=thermalstatic_core(models,inputs); -%THERMALSTATIC_CORE - compute the core temperature and velocity field at thermal steady state. -% -% Usage: -% results=thermalstatic_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 -debug=m_dhu.parameters.debug; -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(debug,'%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,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(debug,'%-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;