Changeset 5820


Ignore:
Timestamp:
09/15/10 09:52:40 (15 years ago)
Author:
seroussi
Message:

modifications in solver_stokescoupling

Location:
issm/trunk/src/m
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • issm/trunk/src/m/solutions/diagnostic_core.m

    r5579 r5820  
    4747        end
    4848                       
    49         if ismacayealpattyn,
     49        if xor(ismacayealpattyn,isstokes), %if macayealpattyn and stokes, we must use stokescoupling
    5050
    5151                displaystring(verbose,'\n%s',['computing horizontal velocities...']);
     
    5454        end
    5555       
     56        if ismacayealpattyn & isstokes,
    5657
    57         if dim==3 & ~isstokes,
     58                displaystring(verbose,'\n%s',['computing coupling velocities for macayealpattyn and stokes...']);
     59                femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
     60                femmodel=solver_stokescoupling_nonlinear(femmodel,conserve_loads);
     61        end
     62
     63        if dim==3 & (ismacayealpattyn | ishutter),
    5864       
    5965                displaystring(verbose,'\n%s',['computing vertical velocities...']);
    6066                femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
    6167                femmodel=solver_linear(femmodel);
    62         end
    63 
    64         if isstokes,
    65 
    66                 displaystring(verbose,'\n%s',['computing stokes velocities and pressure...']);
    67                 femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
    68                 femmodel=solver_diagnostic_nonlinear(femmodel,conserve_loads);
    6968        end
    7069
  • issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m

    r5812 r5820  
    1010        %First get ug=ug_horiz+ug_vert
    1111        femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
    12         ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters);
     12        ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
    1313
    1414        femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum);
    15         ug_vert=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters);
    16 
    17         ug=ug_horiz+ug_vert;
     15        ug_vert=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
    1816
    1917        while(~converged),
    2018
    21                 old_ug=ug;
    22                
    2319                %First compute the horizontal velocity
    2420                femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum);
    2521
    26                 %keep a copy of loads for now
    27                 loads=femmodel.loads;
    28 
    2922                %Start non-linear iteration using input velocity:
    30                 ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, loads, femmodel.materials, femmodel.parameters);
     23                ug_horiz=GetSolutionFromInputs(femmodel.elements, femmodel.nodes, femmodel.vertices, femmodel.loads, femmodel.materials, femmodel.parameters);
    3124                uf_horiz=Reducevectorgtof( ug_horiz, femmodel.nodesets,femmodel.parameters);
    3225
    3326                %Update the solution to make sure that vx and vxold are similar
    34                 [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug_horiz);
     27                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);
    3528
    3629                %save pointer to old velocity
     
    3831                old_uf_horiz=uf_horiz;
    3932
    40                 [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters);
     33                [K_gg,K_ff,K_fs,p_g,p_f,kmax]=SystemMatrices(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters);
    4134
    4235                [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters);
     
    4740                ug_horiz= Mergesolutionfromftog( uf_horiz, femmodel.ys, femmodel.nodesets,femmodel.parameters);
    4841
    49                 [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,loads,femmodel.materials,femmodel.parameters,ug_horiz);
     42                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_horiz);
    5043
    5144                %Figure out if convergence have been reached
    52                 converged=convergence(K_ff,p_f,uf,old_uf,femmodel.parameters);
     45                converged=convergence(K_ff,p_f,uf_horiz,old_uf_horiz,femmodel.parameters);
    5346
    5447                %Then compute vertical velocity
     
    6659                [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_vert);
    6760
    68                 %Finally sum the two velocities and check the convergence
    69                 ug=ug_horiz+ug_vert;
    70                 converged=convergence([],[],ug,old_ug,femmodel.parameters);
    71 
    7261        end
    7362
Note: See TracChangeset for help on using the changeset viewer.