Changeset 5820
- Timestamp:
- 09/15/10 09:52:40 (15 years ago)
- Location:
- issm/trunk/src/m
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
issm/trunk/src/m/solutions/diagnostic_core.m
r5579 r5820 47 47 end 48 48 49 if ismacayealpattyn,49 if xor(ismacayealpattyn,isstokes), %if macayealpattyn and stokes, we must use stokescoupling 50 50 51 51 displaystring(verbose,'\n%s',['computing horizontal velocities...']); … … 54 54 end 55 55 56 if ismacayealpattyn & isstokes, 56 57 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), 58 64 59 65 displaystring(verbose,'\n%s',['computing vertical velocities...']); 60 66 femmodel=SetCurrentConfiguration(femmodel,DiagnosticVertAnalysisEnum); 61 67 femmodel=solver_linear(femmodel); 62 end63 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);69 68 end 70 69 -
issm/trunk/src/m/solvers/solver_stokescoupling_nonlinear.m
r5812 r5820 10 10 %First get ug=ug_horiz+ug_vert 11 11 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); 13 13 14 14 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); 18 16 19 17 while(~converged), 20 18 21 old_ug=ug;22 23 19 %First compute the horizontal velocity 24 20 femmodel=SetCurrentConfiguration(femmodel,DiagnosticHorizAnalysisEnum); 25 21 26 %keep a copy of loads for now27 loads=femmodel.loads;28 29 22 %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); 31 24 uf_horiz=Reducevectorgtof( ug_horiz, femmodel.nodesets,femmodel.parameters); 32 25 33 26 %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); 35 28 36 29 %save pointer to old velocity … … 38 31 old_uf_horiz=uf_horiz; 39 32 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); 41 34 42 35 [K_ff, K_fs] = Reducematrixfromgtof( K_gg, femmodel.nodesets,femmodel.parameters); … … 47 40 ug_horiz= Mergesolutionfromftog( uf_horiz, femmodel.ys, femmodel.nodesets,femmodel.parameters); 48 41 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); 50 43 51 44 %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); 53 46 54 47 %Then compute vertical velocity … … 66 59 [femmodel.elements,femmodel.materials]=InputUpdateFromSolution(femmodel.elements,femmodel.nodes,femmodel.vertices,femmodel.loads,femmodel.materials,femmodel.parameters,ug_vert); 67 60 68 %Finally sum the two velocities and check the convergence69 ug=ug_horiz+ug_vert;70 converged=convergence([],[],ug,old_ug,femmodel.parameters);71 72 61 end 73 62
Note:
See TracChangeset
for help on using the changeset viewer.