Changeset 25427

Timestamp:

08/18/20 17:42:43 (5 years ago)

Author:

adhikari

Message:

CHG: added rotational feedbacks to original SESAW solver

Location:

issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver

Files:

• SESAWslr.m (modified) (4 diffs)
• greensfunctions.m (modified) (3 diffs)
• latelonge.m (modified) (1 diff)
• rotationalfeedback.m (added)

issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/SESAWslr.m

@@ -1 @@
-function [eust,rsl] = SESAWslr(index,greens,para) 
+function [eust,rsl] = SESAWslr(index,lat,long,greens,para) 
 %SESAWslr :: computes GRD slr due to applied surface loads based on SESAW method. 
 % 
@@ -5 +5 @@
 %
 %   Usage:
-%      [eus,rsl]=greensfunctions(md.mesh.elements,greens,para) 
+%      [eus,rsl]=greensfunctions(md.mesh.elements,md.mesh.lat,md.mesh.long,greens,para) 
 %      
 %      greens.Grigid = Gravitational Green's function for the rigid planet. 
@@ -18 +18 @@
 %      para.loads_density = land loads (ice or freshwater) density [kg/m^3] 
 %      para.rel_tol = relative tolerance for iterative SLR solver 
+% 
 %      para.solidearth = rheological model for solid Earth: 'rigid' or 'elastic' 
-%      para.isgrd_R = Do we want to account for Rotational feedback? 0 or 1. 
+% 
+%      para.rotational.flag = Do we want to accoutn for rotational feedback? 1 (yes) or 0 (no) 
+%      para.rotational.earth_radius = mean radius of the Earth [m] 
+%      para.rotational.load_love_k2 = degree 2 load love number k  
+%      para.rotational.tide_love_k2 = degree 2 tide love number k  
+%      para.rotational.tide_love_h2 = degree 2 tide love number h 
+%      para.rotational.tide_love_k2secular = degree 2 secular tide love number k  
+%      para.rotational.moi_p = polar moment of inertia [kg m^2]  
+%      para.rotational.moi_e = mean equatorial moment of inertia [kg m^2]  
+%      para.rotational.omega = mean rotational velocity of earth [rad per second]  
 
 area_element = para.area_element; 
@@ -29 +39 @@
 rel_tol = para.rel_tol; 
 
+rot_flag = para.rotational.flag; 
+
 % total Green's function for elastic earth, i.e. (1+k_l-h_l)
 if strcmpi(para.solidearth,'rigid') 
         Galpha = greens.Grigid; 
 elseif strcmpi(para.solidearth,'elastic') 
-        Galpha = greens.Grigid + greens.Gelast- greens.Uelast;          
+        Galpha = greens.Grigid + greens.Gelast - greens.Uelast;         
 else
         error(['Unknown solidearth model:' para.solidearth, '; Should be rigid or elastic']); 
 end
 
+if rot_flag 
+        % set up parameters for rotationalfeedback function. 
+        para_rot = para.rotational; 
+        para_rot.area_element =  para.area_element; 
+        para_rot.loads_element = para.loads_element 
+        para_rot.loads_density = para.loads_density;  
+end
+
+% densitity ratios 
 density_o_e = rho_o/rho_e; 
 density_l_e = rho_l/rho_e; 
 density_l_o = rho_l/rho_o; 
 
+% ocean and earth's surface areas: 
+ocean_area = sum(ocean_element.*area_element); 
+earth_area = sum(area_element); 
+
 % eustatic term 
-eust=-density_l_o*sum(loads_element.*area_element)/sum(ocean_element.*area_element); 
+eust = -density_l_o*sum(loads_element.*area_element)/ocean_area; 
 
