Index: /issm/trunk/src/m/solutions/ControlInitialization.m =================================================================== --- /issm/trunk/src/m/solutions/ControlInitialization.m (revision 4131) +++ /issm/trunk/src/m/solutions/ControlInitialization.m (revision 4131) @@ -0,0 +1,76 @@ +function models=ControlInitialization(models); + +%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,SlopecomputeAnalysisEnum(),BedXAnalysisEnum()); +slopey=diagnostic_core_linear(m_sl,SlopecomputeAnalysisEnum(),BedYAnalysisEnum()); +slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.vertices,m_sl.loads,m_sl.materials,m_sl.parameters,slopex,'slopex',0); +slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.vertices,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,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); +displaystring(verbose,'\n%s',['extruding horizontal velocities...']); +u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.vertices,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,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,mdh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); +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=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,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/ControlOptions.m =================================================================== --- /issm/trunk/src/m/solutions/ControlOptions.m (revision 4131) +++ /issm/trunk/src/m/solutions/ControlOptions.m (revision 4131) @@ -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/FemModelUpdateInputsFromVector.m =================================================================== --- /issm/trunk/src/m/solutions/FemModelUpdateInputsFromVector.m (revision 4131) +++ /issm/trunk/src/m/solutions/FemModelUpdateInputsFromVector.m (revision 4131) @@ -0,0 +1,13 @@ +function FemModelUpdateInputsFromVector(model, vector, enum, typeenum); +%INPUT FemModelUpdateInputsFromVector(model, vector, enum, typeenum); +% +% Update inputs using a vector, just calls the UpdateInputsFromVector routine for FemModel +% +% Usage: FemModelUpdateInputsFromVector(model, vector, enum, typeenum); +% +% ex: FemModelUpdateInputsFromVector(model, vx, VxEnum, VertexEnum); +% FemModelUpdateInputsFromVector(model, vxelem, VxEnum, ElementEnum); +% +% + + Index: /issm/trunk/src/m/solutions/ModelUpdateInputsFromVector.m =================================================================== --- /issm/trunk/src/m/solutions/ModelUpdateInputsFromVector.m (revision 4131) +++ /issm/trunk/src/m/solutions/ModelUpdateInputsFromVector.m (revision 4131) @@ -0,0 +1,31 @@ +function models=ModelUpdateInputsFromVector(models, vector, enum, typeenum); +%MODELUPDATEINPUTSFROMVECTOR - Update inputs using a vector +% +% Update inputs using a vector, just calls the FemModelUpdateInputsFromVector +% routine for all models in the 'models' structure +% +% Usage: +% ModelUpdateInputsFromVector(models, vector, enum, typeenum); +% +% Example: +% ModelUpdateInputsFromVector(models, vx, VxEnum, VertexEnum); +% ModelUpdateInputsFromVector(models, vxelem, VxEnum, ElementEnum); +% +% + +%Check that vecor is not null +if isempty(vector), + return; %don't bother +end + +%go through models and call UpdateInputsFromVector +modelfields=fields(models); +for i=1:length(modelfields), + field=modelfields(i); field=field{1}; model=models.(field); + + if isstruct(model), %there is an analysis_type model + [model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters] = UpdateInputsFromVector(model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters,vector,enum, typeenum); + models.(field)=model; + end + +end Index: /issm/trunk/src/m/solutions/NewFemModel.m =================================================================== --- /issm/trunk/src/m/solutions/NewFemModel.m (revision 4131) +++ /issm/trunk/src/m/solutions/NewFemModel.m (revision 4131) @@ -0,0 +1,50 @@ +function femmodel=NewFemModel(md,solution_type,analysis_types,nummodels); +%NEWFEMMODEL - create a finite element model out of the matlab base \@model md. +% For each analysis_type contained in analysis_types, create a set of nodes, constraints +% and loads. All analyses rely on the same elements, vertices and parameters. See +% FemModel.cpp in src/c/objects for more information on the FemModel implementation in c++ +% +% Usage: +% femmodel=NewFemModel(md,solution_type,analysis_types,nummodels) +% + + + femmodel.solution_type=solution_type; + femmodel.analysis_counter=nummodels; %point to last analysis_type carried out + + %Dynamically allocate whatever is a list of length nummodels: */ + femmodel.analysis_type_list=analysis_types; + + displaystring(md.verbose,'\n reading data from model %s...',md.name); + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.constraints,femmodel.loads,femmodel.materials,femmodel.parameters]=ModelProcessor(md,solution_type,femmodel.analysis_type_list); + + %now, go through all analyses types and post-process datasets + for i=1:length(analysis_types), + + analysis_type=femmodel.analysis_type_list(i); + displaystring(md.verbose,'%s%s',' dealing with analysis type: ',EnumAsString(analysis_type)); + + displaystring(md.verbose,'%s',' generating degrees of freedofemmodel...'); + if ~isfield(femmodel,'part') [femmodel.vertices,femmodel.part,femmodel.tpart]=VerticesDof(vertices, femmodel.parameters); %do not create partition vector twice! we only have one set of vertices! + + [femmodel.nodes]=NodesDof(femmodel.nodes,femmodel.parameters); + + displaystring(md.verbose,'%s',' generating single point constraints...'); + [femmodel.nodes,femmodel.m_yg(i)]=SpcNodes(femmodel.nodes,femmodel.constraints,analysis_type); + + displaystring(md.verbose,'%s',' generating rigid body constraints...'); + [femmodel.m_Rmg(i),femmodel.nodes]=MpcNodes(femmodel.nodes,femmodel.constraints,analysis_types); + + displaystring(md.verbose,'%s',' generating node sets...'); + femmodel.m_nodesets(i)=BuildNodeSets(femmodel.nodes,analysis_type); + + displaystring(md.verbose,'%s',' reducing single point constraints vector...'); + femmodel.m_ys(i)=Reducevectorgtos(femmodel.m_yg(i).vector,femmodel.m_nodesets(i)); + + displaystring(md.verbose,'%s',' normalizing rigid body constraints matrix...'); + femmodel.m_Gmn(i)= NormalizeConstraints(femmodel.m_Rmg(i),femmodel.m_nodesets(i)); + + displaystring(md.verbose,'%s',' configuring element and loads...'); + [femmodel.elements,femmodel.loads,femmodel.nodes,femmodel.parameters] = ConfigureObjects( femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices,femmodel.materials,femmodel.parameters); + end + end Index: /issm/trunk/src/m/solutions/SetCurrentAnalysis.m =================================================================== --- /issm/trunk/src/m/solutions/SetCurrentAnalysis.m (revision 4131) +++ /issm/trunk/src/m/solutions/SetCurrentAnalysis.m (revision 4131) @@ -0,0 +1,36 @@ +function femmodel=SetCurrentAnalysis(femmodel,analysis_enum) +%SETCURRENTANALYSIS - set current analysis used to configure elements and nodes in our solutions +% +% Usage: +% femmodel=SetCurrentAnalysis(femmodel,analysis_type) +% +% Ex: +% femmodel=SetCurrentAnalysis(femmodel,DiagnosticHorizAnalysisEnum) +% + + + %first, look for analysis: + int found=0; + for i=1:length(femmodel.analysis_type_list), + if analysis_type_list(i)==analysis_enum, + found=i; + break; + end + + if (found!=1), + error('SetCurrentAnalysis error message: could not find analysis_type in list of FemModel analyses'); + end + + + %activate matrices and vectors: + femmodel.Rmg=femmodel.m_Rmg(found); + femmodel.Gmn=femmodel.m_Gmn(found); + femmodel.nodesets=femmodel.m_nodesets(found); + femmodel.yg=femmodel.m_yg(found); + femmodel.ys=femmodel.m_ys(found); + + + %Now, plug analysis_counter and analysis_type inside the parameters: + %set counter and analyse_type + femmodel.parameters.AnalysisCounter=found; + femmodel.parameters.AnalysisType=analysis_enum; Index: /issm/trunk/src/m/solutions/SetCurrentAnalysisAlias.m =================================================================== --- /issm/trunk/src/m/solutions/SetCurrentAnalysisAlias.m (revision 4131) +++ /issm/trunk/src/m/solutions/SetCurrentAnalysisAlias.m (revision 4131) @@ -0,0 +1,37 @@ +function femmodel=SetCurrentAnalysisAlias(femmodel,base_analysis_enum,real_analysis_type) +%SETCURRENTANALYSISALIAS- set current analysis used to configure elements and nodes in our solutions, +% and use another analysis_type instead. +% +% Usage: +% femmodel=SetCurrentAnalysisAlias(femmodel,base_analysis_type,real_analysis_type) +% +% Ex: +% femmodel=SetCurrentAnalysis(femmodel,DiagnosticHorizAnalysisEnum,AdjointAnalysisEnum) +% + + + %first, look for analysis: + int found=0; + for i=1:length(femmodel.analysis_type_list), + if analysis_type_list(i)==base_analysis_enum, + found=i; + break; + end + + if (found!=1), + error('SetCurrentAnalysisAlias error message: could not find analysis_type in list of FemModel analyses'); + end + + + %activate matrices and vectors: + femmodel.Rmg=femmodel.m_Rmg(found); + femmodel.Gmn=femmodel.m_Gmn(found); + femmodel.nodesets=femmodel.m_nodesets(found); + femmodel.yg=femmodel.m_yg(found); + femmodel.ys=femmodel.m_ys(found); + + + %Now, plug analysis_counter and analysis_type inside the parameters: + %set counter and analyse_type + femmodel.parameters.AnalysisCounter=found; + femmodel.parameters.AnalysisType=real_analysis_type; Index: /issm/trunk/src/m/solutions/SpawnCore.m =================================================================== --- /issm/trunk/src/m/solutions/SpawnCore.m (revision 4131) +++ /issm/trunk/src/m/solutions/SpawnCore.m (revision 4131) @@ -0,0 +1,63 @@ +function responses=SpawnCore(models,variables,variabledescriptors,counter); +%SPAWNCORE - for Qmu analysis, using Dakota. Spawn the core solution. +% +% Usage: +% responses=SpawnCore(models,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); + +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/balancedthickness.m =================================================================== --- /issm/trunk/src/m/solutions/balancedthickness.m (revision 4131) +++ /issm/trunk/src/m/solutions/balancedthickness.m (revision 4131) @@ -0,0 +1,31 @@ +function md=balancedthickness(md) +%BALANCEDTHICKNESS - balancedthickness solution sequence. +% +% Usage: +% md=balancedthickness(md) + + %timing + t1=clock; + + analysis_types=[BalancedThicknessAnalysisEnum]; + solution_type=BalancedthicknessAnalysisEnum; + + displaystring(md.verbose,'%s',['create finite element model']); + femmodel=NewFemModel(md,solution_type,analysis_types,1); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %compute solution + if ~qmu_analysis, + displaystring(verbose,'%s',['call computational core']); + femmodel=balancedthickness_core(femmodel); + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/balancedthickness2.m =================================================================== --- /issm/trunk/src/m/solutions/balancedthickness2.m (revision 4131) +++ /issm/trunk/src/m/solutions/balancedthickness2.m (revision 4131) @@ -0,0 +1,28 @@ +function md=balancedthickness2(md) +%BALANCEDTHICKNESS - balancedthickness2 solution sequence. +% +% Usage: +% md=balancedthickness2(md) + + analysis_types=[Balancedthickness2AnalysisEnum]; + solution_type=Balancedthickness2AnalysisEnum; + + displaystring(md.verbose,'%s',['create finite element model']); + femmodel=NewFemModel(md,solution_type,analysis_types,1); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %compute solution + if ~qmu_analysis, + displaystring(verbose,'%s',['call computational core']); + femmodel=balancedthickness2_core(femmodel); + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/balancedthickness2_core.m =================================================================== --- /issm/trunk/src/m/solutions/balancedthickness2_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/balancedthickness2_core.m (revision 4131) @@ -0,0 +1,30 @@ +function femmodel=balancedthickness2_core(femmodel) +%BALANCEDTHICKNESS_CORE - linear solution sequence +% +% Usage: +% femmodel=balancedthickness2_core(femmodel) + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + dim=femmodel.parameters.Dim; + + %Activate formulation + femmodel=SetCurrentAnalysis(femmodel,Balancedthickness2AnalysisEnum); + + displaystring(verbose,'\n%s',['depth averaging velocities...']); + %[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum); + %[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum); + if dim==3, + % [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,VzAverageEnum); + end + + displaystring(verbose,'\n%s',['call computational core...']); + femmodel=solver_linear(femmodel); + + displaystring(verbose,'\n%s',['extude computed thickness on all layers...']); + %femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + +end %end function Index: /issm/trunk/src/m/solutions/balancedthickness_core.m =================================================================== --- /issm/trunk/src/m/solutions/balancedthickness_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/balancedthickness_core.m (revision 4131) @@ -0,0 +1,30 @@ +function femmodel=balancedthickness_core(femmodel) +%BALANCEDTHICKNESS_CORE - linear solution sequence +% +% Usage: +% femmodel=balancedthickness_core(femmode) + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + dim=femmodel.parameters.Dim; + + %Activate formulation + femmodel=SetCurrentAnalysis(femmodel,BalancedthicknessAnalysisEnum); + + displaystring(verbose,'\n%s',['depth averaging velocities...']); + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum); + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum); + if dim==3, + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,VzAverageEnum); + end + + displaystring(verbose,'\n%s',['call computational core...']); + femmodel=solver_linear(femmodel); + + displaystring(verbose,'\n%s',['extude computed thickness on all layers...']); + femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + +end %end function Index: /issm/trunk/src/m/solutions/balancedvelocities.m =================================================================== --- /issm/trunk/src/m/solutions/balancedvelocities.m (revision 4131) +++ /issm/trunk/src/m/solutions/balancedvelocities.m (revision 4131) @@ -0,0 +1,31 @@ +function md=balancedvelocities(md) +%BALANCEDVELOCITIES - balancedvelocities solution sequence. +% +% Usage: +% md=balancedvelocities(md) + + %timing + t1=clock; + + analysis_types=[BalancedvelocitiesAnalysisEnum]; + solution_type=BalancedvelocitiesAnalysisEnum; + + displaystring(md.verbose,'%s',['create finite element model']); + femmodel=NewFemModel(md,solution_type,analysis_types,1); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %compute solution + if ~qmu_analysis, + displaystring(verbose,'%s',['call computational core']); + femmodel=balancedvelocities_core(femmodel); + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/balancedvelocities_core.m =================================================================== --- /issm/trunk/src/m/solutions/balancedvelocities_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/balancedvelocities_core.m (revision 4131) @@ -0,0 +1,32 @@ +function femmodel=balancedvelocities_core(femmdoel) +%BALANCEDVELOCITIES_CORE - linear solution sequence +% +% Usage: +% femmodel=balancedvelocities_core(femmodel) + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + dim=femmodel.parameters.Dim; + + %Activate formulation + femmodel=SetCurrentAnalysis(femmodel,BalancedvelocitiesAnalysisEnum); + + displaystring(verbose,'\n%s',['depth averaging velocities...']); + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum); + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum); + if dim==3, + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,VzAverageEnum); + end + + displaystring(verbose,'\n%s',['call computational core...']); + femmodel=solver_linear(femmodel); + + displaystring(verbose,'\n%s',['extude computed thickness on all layers...']); + femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + +end %end function Index: /issm/trunk/src/m/solutions/bedslope.m =================================================================== --- /issm/trunk/src/m/solutions/bedslope.m (revision 4131) +++ /issm/trunk/src/m/solutions/bedslope.m (revision 4131) @@ -0,0 +1,35 @@ +function md=bedslopecompute(md); +%SLOPECOMPUTE - compute the bed slope of a model +% +% Usage: +% md=bedslope(md) +% + %timing + t1=clock; + + analysis_types=SlopeAnalysisEnum; + solution_type=SlopeAnalysisEnum; + + displaystring(md.verbose,'%s',['create finite element model']); + femmodel=NewFemModel(md,solution_type,analysis_types,1); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %compute solution + if ~qmu_analysis, + displaystring(verbose,'%s',['call computational core']); + femmodel=bedslope_core(femmodel); + + displaystring(verbose,'%s',['write results']); + OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters); + + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/bedslope_core.m =================================================================== --- /issm/trunk/src/m/solutions/bedslope_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/bedslope_core.m (revision 4131) @@ -0,0 +1,30 @@ +function femmodel=bedslope_core(femmodel) +%SURFACESLOPE_CORE - core of the bed slope computation solution +% +% Usage: +% femmodel=bedslope_core(femmodel) +% + + + %Recover some parameters: + verbose=femmodel.parameters.Verbose; + dim=femmodel.parameters.Dim; + + displaystring(verbose,'\n%s',['computing bed slope...']); + + %Call on core computations: + femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum); + femmodel=solver_linear(femmodel); + femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum); + femmodel=solver_linear(femmodel); + + %extrude inputs if we are in 3D: */ + if dim==3, + displaystring(verbose,'\n%s',['extruding bed slopein 3d...']); + femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum); + femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum); + } + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum); + InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum); Index: /issm/trunk/src/m/solutions/control_core.m =================================================================== --- /issm/trunk/src/m/solutions/control_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/control_core.m (revision 4131) @@ -0,0 +1,148 @@ +function femmodel=control_core(femmodel) +%CONTROL_CORE - compute the core inversion +% +% Usage: +% femmodel=control_core(femmodel); +% + + %Preprocess models + femmodel=ControlInitialization(femmodel); + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + nsteps=femmodel.parameters.NSteps; + control_type=femmodel.parameters.ControlType; + fit=femmodel.parameters.Fit; + optscal=femmodel.parameters.OptScal; + maxiter=femmodel.parameters.MaxIter; + cm_jump=femmodel.parameters.CmJump; + eps_cm=femmodel.parameters.EpsCm; + tolx=femmodel.parameters.TolX; + cm_min=femmodel.parameters.CmMin; + cm_max=femmodel.parameters.CmMax; + cm_gradient=femmodel.parameters.CmGradient; + control_steady=femmodel.parameters.ControlSteady; + +%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.ControlSteady; + steadystate_results=steadystate_core(models); t_g=steadystate_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.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes); + + %Update inputs in datasets + [model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters]=UpdateFromInputs(model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters); + + displaystring(verbose,'\n%s',[' computing gradJ...']); + results_grad=gradjcompute_core(models); + 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',4,model.parameters.NumberOfNodes); + else + if model.parameters.ControlSteady; + inputs=add(inputs,'velocity',u_g,'doublevec',3,model.parameters.NumberOfNodes); + else + inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.NumberOfNodes); + end + end + else + inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.NumberOfNodes); + end + + if n>=2 & search_scalar==0, + displaystring(verbose,'\n%s',[' normalizing directions...']); + 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,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.ControlSteady; + inputs=add(inputs,model.parameters.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes); + steadystate_results=steadystate_core(models); t_g=steadystate_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.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes); + results_diag=diagnostic_core(models); + 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.ControlSteady, + results.t_g=t_g; + results.m_g=m_g; +end Index: /issm/trunk/src/m/solutions/convergence.m =================================================================== --- /issm/trunk/src/m/solutions/convergence.m (revision 4131) +++ /issm/trunk/src/m/solutions/convergence.m (revision 4131) @@ -0,0 +1,87 @@ +function converged=convergence(K_ff,p_f,u_f,u_f_old,parameters) + +%Get convergence options +verbose=parameters.Verbose; +yts=parameters.Yts; +eps_res=parameters.EpsRes; +eps_rel=parameters.EpsRel; +eps_abs=parameters.EpsAbs; + +%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 nu, %avoid "dividing by zero" warning + 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 + converged=0; + end + else + if nu, %avoid "dividing by zero" warning + displaystring(verbose,'%-60s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' %'); + end + 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/diagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/diagnostic.m (revision 4131) +++ /issm/trunk/src/m/solutions/diagnostic.m (revision 4131) @@ -0,0 +1,45 @@ +function md=diagnostic(md); +%DIAGNOSTIC - compute the velocity field of a model +% +% Usage: +% md=diagnostic(md) +% + %timing + t1=clock; + + analysis_types=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum]; + solution_type=DiagnosticAnalysisEnum; + + displaystring(md.verbose,'%s',['create finite element model']); + femmodel=NewFemModel(md,solution_type,analysis_types,5); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + control_analysis=femmodel.parameters.ControlAnalysis; + + %compute solution + if ~qmu_analysis, + if ~control_analysis, + + displaystring(verbose,'%s',['call computational core']); + femmodel=diagnostic_core(femmodel); + + else, + + displaystring(verbose,'%s',['call computational core']); + femmodel=control_core(femmodel); + + end + + displaystring(verbose,'%s',['write results']); + OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters); + + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/diagnostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/diagnostic_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/diagnostic_core.m (revision 4131) @@ -0,0 +1,76 @@ +function femmodel=diagnostic_core(femmodel); +%DIAGNOSTIC_CORE - compute the core velocity field +% +% Usage: +% results=diagnostic_core(model); +% + + %some parameters + modify_loads=boolean(1); + conserve_loads=boolean(1); + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + dim=femmodel.parameters.Dim; + ishutter=femmodel.parameters.IsHutter; + ismacayealpattyn=femmodel.parameters.IsMacAyealPattyn; + isstokes=femmodel.parameters.IsStokes; + stokesreconditioning=femmodel.parameters.stokesreconditioning; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %for qmu analysis, be sure the velocity input we are starting from is the one in the parameters: + if qmu_analysis, + femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVxEnum,VxEnum); + femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVyEnum,VyEnum); + femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVzEnum,VzEnum); + femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuPressureEnum,PressureEnum); + end + + %Compute slopes: + if(ishutter)femmodel=surfaceslope_core(femmodel);end + if(isstokes)femmodel=bedslope_core(femmodel);end + + if ishutter, + + displaystring(verbose,'\n%s',['computing hutter velocities...']); + femmodel=SetCurrentAnalysis(femmodel,DiagnosticHutterAnalysisEnum); + femmodel=solver_linear(femmodel); + + if(ismacayealpattyn)femmodel=ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,HorizAnalysisEnum); end + + end + + if ismacayealpattyn, + + displaystring(verbose,'\n%s',['computing horizontal velocities...']); + femmodel=SetCurrentAnalysis(femmodel,DiagnosticHorizAnalysisEnum); + femmodel=solver_diagnostic_nonlinear(femmodel,modify_loads); + end + + + if dim==3, + + displaystring(verbose,'\n%s',['computing vertical velocities...']); + femmodel=SetCurrentAnalysis(femmodel,DiagnosticVertAnalysisEnum); + femmodel=solver_linear(femmodel); + + if isstokes, + + %"recondition" pressure computed previously: + femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum); + InputScale(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning); + + displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']); + femmodel=ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum); + + displaystring(verbose,'\n%s',['computing stokes velocities and pressure...']); + femmodel=SetCurrentAnalysis(femmodel,DiagnosticStokesAnalysisEnum); + femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads); + end + end + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); + if(dim==3) InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); Index: /issm/trunk/src/m/solutions/dofsetgen.m =================================================================== --- /issm/trunk/src/m/solutions/dofsetgen.m (revision 4131) +++ /issm/trunk/src/m/solutions/dofsetgen.m (revision 4131) @@ -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/gradjcompute_core.m =================================================================== --- /issm/trunk/src/m/solutions/gradjcompute_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/gradjcompute_core.m (revision 4131) @@ -0,0 +1,52 @@ +function results=gradjcompute_core(models); + +%recover active models +m=models.active; + +%recover parameters +verbose=m.parameters.Verbose; +dim=m.parameters.Dim; +extrude_param=m.parameters.ExtrudeParam; +ishutter=m.parameters.IsHutter; +ismacayealpattyn=m.parameters.IsMacAyealPattyn; +isstokes=m.parameters.IsStokes; +analysis_type=m.parameters.AnalysisType; +sub_analysis_type=m.parameters.SubAnalysisType; + +%Recover solution for this stiffness and right hand side: +displaystring(verbose,'%s',' computing velocities...'); +[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,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.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); + +%Reduce adjoint load from g-set to f-set +[Du_f] = Reduceloadfromgtof( Du_g, m.Gmn, K_fs0, m.ys, m.nodesets,1);% 1 because we want ys0 + +%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.ys, m.nodesets,1); % 1 because we want ys0 +inputs=add(inputs,'adjoint',lambda_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes); + +%Compute gradJ +grad_g=Gradj(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); +if (dim==3 & extrude_param), + displaystring(verbose,'%s',' extruding gradient...'); + grad_g=FieldExtrude(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,grad_g,'gradj',0); +end + +if m.parameters.ControlSteady, + %compute initial velocity from diagnostic_core (horiz+vertical) + displaystring(verbose,'%s',' compute 3d initial velocity...'); + results_diag=diagnostic_core(models); + u_g=results_diag.u_g; +end + +%Assign output +results.u_g=u_g; +results.grad_g=grad_g; Index: sm/trunk/src/m/solutions/jpl/ControlInitialization.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/ControlInitialization.m (revision 4130) +++ (revision ) @@ -1,76 +1,0 @@ -function models=ControlInitialization(models); - -%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,SlopecomputeAnalysisEnum(),BedXAnalysisEnum()); -slopey=diagnostic_core_linear(m_sl,SlopecomputeAnalysisEnum(),BedYAnalysisEnum()); -slopex=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.vertices,m_sl.loads,m_sl.materials,m_sl.parameters,slopex,'slopex',0); -slopey=FieldExtrude(m_sl.elements,m_sl.nodes,m_sl.vertices,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,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); -displaystring(verbose,'\n%s',['extruding horizontal velocities...']); -u_g_horiz=FieldExtrude(m_dh.elements,m_dh.nodes,m_dh.vertices,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,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,mdh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum()); -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=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,DiagnosticAnalysisEnum(),StokesAnalysisEnum()); -inputs=add(inputs,'velocity',u_g,'doublevec',4,m_ds.parameters.numberofnodes); - -%assign output -models.active=m_ds; Index: sm/trunk/src/m/solutions/jpl/ControlOptions.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/ControlOptions.m (revision 4130) +++ (revision ) @@ -1,13 +1,0 @@ -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: sm/trunk/src/m/solutions/jpl/FemModelUpdateInputsFromVector.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/FemModelUpdateInputsFromVector.m (revision 4130) +++ (revision ) @@ -1,13 +1,0 @@ -function FemModelUpdateInputsFromVector(model, vector, enum, typeenum); -%INPUT FemModelUpdateInputsFromVector(model, vector, enum, typeenum); -% -% Update inputs using a vector, just calls the UpdateInputsFromVector routine for FemModel -% -% Usage: FemModelUpdateInputsFromVector(model, vector, enum, typeenum); -% -% ex: FemModelUpdateInputsFromVector(model, vx, VxEnum, VertexEnum); -% FemModelUpdateInputsFromVector(model, vxelem, VxEnum, ElementEnum); -% -% - - Index: sm/trunk/src/m/solutions/jpl/ModelUpdateInputsFromVector.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/ModelUpdateInputsFromVector.m (revision 4130) +++ (revision ) @@ -1,31 +1,0 @@ -function models=ModelUpdateInputsFromVector(models, vector, enum, typeenum); -%MODELUPDATEINPUTSFROMVECTOR - Update inputs using a vector -% -% Update inputs using a vector, just calls the FemModelUpdateInputsFromVector -% routine for all models in the 'models' structure -% -% Usage: -% ModelUpdateInputsFromVector(models, vector, enum, typeenum); -% -% Example: -% ModelUpdateInputsFromVector(models, vx, VxEnum, VertexEnum); -% ModelUpdateInputsFromVector(models, vxelem, VxEnum, ElementEnum); -% -% - -%Check that vecor is not null -if isempty(vector), - return; %don't bother -end - -%go through models and call UpdateInputsFromVector -modelfields=fields(models); -for i=1:length(modelfields), - field=modelfields(i); field=field{1}; model=models.(field); - - if isstruct(model), %there is an analysis_type model - [model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters] = UpdateInputsFromVector(model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters,vector,enum, typeenum); - models.(field)=model; - end - -end Index: sm/trunk/src/m/solutions/jpl/NewFemModel.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/NewFemModel.m (revision 4130) +++ (revision ) @@ -1,50 +1,0 @@ -function femmodel=NewFemModel(md,solution_type,analysis_types,nummodels); -%NEWFEMMODEL - create a finite element model out of the matlab base \@model md. -% For each analysis_type contained in analysis_types, create a set of nodes, constraints -% and loads. All analyses rely on the same elements, vertices and parameters. See -% FemModel.cpp in src/c/objects for more information on the FemModel implementation in c++ -% -% Usage: -% femmodel=NewFemModel(md,solution_type,analysis_types,nummodels) -% - - - femmodel.solution_type=solution_type; - femmodel.analysis_counter=nummodels; %point to last analysis_type carried out - - %Dynamically allocate whatever is a list of length nummodels: */ - femmodel.analysis_type_list=analysis_types; - - displaystring(md.verbose,'\n reading data from model %s...',md.name); - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.constraints,femmodel.loads,femmodel.materials,femmodel.parameters]=ModelProcessor(md,solution_type,femmodel.analysis_type_list); - - %now, go through all analyses types and post-process datasets - for i=1:length(analysis_types), - - analysis_type=femmodel.analysis_type_list(i); - displaystring(md.verbose,'%s%s',' dealing with analysis type: ',EnumAsString(analysis_type)); - - displaystring(md.verbose,'%s',' generating degrees of freedofemmodel...'); - if ~isfield(femmodel,'part') [femmodel.vertices,femmodel.part,femmodel.tpart]=VerticesDof(vertices, femmodel.parameters); %do not create partition vector twice! we only have one set of vertices! - - [femmodel.nodes]=NodesDof(femmodel.nodes,femmodel.parameters); - - displaystring(md.verbose,'%s',' generating single point constraints...'); - [femmodel.nodes,femmodel.m_yg(i)]=SpcNodes(femmodel.nodes,femmodel.constraints,analysis_type); - - displaystring(md.verbose,'%s',' generating rigid body constraints...'); - [femmodel.m_Rmg(i),femmodel.nodes]=MpcNodes(femmodel.nodes,femmodel.constraints,analysis_types); - - displaystring(md.verbose,'%s',' generating node sets...'); - femmodel.m_nodesets(i)=BuildNodeSets(femmodel.nodes,analysis_type); - - displaystring(md.verbose,'%s',' reducing single point constraints vector...'); - femmodel.m_ys(i)=Reducevectorgtos(femmodel.m_yg(i).vector,femmodel.m_nodesets(i)); - - displaystring(md.verbose,'%s',' normalizing rigid body constraints matrix...'); - femmodel.m_Gmn(i)= NormalizeConstraints(femmodel.m_Rmg(i),femmodel.m_nodesets(i)); - - displaystring(md.verbose,'%s',' configuring element and loads...'); - [femmodel.elements,femmodel.loads,femmodel.nodes,femmodel.parameters] = ConfigureObjects( femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices,femmodel.materials,femmodel.parameters); - end - end Index: sm/trunk/src/m/solutions/jpl/SetCurrentAnalysis.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/SetCurrentAnalysis.m (revision 4130) +++ (revision ) @@ -1,36 +1,0 @@ -function femmodel=SetCurrentAnalysis(femmodel,analysis_enum) -%SETCURRENTANALYSIS - set current analysis used to configure elements and nodes in our solutions -% -% Usage: -% femmodel=SetCurrentAnalysis(femmodel,analysis_type) -% -% Ex: -% femmodel=SetCurrentAnalysis(femmodel,DiagnosticHorizAnalysisEnum) -% - - - %first, look for analysis: - int found=0; - for i=1:length(femmodel.analysis_type_list), - if analysis_type_list(i)==analysis_enum, - found=i; - break; - end - - if (found!=1), - error('SetCurrentAnalysis error message: could not find analysis_type in list of FemModel analyses'); - end - - - %activate matrices and vectors: - femmodel.Rmg=femmodel.m_Rmg(found); - femmodel.Gmn=femmodel.m_Gmn(found); - femmodel.nodesets=femmodel.m_nodesets(found); - femmodel.yg=femmodel.m_yg(found); - femmodel.ys=femmodel.m_ys(found); - - - %Now, plug analysis_counter and analysis_type inside the parameters: - %set counter and analyse_type - femmodel.parameters.AnalysisCounter=found; - femmodel.parameters.AnalysisType=analysis_enum; Index: sm/trunk/src/m/solutions/jpl/SetCurrentAnalysisAlias.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/SetCurrentAnalysisAlias.m (revision 4130) +++ (revision ) @@ -1,37 +1,0 @@ -function femmodel=SetCurrentAnalysisAlias(femmodel,base_analysis_enum,real_analysis_type) -%SETCURRENTANALYSISALIAS- set current analysis used to configure elements and nodes in our solutions, -% and use another analysis_type instead. -% -% Usage: -% femmodel=SetCurrentAnalysisAlias(femmodel,base_analysis_type,real_analysis_type) -% -% Ex: -% femmodel=SetCurrentAnalysis(femmodel,DiagnosticHorizAnalysisEnum,AdjointAnalysisEnum) -% - - - %first, look for analysis: - int found=0; - for i=1:length(femmodel.analysis_type_list), - if analysis_type_list(i)==base_analysis_enum, - found=i; - break; - end - - if (found!=1), - error('SetCurrentAnalysisAlias error message: could not find analysis_type in list of FemModel analyses'); - end - - - %activate matrices and vectors: - femmodel.Rmg=femmodel.m_Rmg(found); - femmodel.Gmn=femmodel.m_Gmn(found); - femmodel.nodesets=femmodel.m_nodesets(found); - femmodel.yg=femmodel.m_yg(found); - femmodel.ys=femmodel.m_ys(found); - - - %Now, plug analysis_counter and analysis_type inside the parameters: - %set counter and analyse_type - femmodel.parameters.AnalysisCounter=found; - femmodel.parameters.AnalysisType=real_analysis_type; Index: sm/trunk/src/m/solutions/jpl/SpawnCore.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/SpawnCore.m (revision 4130) +++ (revision ) @@ -1,63 +1,0 @@ -function responses=SpawnCore(models,variables,variabledescriptors,counter); -%SPAWNCORE - for Qmu analysis, using Dakota. Spawn the core solution. -% -% Usage: -% responses=SpawnCore(models,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); - -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: sm/trunk/src/m/solutions/jpl/balancedthickness.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/balancedthickness.m (revision 4130) +++ (revision ) @@ -1,31 +1,0 @@ -function md=balancedthickness(md) -%BALANCEDTHICKNESS - balancedthickness solution sequence. -% -% Usage: -% md=balancedthickness(md) - - %timing - t1=clock; - - analysis_types=[BalancedThicknessAnalysisEnum]; - solution_type=BalancedthicknessAnalysisEnum; - - displaystring(md.verbose,'%s',['create finite element model']); - femmodel=NewFemModel(md,solution_type,analysis_types,1); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %compute solution - if ~qmu_analysis, - displaystring(verbose,'%s',['call computational core']); - femmodel=balancedthickness_core(femmodel); - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/balancedthickness2.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/balancedthickness2.m (revision 4130) +++ (revision ) @@ -1,28 +1,0 @@ -function md=balancedthickness2(md) -%BALANCEDTHICKNESS - balancedthickness2 solution sequence. -% -% Usage: -% md=balancedthickness2(md) - - analysis_types=[Balancedthickness2AnalysisEnum]; - solution_type=Balancedthickness2AnalysisEnum; - - displaystring(md.verbose,'%s',['create finite element model']); - femmodel=NewFemModel(md,solution_type,analysis_types,1); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %compute solution - if ~qmu_analysis, - displaystring(verbose,'%s',['call computational core']); - femmodel=balancedthickness2_core(femmodel); - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/balancedthickness2_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/balancedthickness2_core.m (revision 4130) +++ (revision ) @@ -1,30 +1,0 @@ -function femmodel=balancedthickness2_core(femmodel) -%BALANCEDTHICKNESS_CORE - linear solution sequence -% -% Usage: -% femmodel=balancedthickness2_core(femmodel) - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - dim=femmodel.parameters.Dim; - - %Activate formulation - femmodel=SetCurrentAnalysis(femmodel,Balancedthickness2AnalysisEnum); - - displaystring(verbose,'\n%s',['depth averaging velocities...']); - %[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum); - %[femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum); - if dim==3, - % [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,VzAverageEnum); - end - - displaystring(verbose,'\n%s',['call computational core...']); - femmodel=solver_linear(femmodel); - - displaystring(verbose,'\n%s',['extude computed thickness on all layers...']); - %femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); - -end %end function Index: sm/trunk/src/m/solutions/jpl/balancedthickness_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/balancedthickness_core.m (revision 4130) +++ (revision ) @@ -1,30 +1,0 @@ -function femmodel=balancedthickness_core(femmodel) -%BALANCEDTHICKNESS_CORE - linear solution sequence -% -% Usage: -% femmodel=balancedthickness_core(femmode) - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - dim=femmodel.parameters.Dim; - - %Activate formulation - femmodel=SetCurrentAnalysis(femmodel,BalancedthicknessAnalysisEnum); - - displaystring(verbose,'\n%s',['depth averaging velocities...']); - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum); - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum); - if dim==3, - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,VzAverageEnum); - end - - displaystring(verbose,'\n%s',['call computational core...']); - femmodel=solver_linear(femmodel); - - displaystring(verbose,'\n%s',['extude computed thickness on all layers...']); - femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); - -end %end function Index: sm/trunk/src/m/solutions/jpl/balancedvelocities.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/balancedvelocities.m (revision 4130) +++ (revision ) @@ -1,31 +1,0 @@ -function md=balancedvelocities(md) -%BALANCEDVELOCITIES - balancedvelocities solution sequence. -% -% Usage: -% md=balancedvelocities(md) - - %timing - t1=clock; - - analysis_types=[BalancedvelocitiesAnalysisEnum]; - solution_type=BalancedvelocitiesAnalysisEnum; - - displaystring(md.verbose,'%s',['create finite element model']); - femmodel=NewFemModel(md,solution_type,analysis_types,1); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %compute solution - if ~qmu_analysis, - displaystring(verbose,'%s',['call computational core']); - femmodel=balancedvelocities_core(femmodel); - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/balancedvelocities_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m (revision 4130) +++ (revision ) @@ -1,32 +1,0 @@ -function femmodel=balancedvelocities_core(femmdoel) -%BALANCEDVELOCITIES_CORE - linear solution sequence -% -% Usage: -% femmodel=balancedvelocities_core(femmodel) - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - dim=femmodel.parameters.Dim; - - %Activate formulation - femmodel=SetCurrentAnalysis(femmodel,BalancedvelocitiesAnalysisEnum); - - displaystring(verbose,'\n%s',['depth averaging velocities...']); - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum); - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum); - if dim==3, - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VzEnum,VzAverageEnum); - end - - displaystring(verbose,'\n%s',['call computational core...']); - femmodel=solver_linear(femmodel); - - displaystring(verbose,'\n%s',['extude computed thickness on all layers...']); - femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); - femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); - -end %end function Index: sm/trunk/src/m/solutions/jpl/bedslope.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/bedslope.m (revision 4130) +++ (revision ) @@ -1,35 +1,0 @@ -function md=bedslopecompute(md); -%SLOPECOMPUTE - compute the bed slope of a model -% -% Usage: -% md=bedslope(md) -% - %timing - t1=clock; - - analysis_types=SlopeAnalysisEnum; - solution_type=SlopeAnalysisEnum; - - displaystring(md.verbose,'%s',['create finite element model']); - femmodel=NewFemModel(md,solution_type,analysis_types,1); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %compute solution - if ~qmu_analysis, - displaystring(verbose,'%s',['call computational core']); - femmodel=bedslope_core(femmodel); - - displaystring(verbose,'%s',['write results']); - OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters); - - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/bedslope_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/bedslope_core.m (revision 4130) +++ (revision ) @@ -1,30 +1,0 @@ -function femmodel=bedslope_core(femmodel) -%SURFACESLOPE_CORE - core of the bed slope computation solution -% -% Usage: -% femmodel=bedslope_core(femmodel) -% - - - %Recover some parameters: - verbose=femmodel.parameters.Verbose; - dim=femmodel.parameters.Dim; - - displaystring(verbose,'\n%s',['computing bed slope...']); - - %Call on core computations: - femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum); - femmodel=solver_linear(femmodel); - femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum); - femmodel=solver_linear(femmodel); - - %extrude inputs if we are in 3D: */ - if dim==3, - displaystring(verbose,'\n%s',['extruding bed slopein 3d...']); - femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum); - femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum); - } - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum); - InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum); Index: sm/trunk/src/m/solutions/jpl/control_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/control_core.m (revision 4130) +++ (revision ) @@ -1,148 +1,0 @@ -function femmodel=control_core(femmodel) -%CONTROL_CORE - compute the core inversion -% -% Usage: -% femmodel=control_core(femmodel); -% - - %Preprocess models - femmodel=ControlInitialization(femmodel); - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - nsteps=femmodel.parameters.NSteps; - control_type=femmodel.parameters.ControlType; - fit=femmodel.parameters.Fit; - optscal=femmodel.parameters.OptScal; - maxiter=femmodel.parameters.MaxIter; - cm_jump=femmodel.parameters.CmJump; - eps_cm=femmodel.parameters.EpsCm; - tolx=femmodel.parameters.TolX; - cm_min=femmodel.parameters.CmMin; - cm_max=femmodel.parameters.CmMax; - cm_gradient=femmodel.parameters.CmGradient; - control_steady=femmodel.parameters.ControlSteady; - -%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.ControlSteady; - steadystate_results=steadystate_core(models); t_g=steadystate_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.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes); - - %Update inputs in datasets - [model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters]=UpdateFromInputs(model.elements,model.nodes,model.vertices,model.loads,model.materials,model.parameters); - - displaystring(verbose,'\n%s',[' computing gradJ...']); - results_grad=gradjcompute_core(models); - 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',4,model.parameters.NumberOfNodes); - else - if model.parameters.ControlSteady; - inputs=add(inputs,'velocity',u_g,'doublevec',3,model.parameters.NumberOfNodes); - else - inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.NumberOfNodes); - end - end - else - inputs=add(inputs,'velocity',u_g,'doublevec',2,model.parameters.NumberOfNodes); - end - - if n>=2 & search_scalar==0, - displaystring(verbose,'\n%s',[' normalizing directions...']); - 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,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.ControlSteady; - inputs=add(inputs,model.parameters.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes); - steadystate_results=steadystate_core(models); t_g=steadystate_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.ControlType,param_g,'doublevec',1,model.parameters.NumberOfNodes); - results_diag=diagnostic_core(models); - 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.ControlSteady, - results.t_g=t_g; - results.m_g=m_g; -end Index: sm/trunk/src/m/solutions/jpl/convergence.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/convergence.m (revision 4130) +++ (revision ) @@ -1,87 +1,0 @@ -function converged=convergence(K_ff,p_f,u_f,u_f_old,parameters) - -%Get convergence options -verbose=parameters.Verbose; -yts=parameters.Yts; -eps_res=parameters.EpsRes; -eps_rel=parameters.EpsRel; -eps_abs=parameters.EpsAbs; - -%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 nu, %avoid "dividing by zero" warning - 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 - converged=0; - end - else - if nu, %avoid "dividing by zero" warning - displaystring(verbose,'%-60s%g%s',' relative convergence criterion: norm(du)/norm(u)',ndu/nu*100,' %'); - end - 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: sm/trunk/src/m/solutions/jpl/diagnostic.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic.m (revision 4130) +++ (revision ) @@ -1,45 +1,0 @@ -function md=diagnostic(md); -%DIAGNOSTIC - compute the velocity field of a model -% -% Usage: -% md=diagnostic(md) -% - %timing - t1=clock; - - analysis_types=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum]; - solution_type=DiagnosticAnalysisEnum; - - displaystring(md.verbose,'%s',['create finite element model']); - femmodel=NewFemModel(md,solution_type,analysis_types,5); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - control_analysis=femmodel.parameters.ControlAnalysis; - - %compute solution - if ~qmu_analysis, - if ~control_analysis, - - displaystring(verbose,'%s',['call computational core']); - femmodel=diagnostic_core(femmodel); - - else, - - displaystring(verbose,'%s',['call computational core']); - femmodel=control_core(femmodel); - - end - - displaystring(verbose,'%s',['write results']); - OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters); - - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/diagnostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/diagnostic_core.m (revision 4130) +++ (revision ) @@ -1,76 +1,0 @@ -function femmodel=diagnostic_core(femmodel); -%DIAGNOSTIC_CORE - compute the core velocity field -% -% Usage: -% results=diagnostic_core(model); -% - - %some parameters - modify_loads=boolean(1); - conserve_loads=boolean(1); - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - dim=femmodel.parameters.Dim; - ishutter=femmodel.parameters.IsHutter; - ismacayealpattyn=femmodel.parameters.IsMacAyealPattyn; - isstokes=femmodel.parameters.IsStokes; - stokesreconditioning=femmodel.parameters.stokesreconditioning; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %for qmu analysis, be sure the velocity input we are starting from is the one in the parameters: - if qmu_analysis, - femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVxEnum,VxEnum); - femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVyEnum,VyEnum); - femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVzEnum,VzEnum); - femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuPressureEnum,PressureEnum); - end - - %Compute slopes: - if(ishutter)femmodel=surfaceslope_core(femmodel);end - if(isstokes)femmodel=bedslope_core(femmodel);end - - if ishutter, - - displaystring(verbose,'\n%s',['computing hutter velocities...']); - femmodel=SetCurrentAnalysis(femmodel,DiagnosticHutterAnalysisEnum); - femmodel=solver_linear(femmodel); - - if(ismacayealpattyn)femmodel=ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,HorizAnalysisEnum); end - - end - - if ismacayealpattyn, - - displaystring(verbose,'\n%s',['computing horizontal velocities...']); - femmodel=SetCurrentAnalysis(femmodel,DiagnosticHorizAnalysisEnum); - femmodel=solver_diagnostic_nonlinear(femmodel,modify_loads); - end - - - if dim==3, - - displaystring(verbose,'\n%s',['computing vertical velocities...']); - femmodel=SetCurrentAnalysis(femmodel,DiagnosticVertAnalysisEnum); - femmodel=solver_linear(femmodel); - - if isstokes, - - %"recondition" pressure computed previously: - femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum); - InputScale(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning); - - displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']); - femmodel=ResetBoundaryConditions(femmodel,DiagnosticAnalysisEnum,StokesAnalysisEnum); - - displaystring(verbose,'\n%s',['computing stokes velocities and pressure...']); - femmodel=SetCurrentAnalysis(femmodel,DiagnosticStokesAnalysisEnum); - femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads); - end - end - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); - if(dim==3) InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); Index: sm/trunk/src/m/solutions/jpl/dofsetgen.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/dofsetgen.m (revision 4130) +++ (revision ) @@ -1,19 +1,0 @@ -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: sm/trunk/src/m/solutions/jpl/gradjcompute_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/gradjcompute_core.m (revision 4130) +++ (revision ) @@ -1,52 +1,0 @@ -function results=gradjcompute_core(models); - -%recover active models -m=models.active; - -%recover parameters -verbose=m.parameters.Verbose; -dim=m.parameters.Dim; -extrude_param=m.parameters.ExtrudeParam; -ishutter=m.parameters.IsHutter; -ismacayealpattyn=m.parameters.IsMacAyealPattyn; -isstokes=m.parameters.IsStokes; -analysis_type=m.parameters.AnalysisType; -sub_analysis_type=m.parameters.SubAnalysisType; - -%Recover solution for this stiffness and right hand side: -displaystring(verbose,'%s',' computing velocities...'); -[u_g K_ff0 K_fs0 ]=diagnostic_core_nonlinear(m,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.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); - -%Reduce adjoint load from g-set to f-set -[Du_f] = Reduceloadfromgtof( Du_g, m.Gmn, K_fs0, m.ys, m.nodesets,1);% 1 because we want ys0 - -%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.ys, m.nodesets,1); % 1 because we want ys0 -inputs=add(inputs,'adjoint',lambda_g,'doublevec',m.parameters.NumberOfDofsPerNode,m.parameters.NumberOfNodes); - -%Compute gradJ -grad_g=Gradj(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); -if (dim==3 & extrude_param), - displaystring(verbose,'%s',' extruding gradient...'); - grad_g=FieldExtrude(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,grad_g,'gradj',0); -end - -if m.parameters.ControlSteady, - %compute initial velocity from diagnostic_core (horiz+vertical) - displaystring(verbose,'%s',' compute 3d initial velocity...'); - results_diag=diagnostic_core(models); - u_g=results_diag.u_g; -end - -%Assign output -results.u_g=u_g; -results.grad_g=grad_g; Index: sm/trunk/src/m/solutions/jpl/modelsize.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/modelsize.m (revision 4130) +++ (revision ) @@ -1,20 +1,0 @@ -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: sm/trunk/src/m/solutions/jpl/objectivefunctionC.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/objectivefunctionC.m (revision 4130) +++ (revision ) @@ -1,25 +1,0 @@ -function J =objectivefunctionC(search_scalar,models,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.ControlType; -analysis_type=m.parameters.AnalysisType; - -%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.ControlType,parameter,'doublevec',1,m.parameters.NumberOfNodes); - -%Run diagnostic with updated parameters. -u_g=diagnostic_core_nonlinear(m,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=CostFunction(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); Index: sm/trunk/src/m/solutions/jpl/plot_direction.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/plot_direction.m (revision 4130) +++ (revision ) @@ -1,1 +1,0 @@ -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: sm/trunk/src/m/solutions/jpl/plot_newdistribution.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/plot_newdistribution.m (revision 4130) +++ (revision ) @@ -1,1 +1,0 @@ -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: sm/trunk/src/m/solutions/jpl/processresults.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/processresults.m (revision 4130) +++ (revision ) @@ -1,157 +1,0 @@ -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; - - elseif strcmpi(resultsname{j},'v_g'), - - v_g=results(i).v_g; - newresults(i).step=results(i).step; - newresults(i).time=results(i).time; - newresults(i).vel=v_g; - - elseif strcmpi(resultsname{j},'sx_g'), - - sy_g=results(i).sy_g; - newresults(i).step=results(i).step; - newresults(i).time=results(i).time; - newresults(i).slopey=sy_g; - - elseif strcmpi(resultsname{j},'sy_g'), - - sx_g=results(i).sx_g; - newresults(i).step=results(i).step; - newresults(i).time=results(i).time; - newresults(i).slopex=sx_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: sm/trunk/src/m/solutions/jpl/prognostic.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic.m (revision 4130) +++ (revision ) @@ -1,32 +1,0 @@ -function md=prognostic(md) -%PROGNOSITC - prognostic solution sequence. -% -% Usage: -% md=prognostic(md) - - %timing - t1=clock; - - %Build all models requested for diagnostic simulation - analysis_types=[PrognosticAnalysisEnum]; - solution_type=PrognosticAnalysisEnum; - - displaystring(md.verbose,'%s',['reading prognostic model data']); - femmodel=NewFemModel(md,solution_type,analysis_types,1); - - % figure out number of dof: just for information purposes. - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %compute solution - if ~qmu_analysis, - displaystring(verbose,'%s',['call computational core']); - femmodel=prognostic_core(femmodel); - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/prognostic2.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic2.m (revision 4130) +++ (revision ) @@ -1,31 +1,0 @@ -function md=prognostic2(md) -%PROGNOSITC2 - prognostic2 solution sequence. -% -% Usage: -% md=prognostic2(md) - - %timing - t1=clock; - - analysis_types=[Prognostic2AnalysisEnum]; - solution_type=Prognostic2AnalysisEnum; - - displaystring(md.verbose,'%s',['create finite element model']); - femmodel=NewFemModel(md,solution_type,analysis_types,1); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %compute solution - if ~qmu_analysis, - displaystring(verbose,'%s',['call computational core']); - femmodel=prognostic2_core(femmodel); - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/prognostic2_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic2_core.m (revision 4130) +++ (revision ) @@ -1,26 +1,0 @@ -function femmodel=prognostic2_core(femmodel) -%PROGNOSTIC2_CORE - linear solution sequence -% -% Usage: -% femmodel=prognostic2_core(femmodel) - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - - %Activate formulation - femmodel=SetCurrentAnalysis(femmodel,Prognostic2AnalysisEnum); - - displaystring(verbose,'\n%s',['depth averaging velocities...']); -% [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum); -% [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum); - - displaystring(verbose,'\n%s',['call computational core...']); - femmodel=solver_linear(femmodel); - - displaystring(verbose,'\n%s',['extude computed thickness on all layers...']); -% femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); - -end %end function Index: sm/trunk/src/m/solutions/jpl/prognostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/prognostic_core.m (revision 4130) +++ (revision ) @@ -1,26 +1,0 @@ -function femmodel=prognostic_core(femmodel) -%PROGNOSTIC_CORE - linear solution sequence -% -% Usage: -% femmodel=prognostic_core(femmodel) - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - - %Activate formulation - femmodel=SetCurrentAnalysis(femmodel,PrognosticAnalysisEnum); - - displaystring(verbose,'\n%s',['depth averaging velocities...']); - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum); - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum); - - displaystring(verbose,'\n%s',['call computational core...']); - femmodel=solver_linear(femmodel); - - displaystring(verbose,'\n%s',['extude computed thickness on all layers...']); - femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); - -end %end function Index: sm/trunk/src/m/solutions/jpl/steadystate.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/steadystate.m (revision 4130) +++ (revision ) @@ -1,42 +1,0 @@ -function md=steadystate(md); -%STEADYSTATE - compute the velocity and temperature field of a model in steady state. -% -% Usage: -% md=steadystate(md) -% - - %timing - t1=clock; - - analysis_types=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum]; - solution_type=SteadyStateAnalysisEnum; - - displaystring(md.verbose,'%s',['create finite element model']); - femmodel=NewFemModel(md,solution_type,analysis_types,7); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - control_analysis=femmodel.parameters.ControlAnalysis; - - %compute solution - if ~qmu_analysis, - if ~control_analysis, - - displaystring(verbose,'%s',['call computational core']); - femmodel=steadystate_core(femmodel); - - else, - - displaystring(verbose,'%s',['call computational core']); - femmodel=control_core(femmodel); - - end - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/steadystate_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/steadystate_core.m (revision 4130) +++ (revision ) @@ -1,53 +1,0 @@ -function femmodel=steadystate_core(femmodel); -%STEADYSTATE_CORE - compute the core temperature and velocity field at thermal steady state. -% -% Usage: -% femmodel=steadystate_core(femmodel); -% - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - dim=femmodel.parameters.Dim; - - %Initialize counter - step=1; - - while true, - - displaystring(verbose,'\n%s%i/\n','computing velocities and temperatures for step: ',step); - femmodel=thermal_core(femmodel); - - displaystring(verbose,'\n%s',['computing depth average temperature...']); - [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,TemperatureEnum,TemperatureAverageEnum); - - displaystring(verbose,'\n%s',['computing new velocity...']); - femmodel=diagnostic_core(femmodel); - - displaystring(verbose,'\n%s',['checking temperature, velocity and pressure convergence...']); - if step>1, - if steadystateconvergence(femmodel), - break; - end - end - - displaystring(verbose,'\n%s',['saving velocities, temperature and pressure for convergence...']); - femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,VxOldEnum); - femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,VyOldEnum); - femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,VzOldEnum); - femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureOlddnum); - femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureOldEnum); - - %Increase counter - step=step+1; - - end - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum); - end - -end %end of function Index: sm/trunk/src/m/solutions/jpl/surfaceslope.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/surfaceslope.m (revision 4130) +++ (revision ) @@ -1,35 +1,0 @@ -function md=surfaceslopecompute(md); -%SLOPECOMPUTE - compute the surface slope of a model -% -% Usage: -% md=surfaceslope(md) -% - %timing - t1=clock; - - analysis_types=SlopeAnalysisEnum; - solution_type=SlopeAnalysisEnum; - - displaystring(md.verbose,'%s',['create finite element model']); - femmodel=NewFemModel(md,solution_type,analysis_types,1); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %compute solution - if ~qmu_analysis, - displaystring(verbose,'%s',['call computational core']); - femmodel=surfaceslope_core(femmodel); - - displaystring(verbose,'%s',['write results']); - OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters); - - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/surfaceslope_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/surfaceslope_core.m (revision 4130) +++ (revision ) @@ -1,31 +1,0 @@ -function femmodel=surfaceslope_core(femmodel) -%SURFACESLOPE_CORE - core of the surface slope computation solution -% -% Usage: -% femmodel=surfaceslope_core(femmodel) -% - - %Recover some parameters: - verbose=femmodel.parameters.Verbose; - dim=femmodel.parameters.Dim; - - displaystring(verbose,'\n%s',['computing surface slope...']); - - %Call on core computations: - femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum); - femmodel=solver_linear(femmodel); - femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum); - femmodel=solver_linear(femmodel); - - %extrude inputs if we are in 3D: */ - if dim==3, - displaystring(verbose,'\n%s',['extruding surface slopein 3d...']); - femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum); - femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum); - } - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum); - InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum); - -end %end function Index: sm/trunk/src/m/solutions/jpl/thermal.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/thermal.m (revision 4130) +++ (revision ) @@ -1,31 +1,0 @@ -function md=thermal(md) -%THERMAL - thermal solution sequence: steady state and transient -% -% Usage: -% md=thermal(md) - - %timing - t1=clock; - - analysis_types=[ThermalAnalysisEnum,MeltingAnalysisEnum]; - solution_type=ThermalAnalysisEnum; - - displaystring(md.verbose,'%s',['reading thermal model data']); - femmodel=NewFemModel(md,solution_type,analysis_types,2); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %compute solution - if ~qmu_analysis, - displaystring(verbose,'%s',['call computational core']); - femmodel=thermal_core(femmodel); - else - %launch dakota driver for diagnostic core solution - Qmu(models,models.t.parameters); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/thermal_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/thermal_core.m (revision 4130) +++ (revision ) @@ -1,34 +1,0 @@ -function femmodel=thermal_core(femmodel) -%THERMAL_CORE - core of thermal solution -% -% Usage: -% femmodel=thermal_core(femmodel) - - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - ndt=femmodel.parameters.Ndt; - dt=femmodel.parameters.Dt; - - %Compute number of timesteps - if (dt==0 | ndt==0), - dt=0; - nsteps=1; - else - nsteps=floor(ndt/dt); - end - - %Loop through time - for i=1:nsteps+1, - displaystring(verbose,'\n%s%i/%i\n','time step: ',i,nsteps); - time=(i+1)*dt; - - displaystring(verbose,'\n%s',['computing temperature...']); - femmodel=thermal_core_step(femmodel,i,time); - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,i,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,i,time); - end - -end %end of function Index: sm/trunk/src/m/solutions/jpl/thermal_core_step.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/thermal_core_step.m (revision 4130) +++ (revision ) @@ -1,18 +1,0 @@ -function femmodel=thermal_core_step(femmodel) -%THERMAL_CORE_STEP - core of the thermal solution for one step -% -% Usage: -% femmodel=thermal_core_step(femmodel) - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - - displaystring(verbose,'\n%s',['computing temperature...']); - femmodel=SetCurrentAnalysis(femmodel,ThermalAnalysisEnum); - femmodel=solver_thermal_nonlineat(femmodel); - - displaystring(verbose,'\n%s',['compute melting...']); - femmodel=SetCurrentAnalysis(femmodel,MeltingAnalysisEnum); - femmodel=solver_linear(femmodel); - -end %end function Index: sm/trunk/src/m/solutions/jpl/transient2d.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/transient2d.m (revision 4130) +++ (revision ) @@ -1,34 +1,0 @@ -function md=transient2d(md); -%TRANSIENT2D - 2d transient solution -% -% Usage: -% md=transient(md) -% -% See also: TRANSIENT3D, TRANSIENT - - analysis_types=[DiagnosticHorizAnalysisEnum,PrognosticAnalysisEnum]; - solution_type=Transient2DAnalysisEnum; - - displaystring(md.verbose,'%s',['create finite element model']); - femmodel=NewFemModel(md,solution_type,analysis_types,2); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %compute solution - if ~qmu_analysis, - displaystring(verbose,'%s',['call computational core']); - femmodel=transient2d_core(femmodel); - - displaystring(verbose,'%s',['write results']); - OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters); - - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/transient2d_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/transient2d_core.m (revision 4130) +++ (revision ) @@ -1,41 +1,0 @@ -function femmodel=transient2d_core(femmodel) -%TRANSIENT2D_CORE - core of transient 2d solution -% -% Usage: -% femmodel=transient2d_core(femmodel) - - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - ndt=femmodel.parameters.Ndt; - dt=femmodel.parameters.Dt; - - %Initialize - time=0; - step=0; - - %Loop through time - while timeelements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=UpdateGeometry(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); - end - -end %end of function Index: sm/trunk/src/m/solutions/jpl/transient3d.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/transient3d.m (revision 4130) +++ (revision ) @@ -1,33 +1,0 @@ -function md=transient3d(md); -% -% Usage: -% md=transient3d(md) -% -% See also: TRANSIENT2D, TRANSIENT - - analysis_types=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,PrognosticAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum]; - solution_type=Transient3DAnalysisEnum; - - displaystring(md.verbose,'%s',['create finite element model']); - femmodel=NewFemModel(md,solution_type,analysis_types,8); - - %retrieve parameters - verbose=femmodel.parameters.Verbose; - qmu_analysis=femmodel.parameters.QmuAnalysis; - - %compute solution - if ~qmu_analysis, - displaystring(verbose,'%s',['call computational core']); - femmodel=transient3d_core(femmodel); - - displaystring(verbose,'%s',['write results']); - OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters); - - else - %launch dakota driver for diagnostic core solution - Qmu(femmodel); - end - - %stop timing - t2=clock; - displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']); Index: sm/trunk/src/m/solutions/jpl/transient3d_core.m =================================================================== --- /issm/trunk/src/m/solutions/jpl/transient3d_core.m (revision 4130) +++ (revision ) @@ -1,53 +1,0 @@ -function femmodel=transient3d_core(femmodel) -%TRANSIENT3D_CORE - core of transient 2d solution -% -% Usage: -% femmodel=transient3d_core(femmodel) - - - %recover parameters common to all solutions - verbose=femmodel.parameters.Verbose; - ndt=femmodel.parameters.Ndt; - dt=femmodel.parameters.Dt; - - %Initialize - time=0; - step=0; - - %Loop through time - while timeelements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=InputDepthAverage(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureAverageEnum); - - displaystring(verbose,'\n%s',['computing new velocities...']); - femmodel=diagnostic_core(femmodel); - - displaystring(verbose,'\n%s',['computing new thickness...']); - femmodel=prognostic_core(femmodel); - - displaystring(verbose,'\n%s',['updating geometry...']); - [femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=UpdateGeometry(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - - displaystring(verbose,'\n%s',['updating vertices position...']); - [femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=UpdateVertexPostions(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); - - displaystring(verbose,'\n%s',['saving results...']); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time); - InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingEnum,step,time); - end - -end %end of function Index: /issm/trunk/src/m/solutions/modelsize.m =================================================================== --- /issm/trunk/src/m/solutions/modelsize.m (revision 4131) +++ /issm/trunk/src/m/solutions/modelsize.m (revision 4131) @@ -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/objectivefunctionC.m =================================================================== --- /issm/trunk/src/m/solutions/objectivefunctionC.m (revision 4131) +++ /issm/trunk/src/m/solutions/objectivefunctionC.m (revision 4131) @@ -0,0 +1,25 @@ +function J =objectivefunctionC(search_scalar,models,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.ControlType; +analysis_type=m.parameters.AnalysisType; + +%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.ControlType,parameter,'doublevec',1,m.parameters.NumberOfNodes); + +%Run diagnostic with updated parameters. +u_g=diagnostic_core_nonlinear(m,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=CostFunction(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type); Index: /issm/trunk/src/m/solutions/plot_direction.m =================================================================== --- /issm/trunk/src/m/solutions/plot_direction.m (revision 4131) +++ /issm/trunk/src/m/solutions/plot_direction.m (revision 4131) @@ -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/plot_newdistribution.m =================================================================== --- /issm/trunk/src/m/solutions/plot_newdistribution.m (revision 4131) +++ /issm/trunk/src/m/solutions/plot_newdistribution.m (revision 4131) @@ -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/processresults.m =================================================================== --- /issm/trunk/src/m/solutions/processresults.m (revision 4131) +++ /issm/trunk/src/m/solutions/processresults.m (revision 4131) @@ -0,0 +1,157 @@ +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; + + elseif strcmpi(resultsname{j},'v_g'), + + v_g=results(i).v_g; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).vel=v_g; + + elseif strcmpi(resultsname{j},'sx_g'), + + sy_g=results(i).sy_g; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).slopey=sy_g; + + elseif strcmpi(resultsname{j},'sy_g'), + + sx_g=results(i).sx_g; + newresults(i).step=results(i).step; + newresults(i).time=results(i).time; + newresults(i).slopex=sx_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/prognostic.m =================================================================== --- /issm/trunk/src/m/solutions/prognostic.m (revision 4131) +++ /issm/trunk/src/m/solutions/prognostic.m (revision 4131) @@ -0,0 +1,32 @@ +function md=prognostic(md) +%PROGNOSITC - prognostic solution sequence. +% +% Usage: +% md=prognostic(md) + + %timing + t1=clock; + + %Build all models requested for diagnostic simulation + analysis_types=[PrognosticAnalysisEnum]; + solution_type=PrognosticAnalysisEnum; + + displaystring(md.verbose,'%s',['reading prognostic model data']); + femmodel=NewFemModel(md,solution_type,analysis_types,1); + + % figure out number of dof: just for information purposes. + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %compute solution + if ~qmu_analysis, + displaystring(verbose,'%s',['call computational core']); + femmodel=prognostic_core(femmodel); + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/prognostic2.m =================================================================== --- /issm/trunk/src/m/solutions/prognostic2.m (revision 4131) +++ /issm/trunk/src/m/solutions/prognostic2.m (revision 4131) @@ -0,0 +1,31 @@ +function md=prognostic2(md) +%PROGNOSITC2 - prognostic2 solution sequence. +% +% Usage: +% md=prognostic2(md) + + %timing + t1=clock; + + analysis_types=[Prognostic2AnalysisEnum]; + solution_type=Prognostic2AnalysisEnum; + + displaystring(md.verbose,'%s',['create finite element model']); + femmodel=NewFemModel(md,solution_type,analysis_types,1); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %compute solution + if ~qmu_analysis, + displaystring(verbose,'%s',['call computational core']); + femmodel=prognostic2_core(femmodel); + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/prognostic2_core.m =================================================================== --- /issm/trunk/src/m/solutions/prognostic2_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/prognostic2_core.m (revision 4131) @@ -0,0 +1,26 @@ +function femmodel=prognostic2_core(femmodel) +%PROGNOSTIC2_CORE - linear solution sequence +% +% Usage: +% femmodel=prognostic2_core(femmodel) + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + + %Activate formulation + femmodel=SetCurrentAnalysis(femmodel,Prognostic2AnalysisEnum); + + displaystring(verbose,'\n%s',['depth averaging velocities...']); +% [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum); +% [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum); + + displaystring(verbose,'\n%s',['call computational core...']); + femmodel=solver_linear(femmodel); + + displaystring(verbose,'\n%s',['extude computed thickness on all layers...']); +% femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + +end %end function Index: /issm/trunk/src/m/solutions/prognostic_core.m =================================================================== --- /issm/trunk/src/m/solutions/prognostic_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/prognostic_core.m (revision 4131) @@ -0,0 +1,26 @@ +function femmodel=prognostic_core(femmodel) +%PROGNOSTIC_CORE - linear solution sequence +% +% Usage: +% femmodel=prognostic_core(femmodel) + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + + %Activate formulation + femmodel=SetCurrentAnalysis(femmodel,PrognosticAnalysisEnum); + + displaystring(verbose,'\n%s',['depth averaging velocities...']); + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VxEnum,VxAverageEnum); + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,VyEnum,VyAverageEnum); + + displaystring(verbose,'\n%s',['call computational core...']); + femmodel=solver_linear(femmodel); + + displaystring(verbose,'\n%s',['extude computed thickness on all layers...']); + femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum); + +end %end function Index: /issm/trunk/src/m/solutions/steadystate.m =================================================================== --- /issm/trunk/src/m/solutions/steadystate.m (revision 4131) +++ /issm/trunk/src/m/solutions/steadystate.m (revision 4131) @@ -0,0 +1,42 @@ +function md=steadystate(md); +%STEADYSTATE - compute the velocity and temperature field of a model in steady state. +% +% Usage: +% md=steadystate(md) +% + + %timing + t1=clock; + + analysis_types=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum]; + solution_type=SteadyStateAnalysisEnum; + + displaystring(md.verbose,'%s',['create finite element model']); + femmodel=NewFemModel(md,solution_type,analysis_types,7); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + control_analysis=femmodel.parameters.ControlAnalysis; + + %compute solution + if ~qmu_analysis, + if ~control_analysis, + + displaystring(verbose,'%s',['call computational core']); + femmodel=steadystate_core(femmodel); + + else, + + displaystring(verbose,'%s',['call computational core']); + femmodel=control_core(femmodel); + + end + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/steadystate_core.m =================================================================== --- /issm/trunk/src/m/solutions/steadystate_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/steadystate_core.m (revision 4131) @@ -0,0 +1,53 @@ +function femmodel=steadystate_core(femmodel); +%STEADYSTATE_CORE - compute the core temperature and velocity field at thermal steady state. +% +% Usage: +% femmodel=steadystate_core(femmodel); +% + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + dim=femmodel.parameters.Dim; + + %Initialize counter + step=1; + + while true, + + displaystring(verbose,'\n%s%i/\n','computing velocities and temperatures for step: ',step); + femmodel=thermal_core(femmodel); + + displaystring(verbose,'\n%s',['computing depth average temperature...']); + [femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters]=InputDepthAverage(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,TemperatureEnum,TemperatureAverageEnum); + + displaystring(verbose,'\n%s',['computing new velocity...']); + femmodel=diagnostic_core(femmodel); + + displaystring(verbose,'\n%s',['checking temperature, velocity and pressure convergence...']); + if step>1, + if steadystateconvergence(femmodel), + break; + end + end + + displaystring(verbose,'\n%s',['saving velocities, temperature and pressure for convergence...']); + femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,VxOldEnum); + femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,VyOldEnum); + femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,VzOldEnum); + femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureOlddnum); + femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureOldEnum); + + %Increase counter + step=step+1; + + end + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum); + end + +end %end of function Index: /issm/trunk/src/m/solutions/surfaceslope.m =================================================================== --- /issm/trunk/src/m/solutions/surfaceslope.m (revision 4131) +++ /issm/trunk/src/m/solutions/surfaceslope.m (revision 4131) @@ -0,0 +1,35 @@ +function md=surfaceslopecompute(md); +%SLOPECOMPUTE - compute the surface slope of a model +% +% Usage: +% md=surfaceslope(md) +% + %timing + t1=clock; + + analysis_types=SlopeAnalysisEnum; + solution_type=SlopeAnalysisEnum; + + displaystring(md.verbose,'%s',['create finite element model']); + femmodel=NewFemModel(md,solution_type,analysis_types,1); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %compute solution + if ~qmu_analysis, + displaystring(verbose,'%s',['call computational core']); + femmodel=surfaceslope_core(femmodel); + + displaystring(verbose,'%s',['write results']); + OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters); + + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/surfaceslope_core.m =================================================================== --- /issm/trunk/src/m/solutions/surfaceslope_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/surfaceslope_core.m (revision 4131) @@ -0,0 +1,31 @@ +function femmodel=surfaceslope_core(femmodel) +%SURFACESLOPE_CORE - core of the surface slope computation solution +% +% Usage: +% femmodel=surfaceslope_core(femmodel) +% + + %Recover some parameters: + verbose=femmodel.parameters.Verbose; + dim=femmodel.parameters.Dim; + + displaystring(verbose,'\n%s',['computing surface slope...']); + + %Call on core computations: + femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeXAnalysisEnum); + femmodel=solver_linear(femmodel); + femmodel=SetCurrentAnalysisAlias(femmodel,SlopeAnalysisEnum,SurfaceSlopeYAnalysisEnum); + femmodel=solver_linear(femmodel); + + %extrude inputs if we are in 3D: */ + if dim==3, + displaystring(verbose,'\n%s',['extruding surface slopein 3d...']); + femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum); + femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum); + } + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum); + InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum); + +end %end function Index: /issm/trunk/src/m/solutions/thermal.m =================================================================== --- /issm/trunk/src/m/solutions/thermal.m (revision 4131) +++ /issm/trunk/src/m/solutions/thermal.m (revision 4131) @@ -0,0 +1,31 @@ +function md=thermal(md) +%THERMAL - thermal solution sequence: steady state and transient +% +% Usage: +% md=thermal(md) + + %timing + t1=clock; + + analysis_types=[ThermalAnalysisEnum,MeltingAnalysisEnum]; + solution_type=ThermalAnalysisEnum; + + displaystring(md.verbose,'%s',['reading thermal model data']); + femmodel=NewFemModel(md,solution_type,analysis_types,2); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %compute solution + if ~qmu_analysis, + displaystring(verbose,'%s',['call computational core']); + femmodel=thermal_core(femmodel); + else + %launch dakota driver for diagnostic core solution + Qmu(models,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/thermal_core.m =================================================================== --- /issm/trunk/src/m/solutions/thermal_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/thermal_core.m (revision 4131) @@ -0,0 +1,34 @@ +function femmodel=thermal_core(femmodel) +%THERMAL_CORE - core of thermal solution +% +% Usage: +% femmodel=thermal_core(femmodel) + + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + ndt=femmodel.parameters.Ndt; + dt=femmodel.parameters.Dt; + + %Compute number of timesteps + if (dt==0 | ndt==0), + dt=0; + nsteps=1; + else + nsteps=floor(ndt/dt); + end + + %Loop through time + for i=1:nsteps+1, + displaystring(verbose,'\n%s%i/%i\n','time step: ',i,nsteps); + time=(i+1)*dt; + + displaystring(verbose,'\n%s',['computing temperature...']); + femmodel=thermal_core_step(femmodel,i,time); + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,i,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,i,time); + end + +end %end of function Index: /issm/trunk/src/m/solutions/thermal_core_step.m =================================================================== --- /issm/trunk/src/m/solutions/thermal_core_step.m (revision 4131) +++ /issm/trunk/src/m/solutions/thermal_core_step.m (revision 4131) @@ -0,0 +1,18 @@ +function femmodel=thermal_core_step(femmodel) +%THERMAL_CORE_STEP - core of the thermal solution for one step +% +% Usage: +% femmodel=thermal_core_step(femmodel) + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + + displaystring(verbose,'\n%s',['computing temperature...']); + femmodel=SetCurrentAnalysis(femmodel,ThermalAnalysisEnum); + femmodel=solver_thermal_nonlineat(femmodel); + + displaystring(verbose,'\n%s',['compute melting...']); + femmodel=SetCurrentAnalysis(femmodel,MeltingAnalysisEnum); + femmodel=solver_linear(femmodel); + +end %end function Index: /issm/trunk/src/m/solutions/transient2d.m =================================================================== --- /issm/trunk/src/m/solutions/transient2d.m (revision 4131) +++ /issm/trunk/src/m/solutions/transient2d.m (revision 4131) @@ -0,0 +1,34 @@ +function md=transient2d(md); +%TRANSIENT2D - 2d transient solution +% +% Usage: +% md=transient(md) +% +% See also: TRANSIENT3D, TRANSIENT + + analysis_types=[DiagnosticHorizAnalysisEnum,PrognosticAnalysisEnum]; + solution_type=Transient2DAnalysisEnum; + + displaystring(md.verbose,'%s',['create finite element model']); + femmodel=NewFemModel(md,solution_type,analysis_types,2); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %compute solution + if ~qmu_analysis, + displaystring(verbose,'%s',['call computational core']); + femmodel=transient2d_core(femmodel); + + displaystring(verbose,'%s',['write results']); + OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters); + + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/transient2d_core.m =================================================================== --- /issm/trunk/src/m/solutions/transient2d_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/transient2d_core.m (revision 4131) @@ -0,0 +1,41 @@ +function femmodel=transient2d_core(femmodel) +%TRANSIENT2D_CORE - core of transient 2d solution +% +% Usage: +% femmodel=transient2d_core(femmodel) + + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + ndt=femmodel.parameters.Ndt; + dt=femmodel.parameters.Dt; + + %Initialize + time=0; + step=0; + + %Loop through time + while timeelements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=UpdateGeometry(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); + end + +end %end of function Index: /issm/trunk/src/m/solutions/transient3d.m =================================================================== --- /issm/trunk/src/m/solutions/transient3d.m (revision 4131) +++ /issm/trunk/src/m/solutions/transient3d.m (revision 4131) @@ -0,0 +1,33 @@ +function md=transient3d(md); +% +% Usage: +% md=transient3d(md) +% +% See also: TRANSIENT2D, TRANSIENT + + analysis_types=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,PrognosticAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum]; + solution_type=Transient3DAnalysisEnum; + + displaystring(md.verbose,'%s',['create finite element model']); + femmodel=NewFemModel(md,solution_type,analysis_types,8); + + %retrieve parameters + verbose=femmodel.parameters.Verbose; + qmu_analysis=femmodel.parameters.QmuAnalysis; + + %compute solution + if ~qmu_analysis, + displaystring(verbose,'%s',['call computational core']); + femmodel=transient3d_core(femmodel); + + displaystring(verbose,'%s',['write results']); + OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters); + + else + %launch dakota driver for diagnostic core solution + Qmu(femmodel); + 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/transient3d_core.m =================================================================== --- /issm/trunk/src/m/solutions/transient3d_core.m (revision 4131) +++ /issm/trunk/src/m/solutions/transient3d_core.m (revision 4131) @@ -0,0 +1,53 @@ +function femmodel=transient3d_core(femmodel) +%TRANSIENT3D_CORE - core of transient 2d solution +% +% Usage: +% femmodel=transient3d_core(femmodel) + + + %recover parameters common to all solutions + verbose=femmodel.parameters.Verbose; + ndt=femmodel.parameters.Ndt; + dt=femmodel.parameters.Dt; + + %Initialize + time=0; + step=0; + + %Loop through time + while timeelements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=InputDepthAverage(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureAverageEnum); + + displaystring(verbose,'\n%s',['computing new velocities...']); + femmodel=diagnostic_core(femmodel); + + displaystring(verbose,'\n%s',['computing new thickness...']); + femmodel=prognostic_core(femmodel); + + displaystring(verbose,'\n%s',['updating geometry...']); + [femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=UpdateGeometry(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + displaystring(verbose,'\n%s',['updating vertices position...']); + [femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters]=UpdateVertexPostions(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters); + + displaystring(verbose,'\n%s',['saving results...']); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,SurfaceEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,BedEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,step,time); + InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingEnum,step,time); + end + +end %end of function