Changeset 4126


Ignore:
Timestamp:
06/22/10 16:02:09 (15 years ago)
Author:
seroussi
Message:

core solutions in matlab

Location:
issm/trunk/src/m/solutions/jpl
Files:
3 added
1 deleted
12 edited

Legend:

Unmodified
Added
Removed

  • issm/trunk/src/m/solutions/jpl/balancedthickness2_core.m

    r3839 r4126  
    1 function results=balancedthickness2_core(models,analysis_type,sub_analysis_type)
     1function femmodel=balancedthickness2_core(femmodel)
    22%BALANCEDTHICKNESS_CORE - linear solution sequence
    33%
    44%   Usage:
    5 %      h_g=balancedthickness2_core(m,analysis_type,sub_analysis_type)
     5%      femmodel=balancedthickness2_core(femmodel)
    66
    7         %get FE model
    8         m=models.p;
    9         results.time=0;
    10         results.step=1;
     7        %recover parameters common to all solutions
     8        verbose=femmodel.parameters.Verbose;
     9        dim=femmodel.parameters.Dim;
    1110
    12         displaystring(m.parameters.Verbose,'\n%s',['depth averaging velocity...']);
    13         %Take only the first two dofs of m.parameters.u_g
    14         vx_g=get(inputs,'vx',1);
    15         vy_g=get(inputs,'vy',1);
     11        %Activate formulation
     12        femmodel=SetCurrentAnalysis(femmodel,Balancedthickness2AnalysisEnum);
    1613
    17         %NOT WORKING YET!!!!!
    18         %vx_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vx_g,'vx');
    19         %vy_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vy_g,'vy');
     14        displaystring(verbose,'\n%s',['depth averaging velocities...']);
     15        %[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);
     16        %[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);
     17        if dim==3,
     18        %       [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);
     19        end
    2020
    21         inputs=add(inputs,'vx_average',vx_g,'doublevec',1,m.parameters.NumberOfVertices);
    22         inputs=add(inputs,'vy_average',vy_g,'doublevec',1,m.parameters.NumberOfVertices);
     21        displaystring(verbose,'\n%s',['call computational core...']);
     22        femmodel=solver_linear(femmodel);
     23       
     24        displaystring(verbose,'\n%s',['extude computed thickness on all layers...']);
     25        %femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
    2326
    24         displaystring(m.parameters.Verbose,'\n%s',['call computational core:']);
    25         results.h_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type);
    26 
    27         displaystring(m.parameters.Verbose,'\n%s',['averaging over vertices']);
    28         results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.h_g);
    29 
    30         displaystring(m.parameters.Verbose,'\n%s',['extrude computed thickness on all layers:']);
    31         %results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness');
    32 
     27        displaystring(verbose,'\n%s',['saving results...']);
     28        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
     29       
    3330end %end function

  • issm/trunk/src/m/solutions/jpl/balancedthickness_core.m

    r3975 r4126  
    1 function results=balancedthickness_core(models,analysis_type,sub_analysis_type)
     1function femmodel=balancedthickness_core(femmodel)
    22%BALANCEDTHICKNESS_CORE - linear solution sequence
    33%
    44%   Usage:
    5 %      h_g=balancedthickness_core(m,analysis_type,sub_analysis_type)
     5%      femmodel=balancedthickness_core(femmode)
    66
    7         %get FE model
    8         verbose=models.bt.parameters.Verbose;
    9         results.time=0;
    10         results.step=1;
     7        %recover parameters common to all solutions
     8        verbose=femmodel.parameters.Verbose;
     9        dim=femmodel.parameters.Dim;
    1110
    12         displaystring(verbose,'\n%s',['depth averaging velocity...']);
    13         [models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters]=DepthAverageInput(models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters,VxEnum);
    14         [models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters]=DepthAverageInput(models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters,VyEnum);
     11        %Activate formulation
     12        femmodel=SetCurrentAnalysis(femmodel,BalancedthicknessAnalysisEnum);
    1513
    16         displaystring(verbose,'\n%s',['call computational core:']);
    17         h_g=diagnostic_core_linear(models.bt,analysis_type,sub_analysis_type);
     14        displaystring(verbose,'\n%s',['depth averaging velocities...']);
     15        [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);
     16        [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);
     17        if dim==3,
     18                [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);
     19        end
    1820
    19         %Update
    20         models.bt.elements=UpdateInputsFromSolution(models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters,h_g,BalancedthicknessAnalysisEnum,NoneAnalysisEnum);
     21        displaystring(verbose,'\n%s',['call computational core...']);
     22        femmodel=solver_linear(femmodel);
     23       
     24        displaystring(verbose,'\n%s',['extude computed thickness on all layers...']);
     25        femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
    2126
    22         displaystring(verbose,'\n%s',['extrude computed thickness on all layers:']);
    23         models.bt.elements=ExtrudeInput(models.bt.elements,models.bt.nodes,models.bt.vertices,models.bt.loads,models.bt.materials,models.bt.parameters,ThicknessEnum);
    24 
    25         %NEED TO BE CLEANED
    26         results.h_g=h_g;
     27        displaystring(verbose,'\n%s',['saving results...']);
     28        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
     29       
    2730end %end function

  • issm/trunk/src/m/solutions/jpl/balancedvelocities_core.m

    r3839 r4126  
    1 function results=balancedvelocities_core(models,analysis_type,sub_analysis_type)
     1function femmodel=balancedvelocities_core(femmdoel)
    22%BALANCEDVELOCITIES_CORE - linear solution sequence
    33%
    44%   Usage:
    5 %      v_g=balancedvelocities_core(m,analysis_type,sub_analysis_type)
     5%      femmodel=balancedvelocities_core(femmodel)
    66
    7         %get FE model
    8         m=models.p;
    9         results.time=0;
    10         results.step=1;
     7        %recover parameters common to all solutions
     8        verbose=femmodel.parameters.Verbose;
     9        dim=femmodel.parameters.Dim;
    1110
    12         displaystring(m.parameters.Verbose,'\n%s',['depth averaging velocity...']);
    13         %Take only the first two dofs of m.parameters.u_g
    14         u_g=get(inputs,'velocity',[1 1 0 0]);
    15         u_g=FieldDepthAverage(m.elements,m.nodes,m.loads,m.materials,m.parameters,u_g,'velocity');
    16         inputs=add(inputs,'velocity_average',u_g,'doublevec',2,m.parameters.NumberOfNodes);
     11        %Activate formulation
     12        femmodel=SetCurrentAnalysis(femmodel,BalancedvelocitiesAnalysisEnum);
    1713
    18         displaystring(m.parameters.Verbose,'\n%s',['call computational core:']);
    19         results.h_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type);
     14        displaystring(verbose,'\n%s',['depth averaging velocities...']);
     15        [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);
     16        [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);
     17        if dim==3,
     18                [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);
     19        end
    2020
    21         displaystring(m.parameters.Verbose,'\n%s',['extrude computed thickness on all layers:']);
    22         results.v_g=FieldExtrude(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness',0);
     21        displaystring(verbose,'\n%s',['call computational core...']);
     22        femmodel=solver_linear(femmodel);
     23       
     24        displaystring(verbose,'\n%s',['extude computed thickness on all layers...']);
     25        femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
     26        femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
     27
     28        displaystring(verbose,'\n%s',['saving results...']);
     29        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
     30        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
    2331
    2432end %end function

  • issm/trunk/src/m/solutions/jpl/bedslope_core.m

    r4113 r4126  
    2222        if dim==3,
    2323                displaystring(verbose,'\n%s',['extruding bed slopein 3d...']);
    24                 InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
    25                 InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);
     24                femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
     25                femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);
    2626        }
    2727       

  • issm/trunk/src/m/solutions/jpl/diagnostic_core.m

    r4113 r4126  
    2121        %for qmu analysis, be sure the velocity input we are starting from  is the one in the parameters:
    2222        if qmu_analysis,
    23                 InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVxEnum,VxEnum);
    24                 InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVyEnum,VyEnum);
    25                 InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVzEnum,VzEnum);
    26                 InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuPressureEnum,PressureEnum);
     23                femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVxEnum,VxEnum);
     24                femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVyEnum,VyEnum);
     25                femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuVzEnum,VzEnum);
     26                femmodel.elements=InputDuplicate(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,QmuPressureEnum,PressureEnum);
    2727        end
    2828
    … …  
    5858
    5959                        %"recondition" pressure computed previously:
    60                         InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum);
     60                        femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureStokesEnum);
    6161                        InputScale(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureStokesEnum,1.0/stokesreconditioning);
    6262

  • issm/trunk/src/m/solutions/jpl/prognostic2_core.m

    r3839 r4126  
    1 function results=prognostic2_core(models,analysis_type,sub_analysis_type)
     1function femmodel=prognostic2_core(femmodel)
    22%PROGNOSTIC2_CORE - linear solution sequence
    33%
    44%   Usage:
    5 %      h_g=prognostic2_core(m,analysis_type,sub_analysis_type)
     5%      femmodel=prognostic2_core(femmodel)
    66
    7         %get FE model
    8         m=models.p;
    9         results.time=0;
    10         results.step=1;
     7        %recover parameters common to all solutions
     8        verbose=femmodel.parameters.Verbose;
    119
    12         displaystring(m.parameters.Verbose,'\n%s',['depth averaging velocity...']);
    13         %Take only the first two dofs of m.parameters.u_g
    14         vx_g=get(inputs,'vx',1);
    15         vy_g=get(inputs,'vy',1);
     10        %Activate formulation
     11        femmodel=SetCurrentAnalysis(femmodel,Prognostic2AnalysisEnum);
    1612
    17         %NOT WORKING YET!!!!!
    18         %vx_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vx_g,'vx');
    19         %vy_g=FieldDepthAverage(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,vy_g,'vy');
     13        displaystring(verbose,'\n%s',['depth averaging velocities...']);
     14%       [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);
     15%       [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);
    2016
    21         inputs=add(inputs,'vx_average',vx_g,'doublevec',1,m.parameters.NumberOfVertices);
    22         inputs=add(inputs,'vy_average',vy_g,'doublevec',1,m.parameters.NumberOfVertices);
     17        displaystring(verbose,'\n%s',['call computational core...']);
     18        femmodel=solver_linear(femmodel);
     19       
     20        displaystring(verbose,'\n%s',['extude computed thickness on all layers...']);
     21%       femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
    2322
    24         displaystring(m.parameters.Verbose,'\n%s',['call computational core:']);
    25         results.h_g=diagnostic_core_linear(m,analysis_type,sub_analysis_type);
    26 
    27         displaystring(m.parameters.Verbose,'\n%s',['averaging over vertices']);
    28         results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,results.h_g);
    29 
    30         displaystring(m.parameters.Verbose,'\n%s',['extrude computed thickness on all layers:']);
    31         %results.h_g=FieldAverageOntoVertices(m.elements,m.nodes,m.loads,m.materials,m.parameters,results.h_g,'thickness');
    32 
     23        displaystring(verbose,'\n%s',['saving results...']);
     24        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
     25       
    3326end %end function

  • issm/trunk/src/m/solutions/jpl/prognostic_core.m

    r3885 r4126  
    1 function results=prognostic_core(models,analysis_type,sub_analysis_type)
     1function femmodel=prognostic_core(femmodel)
    22%PROGNOSTIC_CORE - linear solution sequence
    33%
    44%   Usage:
    5 %      results=prognostic_core(m,analysis_type,sub_analysis_type)
     5%      femmodel=prognostic_core(femmodel)
    66
    7         %get FE model
    8         verbose=models.p.parameters.Verbose;
    9         results.time=0;
    10         results.step=1;
     7        %recover parameters common to all solutions
     8        verbose=femmodel.parameters.Verbose;
    119
    12         displaystring(verbose,'\n%s',['depth averaging velocity...']);
    13         [models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters]=DepthAverageInput(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,VxEnum);
    14         [models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters]=DepthAverageInput(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,VyEnum);
     10        %Activate formulation
     11        femmodel=SetCurrentAnalysis(femmodel,PrognosticAnalysisEnum);
    1512
    16         displaystring(verbose,'\n%s',['call computational core:']);
    17         h_g=diagnostic_core_linear(models.p,analysis_type,sub_analysis_type);
     13        displaystring(verbose,'\n%s',['depth averaging velocities...']);
     14        [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);
     15        [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);
    1816
    19         %Update
    20         models.p.elements=UpdateInputsFromSolution(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,h_g,PrognosticAnalysisEnum,NoneAnalysisEnum);
     17        displaystring(verbose,'\n%s',['call computational core...']);
     18        femmodel=solver_linear(femmodel);
     19       
     20        displaystring(verbose,'\n%s',['extude computed thickness on all layers...']);
     21        femmodel.elements=InputExtrude(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
    2122
    22         displaystring(verbose,'\n%s',['extrude computed thickness on all layers:']);
    23         models.p.elements=ExtrudeInput(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,ThicknessEnum);
    24 
    25         %NEED TO BE CLEANED
    26         results.h_g=h_g;
     23        displaystring(verbose,'\n%s',['saving results...']);
     24        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,ThicknessEnum);
    2725
    2826end %end function

  • issm/trunk/src/m/solutions/jpl/steadystate_core.m

    r3839 r4126  
    1 function results=steadystate_core(models);
     1function femmodel=steadystate_core(femmodel);
    22%STEADYSTATE_CORE - compute the core temperature and velocity field  at thermal steady state.
    33%
    44%   Usage:
    5 %      results=steadystate_core(models);
     5%      femmodel=steadystate_core(femmodel);
    66%
    77
    8 %recover models
    9 m_dh=models.dh;
    10 m_dv=models.dv;
    11 m_ds=models.ds;
    12 m_dhu=models.dhu;
    13 m_sl=models.sl;
    14 m_t=models.t;
    15 m_m=models.m;
     8        %recover parameters common to all solutions
     9        verbose=femmodel.parameters.Verbose;
     10        dim=femmodel.parameters.Dim;
    1611
    17 %recover parameters common to all solutions
    18 verbose=m_dhu.parameters.Verbose;
    19 dim=m_dhu.parameters.Dim;
    20 eps_rel=m_dhu.parameters.EpsRel;
    21 ishutter=m_dhu.parameters.IsHutter;
    22 ismacayealpattyn=m_dh.parameters.IsMacAyealPattyn;
    23 isstokes=m_ds.parameters.IsStokes;
     12        %Initialize counter
     13        step=1;
    2414
    25 %convergence
    26 converged=0;
     15        while true,
    2716
    28 step=1;
     17                displaystring(verbose,'\n%s%i/\n','computing velocities and temperatures for step: ',step);
     18                femmodel=thermal_core(femmodel);
    2919
    30 %initialize:
    31 t_g_old=0;
    32 u_g_old=0;
     20                displaystring(verbose,'\n%s',['computing depth average temperature...']);
     21                [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);
    3322
    34 if isstokes, ndof=4; else ndof=3; end
     23                displaystring(verbose,'\n%s',['computing new velocity...']);
     24                femmodel=diagnostic_core(femmodel);
    3525
    36 while(~converged),
    37        
    38         displaystring(verbose,'%s%i','computing temperature and velocity for step ',step);
    39        
    40         %first compute temperature at steady state.
    41         if (step>1),
    42                 inputs=add(inputs,'velocity',results_diagnostic.u_g,'doublevec',ndof,m_t.parameters.NumberOfNodes);
    43         end
    44         results_thermal=thermal_core(models);
    45        
    46         %add temperature to inputs.
    47         %Compute depth averaged temperature
    48         temperature_average=FieldDepthAverage(m_t.elements,m_t.nodes,m_t.vertices,m_t.loads,m_t.materials,m_t.parameters,results_thermal.t_g,'temperature');
    49         inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,m_t.parameters.NumberOfNodes);
    50         inputs=add(inputs,'temperature',results_thermal.t_g,'doublevec',1,m_t.parameters.NumberOfNodes);
    51        
    52         %now compute diagnostic velocity using the steady state temperature.
    53         results_diagnostic=diagnostic_core(models);
    54        
    55         %convergence?
    56         du_g=results_diagnostic.u_g-u_g_old;
    57         ndu=norm(du_g,2);
    58         nu=norm(u_g_old,2);
    59        
    60         dt_g=results_thermal.t_g-t_g_old;
    61         ndt=norm(dt_g,2);
    62         nt=norm(t_g_old,2);
     26                displaystring(verbose,'\n%s',['checking temperature, velocity and pressure convergence...']);
     27                if step>1,
     28                        if steadystateconvergence(femmodel),
     29                                break;
     30                        end
     31                end
    6332
    64         u_g_old=results_diagnostic.u_g;
    65         t_g_old=results_thermal.t_g;
    66        
    67         displaystring(verbose,'%-60s%g\n                                     %s%g\n                                     %s%g%s',...
    68                       '      relative convergence criterion: velocity -> norm(du)/norm(u)=   ',ndu/nu*100,' temperature -> norm(dt)/norm(t)=',ndt/nt*100,' eps_rel:                        ',eps_rel*100,' %');
    69         if (ndu/nu<=eps_rel)  & (ndt/nt<=eps_rel),
    70                 converged=1;
    71         else
    72                 converged=0;
     33                displaystring(verbose,'\n%s',['saving velocities, temperature and pressure for convergence...']);
     34                femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum,VxOldEnum);
     35                femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum,VyOldEnum);
     36                femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum,VzOldEnum);
     37                femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum,PressureOlddnum);
     38                femmodel.elements=InputDuplicate(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,TemperatureOldEnum);
     39
     40                %Increase counter
     41                step=step+1;
     42
    7343        end
    7444
    75         step=step+1;
    76         if converged,
    77                 break;
     45        displaystring(verbose,'\n%s',['saving results...']);
     46        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VxEnum);
     47        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VyEnum);
     48        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,VzEnum);
     49        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,PressureEnum);
     50        InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum);
    7851        end
    79 end
    8052
    81 %save results from thermal and diagnostic
    82 results.u_g=results_diagnostic.u_g;
    83 results.p_g=results_diagnostic.p_g;
    84 results.t_g=results_thermal.t_g;
    85 results.m_g=results_thermal.m_g;
    86 results.step=step;
    87 results.time=0;
     53end %end of function

  • issm/trunk/src/m/solutions/jpl/surfaceslope_core.m

    r4113 r4126  
    55%      femmodel=surfaceslope_core(femmodel)
    66%
    7 
    87
    98        %Recover some parameters:
    … …  
    2221        if dim==3,
    2322                displaystring(verbose,'\n%s',['extruding surface slopein 3d...']);
    24                 InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
    25                 InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);
     23                femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
     24                femmodel.elements=InputExtrude(femmodel.elements,femmodel.nodes, femmodel.vertices,femmodel.loads, femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);
    2625        }
    2726       
    … …  
    2928        InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeXEnum);
    3029        InputToResult(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,SurfaceSlopeYEnum);
     30
     31end %end function

  • issm/trunk/src/m/solutions/jpl/thermal_core.m

    r3915 r4126  
    1 function results=thermal_core(models)
     1function femmodel=thermal_core(femmodel)
    22%THERMAL_CORE - core of thermal solution
    33%
    44%   Usage:
    5 %      solution=thermal_core(models)
     5%      femmodel=thermal_core(femmodel)
    66
    7 %recover parameters common to all solutions
    8 verbose=models.t.parameters.Verbose;
    97
    10 if models.t.parameters.Dt==0,
     8        %recover parameters common to all solutions
     9        verbose=femmodel.parameters.Verbose;
     10        ndt=femmodel.parameters.Ndt;
     11        dt=femmodel.parameters.Dt;
    1112
    12         results.time=0;
    13         results.step=1;
    14 
    15         displaystring(verbose,'\n%s',['computing temperatures...']);
    16         [t_g models.t.loads melting_offset]=thermal_core_nonlinear(models.t,ThermalAnalysisEnum(),NoneAnalysisEnum());
    17 
    18         displaystring(verbose,'\n%s',['computing melting...']);
    19         models=ModelUpdateInputsFromVector(models,t_g,TemperatureEnum,VertexEnum);
    20         [models.m.elements models.m.loads]=UpdateInputsFromConstant(models.m.elements,models.m.nodes,models.m.vertices,models.m.loads,models.m.materials,models.m.parameters,melting_offset,MeltingOffsetEnum);
    21         m_g=diagnostic_core_linear(models.m,MeltingAnalysisEnum(),NoneAnalysisEnum());
    22 
    23         %NEED TO BE CLEANED
    24         results.t_g=t_g;
    25         results.m_g=m_g;
    26 
    27 else
    28 
    29         %initialize temperature and melting
    30         nsteps=models.t.parameters.Ndt/models.t.parameters.Dt;
    31         results.step=1;
    32         results.time=0;
    33         results.t_g=[]; %FIRST VALUE MISSING
    34         results.m_g=[];
    35 
    36         for n=1:nsteps,
    37 
    38                 displaystring(verbose,'\n%s%i/%i\n','time step: ',n,nsteps);
    39                 results(n+1).step=n+1;
    40                 results(n+1).time=n*models.t.parameters.Dt;
    41 
    42                 displaystring(verbose,'\n%s',['    computing temperatures...']);
    43                 [t_g models.t.loads melting_offset]=thermal_core_nonlinear(models.t,ThermalAnalysisEnum(),NoneAnalysisEnum());
    44 
    45                 displaystring(verbose,'\n%s',['    computing melting...']);
    46                 models=ModelUpdateInputsFromVector(models,t_g,TemperatureEnum,VertexEnum);
    47                 [models.m.elements models.m.loads]=UpdateInputsFromConstant(models.m.elements,models.m.nodes,models.m.vertices,models.m.loads,models.m.materials,models.m.parameters,melting_offset,MeltingOffsetEnum);
    48                 m_g=diagnostic_core_linear(models.m,MeltingAnalysisEnum(),NoneAnalysisEnum());
    49 
    50                 %NEED TO BE CLEANED
    51                 results(n+1).t_g=t_g;
    52                 results(n+1).m_g=m_g;
    53 
     13        %Compute number of timesteps
     14        if (dt==0 | ndt==0),
     15                dt=0;
     16                nsteps=1;
     17        else
     18                nsteps=floor(ndt/dt);
    5419        end
    5520
    56 end
     21        %Loop through time
     22        for i=1:nsteps+1,
     23                displaystring(verbose,'\n%s%i/%i\n','time step: ',i,nsteps);
     24                time=(i+1)*dt;
     25
     26                displaystring(verbose,'\n%s',['computing temperature...']);
     27                femmodel=thermal_core_step(femmodel,i,time);
     28
     29                displaystring(verbose,'\n%s',['saving results...']);
     30                InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,TemperatureEnum,i,time);
     31                InputToResult(femmodel->elements,femmodel->nodes,femmodel->vertices,femmodel->loads,femmodel->materials,femmodel->parameters,MeltingRateEnum,i,time);
     32        end
     33
     34end %end of function

  • issm/trunk/src/m/solutions/jpl/transient2d.m

    r4102 r4126  
    77%   See also: TRANSIENT3D, TRANSIENT
    88
    9 %Build all models requested
    10 models.analysis_type='transient'; %needed for processresults
     9        analysis_types=[DiagnosticHorizAnalysisEnum,PrognosticAnalysisEnum];
     10        solution_type=Transient2DAnalysisEnum;
    1111
    12 displaystring(md.verbose,'%s',['reading diagnostic horiz model data']);
    13 md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=HorizAnalysisEnum(); models.dh=CreateFemModel(md);
     12        displaystring(md.verbose,'%s',['create finite element model']);
     13        femmodel=NewFemModel(md,solution_type,analysis_types,2);
    1414
    15 displaystring(md.verbose,'\n%s',['reading diagnostic vert model data']);
    16 md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=VertAnalysisEnum(); models.dv=CreateFemModel(md);
     15        %retrieve parameters
     16        verbose=femmodel.parameters.Verbose;
     17        qmu_analysis=femmodel.parameters.QmuAnalysis;
    1718
    18 displaystring(md.verbose,'\n%s',['reading diagnostic stokes model data']);
    19 md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=StokesAnalysisEnum(); models.ds=CreateFemModel(md);
     19        %compute solution
     20        if ~qmu_analysis,
     21                displaystring(verbose,'%s',['call computational core']);
     22                femmodel=transient2d_core(femmodel);
     23               
     24                displaystring(verbose,'%s',['write results']);
     25                OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters);
    2026
    21 displaystring(md.verbose,'\n%s',['reading diagnostic hutter model data']);
    22 md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=HutterAnalysisEnum(); models.dhu=CreateFemModel(md);
     27        else
     28                %launch dakota driver for diagnostic core solution
     29                Qmu(femmodel);
     30        end
    2331
    24 displaystring(md.verbose,'\n%s',['reading surface and bed slope computation model data']);
    25 md.analysis_type=SlopecomputeAnalysisEnum(); md.sub_analysis_type=NoneAnalysisEnum(); models.sl=CreateFemModel(md);
    26 
    27 displaystring(md.verbose,'%s',['reading prognostic model data']);
    28 md.analysis_type=PrognosticAnalysisEnum(); models.p=CreateFemModel(md);
    29 
    30 %initialize solution
    31 solution=struct('step',[],'time',[],'u_g',[],'p_g',[],'h_g',[],'s_g',[],'b_g',[]);
    32 solution.step=1;
    33 solution.time=0;
    34 %solution.u_g=models.dh.parameters.u_g(dofsetgen([1,2],3,length(models.dh.parameters.u_g)));
    35 solution.u_g=[];
    36 solution.p_g=[];
    37 %solution.h_g=models.p.parameters.h_g;
    38 solution.h_g=md.thickness;
    39 solution.s_g=md.surface;
    40 %solution.s_g=models.p.parameters.s_g;
    41 solution.b_g=md.bed;
    42 %solution.b_g=models.p.parameters.b_g;
    43 
    44 %initialize inputs
    45 displaystring(md.verbose,'\n%s',['setup inputs...']);
    46 inputs=inputlist;
    47 inputs=add(inputs,'velocity',solution.u_g,'doublevec',2,models.p.parameters.NumberOfNodes);
    48 inputs=add(inputs,'melting',models.p.parameters.m_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    49 inputs=add(inputs,'accumulation',models.p.parameters.a_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    50 inputs=add(inputs,'dt',models.p.parameters.Dt*models.p.parameters.Yts,'double'); %change time from yrs to sec
    51 
    52 % figure out number of dof: just for information purposes.
    53 md.dof=modelsize(models);
    54 
    55 %first time step is given by model.
    56 dt=models.p.parameters.Dt;
    57 finaltime=models.p.parameters.Ndt;
    58 time=dt;
    59 yts=models.p.parameters.Yts;
    60 n=1; %counter
    61 
    62 while  time<finaltime+dt, %make sure we run up to finaltime.
    63 
    64         disp(sprintf('\n%s%g%s%g%s%g\n','time [yr]: ',time,'    iteration number: ',n,'/',floor(finaltime/dt)));
    65 
    66         solution(n+1).step=n+1;
    67         solution(n+1).time=time;
    68 
    69         %update inputs
    70         inputs=add(inputs,'thickness',solution(n).h_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    71         inputs=add(inputs,'surface',solution(n).s_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    72         inputs=add(inputs,'bed',solution(n).b_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    73 
    74         %Deal with velocities.
    75 
    76         %Get horizontal solution.
    77         rawresults=diagnostic_core(models);
    78         solution(n+1).u_g=rawresults.u_g; solution(n+1).p_g=rawresults.p_g;
    79 
    80         %compute new thickness
    81         disp(sprintf('%s','   computing new thickness...'));
    82         inputs=add(inputs,'velocity',solution(n+1).u_g,'doublevec',2,models.p.parameters.NumberOfNodes);
    83         rawresults=prognostic_core(models,PrognosticAnalysisEnum(),NoneAnalysisEnum());
    84         new_thickness=rawresults.h_g;
    85 
    86         %update surface and bed using the new thickness
    87         disp(sprintf('%s','   updating geometry...'));
    88         [new_thickness,new_bed,new_surface]=UpdateGeometry(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,new_thickness,solution(n).b_g,solution(n).s_g);
    89 
    90         %Record bed surface and thickness in the solution
    91         solution(n+1).h_g=new_thickness;
    92         solution(n+1).b_g=new_bed;
    93         solution(n+1).s_g=new_surface;
    94 
    95         %figure out if time stepping is good
    96         %disp(sprintf('%s','   checking time stepping...'));
    97         %[back,dt,time]=TimeStepping(md,solution,dt,time);
    98         %if back,
    99         %       continue;
    100         %end
    101 
    102         %update time and counter
    103         time=time+dt;
    104         n=n+1;
    105 end
    106 
    107 %Load results onto model
    108 results=struct('step',[],'time',[],'vx',[],'vy',[],'vel',[],'pressure',[],'thickness',[],'surface',[],'bed',[]);
    109 for i=1:length(solution),
    110         results(i).step=solution(i).step;
    111         results(i).time=solution(i).time/yts;
    112         results(i).vx=solution(i).u_g(1:2:end)*yts;
    113         results(i).vy=solution(i).u_g(2:2:end)*yts;
    114         results(i).vel=sqrt(solution(i).u_g(1:2:end).^2+solution(i).u_g(2:2:end).^2)*yts;
    115         results(i).bed=solution(i).b_g;
    116         results(i).surface=solution(i).s_g;
    117         results(i).thickness=solution(i).h_g;
    118 end
    119 if ~isstruct(md.results), md.results=struct(); end
    120 md.results.TransientAnalysis=results;
     32        %stop timing
     33        t2=clock;
     34        displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);

  • issm/trunk/src/m/solutions/jpl/transient3d.m

    r3839 r4126  
    66%   See also: TRANSIENT2D, TRANSIENT
    77
    8 %Build all models related to diagnostic
    9 models.analysis_type=TransientAnalysisEnum(); %needed for processresults
     8        analysis_types=[DiagnosticHorizAnalysisEnum,DiagnosticVertAnalysisEnum,DiagnosticStokesAnalysisEnum,DiagnosticHutterAnalysisEnum,SlopeAnalysisEnum,PrognosticAnalysisEnum,ThermalAnalysisEnum,MeltingAnalysisEnum];
     9        solution_type=Transient3DAnalysisEnum;
    1010
    11 displaystring(md.verbose,'%s',['reading diagnostic horiz model data']);
    12 md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=HorizAnalysisEnum(); models.dh=CreateFemModel(md);
     11        displaystring(md.verbose,'%s',['create finite element model']);
     12        femmodel=NewFemModel(md,solution_type,analysis_types,8);
    1313
    14 displaystring(md.verbose,'\n%s',['reading diagnostic vert model data']);
    15 md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=VertAnalysisEnum(); models.dv=CreateFemModel(md);
     14        %retrieve parameters
     15        verbose=femmodel.parameters.Verbose;
     16        qmu_analysis=femmodel.parameters.QmuAnalysis;
    1617
    17 displaystring(md.verbose,'\n%s',['reading diagnostic stokes model data']);
    18 md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=StokesAnalysisEnum(); models.ds=CreateFemModel(md);
     18        %compute solution
     19        if ~qmu_analysis,
     20                displaystring(verbose,'%s',['call computational core']);
     21                femmodel=transient3d_core(femmodel);
     22               
     23                displaystring(verbose,'%s',['write results']);
     24                OutputResults(femmodel.elements, femmodel.loads, femmodel.nodes, femmodel.vertices, femmodel.materials, femmodel.parameters);
    1925
    20 displaystring(md.verbose,'\n%s',['reading diagnostic hutter model data']);
    21 md.analysis_type=DiagnosticAnalysisEnum(); md.sub_analysis_type=HutterAnalysisEnum(); models.dhu=CreateFemModel(md);
     26        else
     27                %launch dakota driver for diagnostic core solution
     28                Qmu(femmodel);
     29        end
    2230
    23 displaystring(md.verbose,'\n%s',['reading surface and bed slope computation model data']);
    24 md.analysis_type=SlopecomputeAnalysisEnum(); md.sub_analysis_type=NoneAnalysisEnum(); models.sl=CreateFemModel(md);
    25 
    26 displaystring(md.verbose,'%s',['reading prognostic model data']);
    27 md.analysis_type=PrognosticAnalysisEnum(); models.p=CreateFemModel(md);
    28 
    29 %Build all models related to thermal
    30 displaystring(md.verbose,'%s',['reading thermal model data']);
    31 md.analysis_type=ThermalAnalysisEnum(); models.t=CreateFemModel(md);
    32 
    33 displaystring(md.verbose,'%s',['reading melting model data']);
    34 md.analysis_type=MeltingAnalysisEnum(); models.m=CreateFemModel(md);
    35 
    36 
    37 %initialize solution
    38 results=struct('step',[],'time',[],'u_g',[],'p_g',[],'h_g',[],'s_g',[],'b_g',[]);
    39 results.step=1;
    40 results.time=0;
    41 results.u_g=models.dh.parameters.u_g;
    42 %results.u_g=models.dh.parameters.u_g(dofsetgen([1,2],3,length(models.dh.parameters.u_g)));
    43 results.p_g=models.p.parameters.p_g;
    44 results.h_g=models.p.parameters.h_g;
    45 results.s_g=models.p.parameters.s_g;
    46 results.b_g=models.p.parameters.b_g;
    47 results.t_g=models.p.parameters.t_g;
    48 results.m_g=models.p.parameters.m_g;
    49 results.a_g=models.p.parameters.a_g;
    50 
    51 %initialize inputs
    52 displaystring(md.verbose,'\n%s',['setup inputs...']);
    53 inputs=inputlist;
    54 inputs=add(inputs,'velocity',results.u_g,'doublevec',3,models.p.parameters.NumberOfNodes);
    55 inputs=add(inputs,'melting',results.m_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    56 inputs=add(inputs,'accumulation',results.a_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    57 inputs=add(inputs,'dt',models.p.parameters.Dt*models.p.parameters.Yts,'double');
    58 
    59 % figure out number of dof: just for information purposes.
    60 md.dof=modelsize(models);
    61 
    62 %first time step is given by model.
    63 dt=models.p.parameters.Dt;
    64 finaltime=models.p.parameters.Ndt;
    65 time=dt;
    66 yts=models.p.parameters.Yts;
    67 n=1; %counter
    68 
    69 while  time<finaltime+dt, %make sure we run up to finaltime.
    70 
    71         displaystring(md.verbose,'\n%s%g%s%g%s%g\n','time [yr]: ',time,'    iteration number: ',n,'/',floor(finaltime/dt));
    72 
    73         results(n+1).step=n+1;
    74         results(n+1).time=time;
    75 
    76         %update inputs
    77         inputs=add(inputs,'thickness',results(n).h_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    78         inputs=add(inputs,'surface',results(n).s_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    79         inputs=add(inputs,'bed',results(n).b_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    80         inputs=add(inputs,'velocity',results(n).u_g,'doublevec',3,models.p.parameters.NumberOfNodes);
    81         inputs=add(inputs,'pressure',results(n).p_g,'doublevec',1,models.p.parameters.NumberOfNodes);
    82         inputs=add(inputs,'temperature',results(n).t_g,'doublevec',1,models.t.parameters.NumberOfNodes);
    83 
    84         %Deal with temperature first
    85         displaystring(md.verbose,'\n%s',['    computing temperatures...']);
    86         [results(n+1).t_g models.t.loads melting_offset]=thermal_core_nonlinear(models.t,ThermalAnalysisEnum(),TransientAnalysisEnum());
    87         inputs=add(inputs,'temperature',results(n+1).t_g,'doublevec',1,models.t.parameters.NumberOfNodes);
    88        
    89         displaystring(md.verbose,'\n%s',['    computing melting...']);
    90         inputs=add(inputs,'melting_offset',melting_offset,'double');
    91         results(n+1).m_g=diagnostic_core_linear(models.m,MeltingAnalysisEnum(),TransientAnalysisEnum());
    92 
    93         %Compute depth averaged temperature
    94         temperature_average=FieldDepthAverage(models.t.elements,models.t.nodes,models.t.vertices,models.t.loads,models.t.materials,models.t.parameters,results(n+1).t_g,'temperature');
    95         inputs=add(inputs,'temperature_average',temperature_average,'doublevec',1,models.t.parameters.NumberOfNodes);
    96 
    97         %Deal with velocities.
    98         rawresults=diagnostic_core(models);
    99         results(n+1).u_g=rawresults.u_g; results(n+1).p_g=rawresults.p_g;
    100 
    101         %compute new thickness
    102         displaystring(md.verbose,'\n%s',['    computing new thickness...']);
    103         inputs=add(inputs,'velocity',results(n+1).u_g,'doublevec',3,models.p.parameters.NumberOfNodes);
    104         rawresults=prognostic_core(models,PrognosticAnalysisEnum(),NoneAnalysisEnum());
    105         new_thickness=rawresults.h_g;
    106 
    107         %update surface and bed using the new thickness
    108         displaystring(md.verbose,'\n%s',['    updating geometry...']);
    109         [new_thickness,new_bed,new_surface]=UpdateGeometry(models.p.elements,models.p.nodes,models.p.vertices,models.p.loads,models.p.materials,models.p.parameters,new_thickness,results(n).b_g,results(n).s_g);
    110 
    111         %Record bed surface and thickness in the results
    112         results(n+1).h_g=new_thickness;
    113         results(n+1).b_g=new_bed;
    114         results(n+1).s_g=new_surface;
    115 
    116         %figure out if time stepping is good
    117         %displaystring(md.verbose,'   checking time stepping...'));
    118         %[back,dt,time]=TimeStepping(md,results,dt,time);
    119         %if back,
    120         %       continue;
    121         %end
    122 
    123         %update node positions
    124         displaystring(md.verbose,'\n%s',['    updating node positions...']);
    125         models.dh.vertices=UpdateVertexPositions(models.dh.vertices,new_bed,new_thickness);
    126         models.dv.vertices=UpdateVertexPositions(models.dv.vertices,new_bed,new_thickness);
    127         models.ds.vertices=UpdateVertexPositions(models.ds.vertices,new_bed,new_thickness);
    128         models.sl.vertices=UpdateVertexPositions(models.sl.vertices,new_bed,new_thickness);
    129         models.p.vertices =UpdateVertexPositions(models.p.vertices,new_bed,new_thickness);
    130         models.t.vertices =UpdateVertexPositions(models.t.vertices,new_bed,new_thickness);
    131         models.m.vertices =UpdateVertexPositions(models.m.vertices,new_bed,new_thickness);
    132        
    133         %update time and counter
    134         time=time+dt;
    135         n=n+1;
    136 end
    137 
    138 %process results
    139 if ~isstruct(md.results), md.results=struct(); end
    140 md.results.transient=processresults(models, results);
     31        %stop timing
     32        t2=clock;
     33        displaystring(md.verbose,'\n%s\n',['solution converged in ' num2str(etime(t2,t1)) ' seconds']);
Note: See TracChangeset for help on using the changeset viewer.