Changeset 18271

Timestamp:

07/21/14 14:32:09 (11 years ago)

Author:

schlegel

Message:

CHG: add runme2

Location:

issm/trunk-jpl/examples/Inversion

Files:

• runme.m (modified) (4 diffs)
• runme2.m (added)

issm/trunk-jpl/examples/Inversion/runme.m

@@ -1 @@
-step=3;
+step=1;
 if step==1
         %Generate observation
         md = model;
         md = triangle(md,'DomainOutline.exp',100000);
-    
-        md = setmask(md,'','');
-    
+        md = setmask(md,'all','');
         md = parameterize(md,'Square.par');
         md = setflowequation(md,'SSA','all');
         md.cluster = generic('np',2);
-    
-    md.geometry.surface=md.geometry.surface+100;
-    md.geometry.base = md.geometry.base+100;
-    md.materials.rheology_B(:)=1.8*10^8;
-    md.friction.coefficient(:)=50;
-    md.friction.coefficient(find(md.mesh.x>400000 & md.mesh.x<600000))=10;
-    
         md = solve(md,StressbalanceSolutionEnum);
-        plotmodel(md,'data',md.friction.coefficient,'caxis',[0 100],'figure',1);
+        plotmodel(md,'data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'figure',1);
         plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
         save model1 md
@@ -25 +16 @@
         %Modify rheology, now constant
         loadmodel('model1.mat');
-        md.friction.coefficient(:)=50;
+        md.materials.rheology_B(:) = 1.8*10^8;
 
         %results of previous run are taken as observations
@@ -34 +25 @@
 
         md = solve(md,StressbalanceSolutionEnum);
-        plotmodel(md,'data',md.friction.coefficient,'caxis',[0 100],'figure',1);
+        plotmodel(md,'data',md.materials.rheology_B,'caxis',[ 1.3 1.9]*10^8,'figure',1);
         plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
         save model2 md
@@ -45 +36 @@
         maxsteps = 20;
         md.inversion.iscontrol = 1;
-        md.inversion.control_parameters = {'FrictionCoefficient'};
+        md.inversion.control_parameters = {'MaterialsRheologyBbar'};
         md.inversion.maxsteps = maxsteps;
-        md.inversion.dxmin=10^-6;
-        md.inversion.cost_functions = [103 501];
-        md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,2);
-        md.inversion.cost_functions_coefficients(:,2)=0;
-        md.inversion.min_parameters    = 10^-5*ones(md.mesh.numberofvertices,1);
-        md.inversion.max_parameters    = 100*ones(md.mesh.numberofvertices,1);
+        md.inversion.cost_functions = 101;
+        md.inversion.cost_functions_coefficients = ones(md.mesh.numberofvertices,1);
+        md.inversion.min_parameters    = paterson(273)*ones(md.mesh.numberofvertices,1);
+        md.inversion.max_parameters    = paterson(200)*ones(md.mesh.numberofvertices,1);
 
         %Go solve!
         md.verbose=verbose(0);
         md=solve(md,StressbalanceSolutionEnum);
-        plotmodel(md,'data',md.results.StressbalanceSolution.FrictionCoefficient,'caxis',[0 100],'figure',1);
+        plotmodel(md,'data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'figure',1);
         plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
 end
+if step==4
+        %invert for ice rigidity
+        loadmodel('model2.mat');
+
+        %Set up inversion parameters
+        maxsteps = 20;
+        md.inversion.iscontrol = 1;
+        md.inversion.control_parameters = {'MaterialsRheologyBbar'};
+        md.inversion.maxsteps = maxsteps;
+        md.inversion.cost_functions = [101 502];
+        md.inversion.cost_functions_coefficients      = ones(md.mesh.numberofvertices,1);
+        md.inversion.cost_functions_coefficients(:,2) = 10^-16*ones(md.mesh.numberofvertices,1);
+        md.inversion.min_parameters    = paterson(273)*ones(md.mesh.numberofvertices,1);
+        md.inversion.max_parameters    = paterson(200)*ones(md.mesh.numberofvertices,1);
+
+        %Go solve!
+        md.verbose=verbose(0);
+        md=solve(md,StressbalanceSolutionEnum);
+        plotmodel(md,'data',md.results.StressbalanceSolution.MaterialsRheologyBbar,'caxis',[ 1.3 1.9]*10^8,'figure',1);
+        plotmodel(md,'data',md.results.StressbalanceSolution.Vel,'figure',2);  
+end