Changeset 130

Timestamp:

04/29/09 14:54:53 (16 years ago)

Author:

Mathieu Morlighem

Message:

Initialize viscosity with current velocity in MacAyeal

File:

• issm/trunk/src/m/solutions/macayeal/diagnostic.m (modified) (5 diffs)

issm/trunk/src/m/solutions/macayeal/diagnostic.m

@@ -25 +25 @@
 x=md.x;
 y=md.y;
-
 index=md.elements;index=sort(index,2); %necessary 
 nods=md.numberofgrids;
@@ -47 +46 @@
 drag_type=md.drag_type;
 drag=md.drag;
-
 criterion_rel=md.eps_rel;
 criterion_abs=md.eps_abs;
@@ -53 +51 @@
 B=md.B; B_bar=(B(index(:,1))+B(index(:,2))+B(index(:,3)))/3;
 glen_coeff=md.n;
+
+%initialize velocities if any
+if (~isnan(md.vx) & ~isnan(md.vy)),
+        u=md.vx/yts; v=md.vy/yts;
+        velocity_is_present=1;
+else
+        velocity_is_present=0;
+end
 
 %average of p and q over the grids (size nel->nods)
@@ -266 +272 @@
 
         %Compute viscosity (as in ICE and CIELO)
-        if convergence_count==1;
+        if (convergence_count==1),
                 %Initialize viscosity
-                nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff);
-        elseif convergence_count==2,
+                if ~velocity_is_present;
+                        nu_bar=viscosity(index,nel,alpha,beta,{},{},B_bar,glen_coeff);
+                else
+                        nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff);
+                        nu_bar(find(nu_bar<=0))=-nu_bar(find(nu_bar<=0));
+                        u_old=u; v_old=v;
+                end
+        elseif (convergence_count==2),
                 nu_bar_oldvalue=nu_bar;
                 nu_bar=viscosity(index,nel,alpha,beta,u,v,B_bar,glen_coeff);
@@ -421 +433 @@
 
         %Test for direct shooting convergence
-        if convergence_count>1,
+        if (convergence_count>1 | velocity_is_present),
 
                 ug=[u_old;v_old];