Changeset 3873


Ignore:
Timestamp:
05/21/10 11:04:58 (15 years ago)
Author:
seroussi
Message:

fixed vertical velocity

Location:
issm/trunk/src/m/solutions/jpl
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/m/solutions/jpl/ModelUpdateInputsFromVector.m

    r3853 r3873  
    1 function ModelUpdateInputsFromVector(models, vector, enum, typeenum);
     1function models=ModelUpdateInputsFromVector(models, vector, enum, typeenum);
    22%MODELUPDATEINPUTSFROMVECTOR - Update inputs using a vector
    33%
     
    2626        if isstruct(model), %there is an analysis_type model
    2727                [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);
     28                models.(field)=model;
    2829        end
    2930
  • issm/trunk/src/m/solutions/jpl/diagnostic_core.m

    r3853 r3873  
    66%
    77
    8 %recover fem models
    9 m_dh=model.dh;
    10 m_dv=model.dv;
    11 m_ds=model.ds;
    12 m_dhu=model.dhu;
    13 m_sl=model.sl;
    14 
    158%recover parameters common to all solutions
    16 verbose=m_dhu.parameters.Verbose;
    17 dim=m_dh.parameters.Dim;
    18 ishutter=m_dhu.parameters.IsHutter;
    19 ismacayealpattyn=m_dh.parameters.IsMacAyealPattyn;
    20 isstokes=m_ds.parameters.IsStokes;
    21 numrifts=m_dhu.parameters.NumRifts;
    22 qmu_analysis=m_dh.parameters.QmuAnalysis;
     9verbose=model.dhu.parameters.Verbose;
     10dim=model.dh.parameters.Dim;
     11ishutter=model.dhu.parameters.IsHutter;
     12ismacayealpattyn=model.dh.parameters.IsMacAyealPattyn;
     13isstokes=model.ds.parameters.IsStokes;
     14numrifts=model.dhu.parameters.NumRifts;
     15qmu_analysis=model.dh.parameters.QmuAnalysis;
    2316
    2417       
    2518%for qmu analysis, be sure the velocity input we are starting from  is the one in the parameters:
    2619if qmu_analysis,
    27         ModelUpdateInputsFromVector(model,m_dh.vx,VxEnum,VertexEnum);
    28         ModelUpdateInputsFromVector(model,m_dh.vy,VyEnum,VertexEnum);
    29         ModelUpdateInputsFromVector(model,m_dh.vz,VzEnum,VertexEnum);
     20        model=ModelUpdateInputsFromVector(model,m_dh.vx,VxEnum,VertexEnum);
     21        model=ModelUpdateInputsFromVector(model,m_dh.vy,VyEnum,VertexEnum);
     22        model=ModelUpdateInputsFromVector(model,m_dh.vz,VzEnum,VertexEnum);
    3023end
    3124
    3225%Compute slopes:
    33 [surfaceslopex,surfaceslopey]=slope_core(m_sl,SurfaceAnalysisEnum);
    34 [bedslopex,bedslopey]=slope_core(m_sl,BedAnalysisEnum);
     26[surfaceslopex,surfaceslopey]=slope_core(model.sl,SurfaceAnalysisEnum);
     27[bedslopex,bedslopey]=slope_core(model.sl,BedAnalysisEnum);
    3528
    3629%Update:
    37 ModelUpdateInputsFromVector(model,surfaceslopex,SurfaceSlopexEnum,VertexEnum);
    38 ModelUpdateInputsFromVector(model,surfaceslopey,SurfaceSlopeyEnum,VertexEnum);
    39 ModelUpdateInputsFromVector(model,bedslopex,BedSlopexEnum,VertexEnum);
    40 ModelUpdateInputsFromVector(model,bedslopey,BedSlopeyEnum,VertexEnum);
     30model=ModelUpdateInputsFromVector(model,surfaceslopex,SurfaceSlopexEnum,VertexEnum);
     31model=ModelUpdateInputsFromVector(model,surfaceslopey,SurfaceSlopeyEnum,VertexEnum);
     32model=ModelUpdateInputsFromVector(model,bedslopex,BedSlopexEnum,VertexEnum);
     33model=ModelUpdateInputsFromVector(model,bedslopey,BedSlopeyEnum,VertexEnum);
    4134
    4235if ishutter,
    4336
    4437        displaystring(verbose,'\n%s',['computing hutter velocities...']);
    45         u_g=diagnostic_core_linear(m_dhu,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
     38        u_g=diagnostic_core_linear(model.dhu,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
    4639
    4740        displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']);
    48         p_g=ComputePressure(m_dhu.elements,m_dhu.nodes,m_dhu.vertices,m_dhu.loads,m_dhu.materials,m_dhu.parameters,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
     41        p_g=ComputePressure(model.dhu.elements,model.dhu.nodes,model.dhu.vertices,model.dhu.loads,model.dhu.materials,model.dhu.parameters,DiagnosticAnalysisEnum(),HutterAnalysisEnum());
    4942
    5043        displaystring(verbose,'\n%s',['update boundary conditions for macyeal pattyn using hutter results...']);
    5144        if ismacayealpattyn,
    52                 m_dh.y_g=u_g;
    53                 [m_dh.ys m_dh.ys0]=Reducevectorgtos(m_dh.y_g,m_dh.nodesets);
     45                model.dh.y_g=u_g;
     46                [model.dh.ys model.dh.ys0]=Reducevectorgtos(model.dh.y_g,model.dh.nodesets);
    5447        end
    5548
     
    5952
    6053        displaystring(verbose,'\n%s',['computing horizontal velocities...']);
    61         [u_g m_dh.loads]=diagnostic_core_nonlinear(m_dh,m_dh.loads,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
     54        [u_g model.dh.loads]=diagnostic_core_nonlinear(model.dh,model.dh.loads,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
    6255
    6356        if dim==2,
    6457                displaystring(verbose,'\n%s',['computing pressure according to MacAyeal...']);
    65                 p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
     58                p_g=ComputePressure(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
    6659        end
    6760end
     
    7063
    7164        displaystring(verbose,'\n%s',['extruding horizontal velocities...']);
    72         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);
     65        u_g_horiz=FieldExtrude(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,u_g,'velocity',1);
    7366
    74         [vx,vy]=SplitSolutionVector(u_g_horiz,m_dh.parameters.NumberOfNodes,m_dh.parameters.NumberOfDofsPerNode);
    75         ModelUpdateInputsFromVector(model,vx,VxEnum,VertexEnum);
    76         ModelUpdateInputsFromVector(model,vy,VyEnum,VertexEnum);
     67        [vx,vy]=SplitSolutionVector(u_g_horiz,model.dh.parameters.NumberOfNodes,model.dh.parameters.NumberOfDofsPerNode);
     68        model=ModelUpdateInputsFromVector(model,vx,VxEnum,VertexEnum);
     69        model=ModelUpdateInputsFromVector(model,vy,VyEnum,VertexEnum);
    7770               
    7871        displaystring(verbose,'\n%s',['computing vertical velocities...']);
    79         u_g_vert=diagnostic_core_linear(m_dv,DiagnosticAnalysisEnum(),VertAnalysisEnum());
    80         ModelUpdateInputsFromVector(model,u_g_vert,VzEnum,VertexEnum);
     72        u_g_vert=diagnostic_core_linear(model.dv,DiagnosticAnalysisEnum(),VertAnalysisEnum());
     73        model=ModelUpdateInputsFromVector(model,u_g_vert,VzEnum,VertexEnum);
    8174
    8275        displaystring(verbose,'\n%s',['computing pressure according to Pattyn...']);
    83         p_g=ComputePressure(m_dh.elements,m_dh.nodes,m_dh.vertices,m_dh.loads,m_dh.materials,m_dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
    84         ModelUpdateInputsFromVector(model,p_g,PressureEnum,VertexEnum);
     76        p_g=ComputePressure(model.dh.elements,model.dh.nodes,model.dh.vertices,model.dh.loads,model.dh.materials,model.dh.parameters,DiagnosticAnalysisEnum(),HorizAnalysisEnum());
     77        model=ModelUpdateInputsFromVector(model,p_g,PressureEnum,VertexEnum);
     78        u_g=zeros(model.dh.parameters.NumberOfNodes*3,1); %%%%%%%%%%%%%%%%%%%%%%% NEED TO BE CLEANED
     79        u_g(1:3:end)=vx;
     80        u_g(2:3:end)=vy;
     81        u_g(3:3:end)=u_g_vert;
    8582       
    8683        if isstokes,
     
    9087
    9188                displaystring(verbose,'\n%s',['update boundary conditions for stokes using velocities previously computed...']);
    92                 m_ds.y_g=zeros(m_ds.nodesets.gsize,1);
    93                 m_ds.y_g(dofsetgen([1,2],4,m_ds.nodesets.gsize))=u_g;
    94                 m_ds.y_g(dofsetgen([3],4,m_ds.nodesets.gsize))=u_g_vert;
    95                 [m_ds.ys m_ds.ys0]=Reducevectorgtos(m_ds.y_g,m_ds.nodesets);
     89                model.ds.y_g=zeros(model.ds.nodesets.gsize,1);
     90                model.ds.y_g(dofsetgen([1,2],4,model.ds.nodesets.gsize))=u_g;
     91                model.ds.y_g(dofsetgen([3],4,model.ds.nodesets.gsize))=u_g_vert;
     92                [model.ds.ys model.ds.ys0]=Reducevectorgtos(model.ds.y_g,model.ds.nodesets);
    9693
    9794                displaystring(verbose,'\n%s',['computing stokes velocities and pressure ...']);
    98                 u_g=diagnostic_core_nonlinear(m_ds,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
     95                u_g=diagnostic_core_nonlinear(model.ds,DiagnosticAnalysisEnum(),StokesAnalysisEnum());
    9996       
    10097                %"decondition" pressure
    101                 p_g=u_g(4:4:end)*m_dh.parameters.stokesreconditioning;
     98                p_g=u_g(4:4:end)*model.dh.parameters.stokesreconditioning;
    10299        end
    103100end
     
    109106
    110107if numrifts,
    111         results.riftproperties=OutputRifts(m_dh.loads,m_dh.parameters);
     108        results.riftproperties=OutputRifts(model.dh.loads,model.dh.parameters);
    112109end
  • issm/trunk/src/m/solutions/jpl/diagnostic_core_linear.m

    r3839 r3873  
    88        m.parameters.Kflag=1; m.parameters.Pflag=1;
    99
    10         %Update inputs in datasets
    11         [m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters]=UpdateFromInputs(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters);
    12        
    1310        %system matrices
    1411        [K_gg, p_g]=SystemMatrices(m.elements,m.nodes,m.vertices,m.loads,m.materials,m.parameters,analysis_type,sub_analysis_type);
     
    2118        %Reduce load from g size to f size
    2219        [p_f] = Reduceloadfromgtof( p_g, m.Gmn, K_fs, m.ys, m.nodesets);
    23        
     20
    2421        %Solve 
    2522        [u_f]=Solver(K_ff,p_f,[],m.parameters);
Note: See TracChangeset for help on using the changeset viewer.