-term1=3*density_l_e.*sum(bsxfun(@times,Galpha,(loads_element.*area_element)'),2)./sum(area_element); 
-func3=mean(term1(index),2).*ocean_element;
-term3=sum(func3.*area_element)./sum(ocean_element.*area_element); 
+term1 = 3*density_l_e.*sum(bsxfun(@times,Galpha,(loads_element.*area_element)'),2)./earth_area; 
+if rot_flag
+        term1 = term1 + rotationalfeedback(index,lat,long,para_rot); 
+end
+func3 = mean(term1(index),2).*ocean_element;
+term3 = sum(func3.*area_element)./ocean_area; 
 
 % computation of sea level change 
-rsl=eust+term1-term3;
+rsl = eust+term1-term3;
 
-norm_diff = 10;                 p = 0; 
-
+norm_diff = 10; p = 0; 
 while norm_diff > rel_tol
         norm_old = sqrt(sum(rsl.^2)); 
+        %
+        term2 = 3*density_o_e.*sum(bsxfun(@times,Galpha,(mean(rsl(index),2).*ocean_element.*area_element)'),2)./earth_area; 
+        if rot_flag
+                para_rot.loads_element = mean(rsl(index),2).*ocean_element; 
+                para_rot.loads_density = para.ocean_density; 
+                term2 = term2 + rotationalfeedback(index,lat,long,para_rot); 
+        end
+        func4 = mean(term2(index),2).*ocean_element;
+        term4 = sum(func4.*area_element)./ocean_area; 
         % 
-        term2=3*density_o_e.*sum(bsxfun(@times,Galpha,(mean(rsl(index),2).*ocean_element.*area_element)'),2)./sum(area_element); 
-        func4=mean(term2(index),2).*ocean_element;
-        term4=sum(func4.*area_element)./sum(ocean_element.*area_element); 
-        % 
-        rsl=eust+term1+term2-term3-term4; 
+        rsl = eust+term1+term2-term3-term4; 
         norm_new  = sqrt(sum(rsl.^2));
         norm_diff = abs(norm_new-norm_old)./norm_old;

issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/greensfunctions.m

@@ -1 @@
-function [Grigid,Gelastic,Uelastic] = greensfunctions(md_mesh,love_numbers) 
+function [Grigid,Gelastic,Uelastic] = greensfunctions(mesh_elements,mesh_lat,mesh_long,love_numbers) 
 %greensfunctions :: computes Greens funtions for rigid and elastic Earth. 
 %
 %   Usage:
-%      [Grigid,Gelastic,Uelastic]=greensfunctions(md_mesh,love_numbers)
+%      [Grigid,Gelastic,Uelastic]=greensfunctions(md.mesh.elements,md.mesh.lat,md.mesh.long,md.solidearth.lovenumbers); 
 
-nv=md_mesh.numberofvertices; 
-ne=md_mesh.numberofelements; 
+ne = length(mesh_elements); 
+nv = length(mesh_lat); 
 
 % compute lat,long at elemental centroids. 
-[late,longe] = latelonge(md_mesh); 
+[late,longe] = latelonge(mesh_elements,mesh_lat,mesh_long); 
 
 % love numbers. 
@@ -36 +36 @@
 qq=1; 
 for j=1:nv 
-        phi1(:,1)=md_mesh.lat(j)./180.*pi;      lambda1(:,1)=md_mesh.long(j)./180.*pi; % size #vertices 
+        phi1(:,1)=mesh_lat(j)./180.*pi; lambda1(:,1)=mesh_long(j)./180.*pi; % size #vertices 
         delPhi=abs(phi2-phi1);                                  delLambda=abs(lambda2-lambda1);
         alpha=2.*asin(sqrt(sin(delPhi./2).^2+cos(phi1).*cos(phi2).*sin(delLambda./2).^2)); 
@@ -43 +43 @@
         Uelastic(j,:)=interp1(xx,elast_loveH,cos(alpha)); 
         if (j==500*qq)
-                display([num2str(j),' of ', num2str(md_mesh.numberofvertices),' vertices done!']);
+                display([num2str(j),' of ', num2str(nv),' vertices done!']);
                 qq=qq+1; 
         end 

issm/trunk-jpl/src/m/contrib/adhikari/adhikari2016GMD_slr_solver/latelonge.m

@@ -1 @@
-function [late,longe] = latelonge(md_mesh) 
+function [late,longe] = latelonge(mesh_elements,mesh_lat,mesh_long) 
 %latelonge :: computes lat,long at elemental centroids on spherical surface
 %
 %   Usage:
-%      [late,longe]=latelonge(md.mesh)
+%      [late,longe]=latelonge(md.mesh.elements,md.mesh.lat,md.mesh.long) 
 
-nv = md_mesh.numberofvertices; 
-if length(md_mesh.lat)~=nv | length(md_mesh.long)~=nv 
-        error('lat,long not defined properly at vertices.');
-end
+ne = length(mesh_elements); 
+nv = length(mesh_lat); 
 
 % lat -> [0,180]; long -> [0,360] to compute centroids 
-lat=90-md_mesh.lat;             lon=md_mesh.long; 
+lat=90-mesh_lat;                lon=mesh_long; 
 lon(lon<0)=180+(180+lon(lon<0)); 
 
-ax_0=lat(md_mesh.elements(:,1)); ay_0=lon(md_mesh.elements(:,1)); 
-bx_0=lat(md_mesh.elements(:,2)); by_0=lon(md_mesh.elements(:,2)); 
-cx_0=lat(md_mesh.elements(:,3)); cy_0=lon(md_mesh.elements(:,3)); 
+ax_0=lat(mesh_elements(:,1)); ay_0=lon(mesh_elements(:,1)); 
+bx_0=lat(mesh_elements(:,2)); by_0=lon(mesh_elements(:,2)); 
+cx_0=lat(mesh_elements(:,3)); cy_0=lon(mesh_elements(:,3)); 
 % find whether long is 0 or 360! This is important to compute centroids as well as elemental area 
-for ii=1:md_mesh.numberofelements
+for ii=1:ne 
         if (min([ay_0(ii),by_0(ii),cy_0(ii)])==0 && max([ay_0(ii),by_0(ii),cy_0(ii)])>180)
                 if ay_0(ii)==0