Skip to content

Commit

Permalink
less and less is working, must be close to a breakthrough
Browse files Browse the repository at this point in the history
  • Loading branch information
@GHilmarG
GHilmarG committed Oct 24, 2023
1 parent a60c470 commit 6530a60
Show file tree
Hide file tree
Showing 3 changed files with 146 additions and 103 deletions.
2 changes: 1 addition & 1 deletion HelmholtzEquation.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,7 +12,7 @@
% If RHS is given as input, c and d are not used.
%
% The in-homogeneous Helmholtz equation with variable coefficents in two spatial
% dimentions is
% dimensions is
%
% $$ a(x,y) f(x,y) - \nabla \cdot (b(x,y) \nabla f(x,y)) = c(x,y) - \nabla \cdot \nabla d(x,y) $$
%
Expand Down
13 changes: 7 additions & 6 deletions PlotResultsFromThicknessInversion.m
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,25 +10,25 @@ function PlotResultsFromThicknessInversion(CtrlVar,MUA,F,BCs,Priors,Meas,hest,ht



UaPlots(CtrlVar,MUA,F,hest,FigureTitle="hest") ;
UaPlots(CtrlVar,MUA,F,hest,FigureTitle="hest",CreateNewFigure=true) ;
title("h estimated")
hold on ; PlotGroundingLines();



if ~isempty(Priors.h)
UaPlots(CtrlVar,MUA,F,Priors.h,FigureTitle="hprior") ; title("h prior")
UaPlots(CtrlVar,MUA,F,Priors.h,FigureTitle="hprior",CreateNewFigure=true) ; title("h prior")
hold on ; PlotGroundingLines();
end

if ~isempty(Meas.h)

UaPlots(CtrlVar,MUA,F,Meas.h,FigureTitle="hmeas") ;
UaPlots(CtrlVar,MUA,F,Meas.h,FigureTitle="hmeas",CreateNewFigure=true) ;
title("$h$ measured",Interpreter="latex")
hold on ; PlotGroundingLines();


UaPlots(CtrlVar,MUA,F,hest-Meas.h,FigureTitle="hest-hmeas") ;
UaPlots(CtrlVar,MUA,F,hest-Meas.h,FigureTitle="hest-hmeas",CreateNewFigure=true) ;
title("$h$: estimated - measured ice thickness",Interpreter="latex")
hold on ; PlotGroundingLines();

Expand All @@ -46,7 +46,7 @@ function PlotResultsFromThicknessInversion(CtrlVar,MUA,F,BCs,Priors,Meas,hest,ht
UaPlots(CtrlVar,MUA,F,htrue,FigureTitle="h true") ; title("$h$ true",Interpreter="latex")
axis tight

FindOrCreateFigure("compare")
FigCompare=FindOrCreateFigure("compare") ; clf(FigCompare)
yyaxis left
plot(htrue,hest,'.')
hold on
Expand All @@ -61,7 +61,8 @@ function PlotResultsFromThicknessInversion(CtrlVar,MUA,F,BCs,Priors,Meas,hest,ht
title(sprintf("R2=%f",lm.Rsquared.Ordinary));
xlabel("h true") ;

UaPlots(CtrlVar,MUA,F,hest-htrue,FigureTitle="hest-htrue") ; title("hest-htrue")
UaPlots(CtrlVar,MUA,F,hest-htrue,FigureTitle="hest-htrue",CreateNewFigure=true) ;
title("hest-htrue")

end

Expand Down
234 changes: 138 additions & 96 deletions html/HelmholtzEquation.html
View file
Original file line number Diff line number Diff line change
@@ -1,103 +1,143 @@

<!DOCTYPE html
PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<!--
<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
<head>
<META http-equiv="Content-Type" content="text/html; charset=UTF-8">
<!--
This HTML was auto-generated from MATLAB code.
To make changes, update the MATLAB code and republish this document.
--><title>HelmholtzEquation</title><meta name="generator" content="MATLAB 9.11"><link rel="schema.DC" href="http://purl.org/dc/elements/1.1/"><meta name="DC.date" content="2021-11-17"><meta name="DC.source" content="HelmholtzEquation.m"><style type="text/css">
html,body,div,span,applet,object,iframe,h1,h2,h3,h4,h5,h6,p,blockquote,pre,a,abbr,acronym,address,big,cite,code,del,dfn,em,font,img,ins,kbd,q,s,samp,small,strike,strong,tt,var,b,u,i,center,dl,dt,dd,ol,ul,li,fieldset,form,label,legend,table,caption,tbody,tfoot,thead,tr,th,td{margin:0;padding:0;border:0;outline:0;font-size:100%;vertical-align:baseline;background:transparent}body{line-height:1}ol,ul{list-style:none}blockquote,q{quotes:none}blockquote:before,blockquote:after,q:before,q:after{content:'';content:none}:focus{outine:0}ins{text-decoration:none}del{text-decoration:line-through}table{border-collapse:collapse;border-spacing:0}

html { min-height:100%; margin-bottom:1px; }
html body { height:100%; margin:0px; font-family:Arial, Helvetica, sans-serif; font-size:10px; color:#000; line-height:140%; background:#fff none; overflow-y:scroll; }
html body td { vertical-align:top; text-align:left; }

h1 { padding:0px; margin:0px 0px 25px; font-family:Arial, Helvetica, sans-serif; font-size:1.5em; color:#d55000; line-height:100%; font-weight:normal; }
h2 { padding:0px; margin:0px 0px 8px; font-family:Arial, Helvetica, sans-serif; font-size:1.2em; color:#000; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }
h3 { padding:0px; margin:0px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:1.1em; color:#000; font-weight:bold; line-height:140%; }

a { color:#005fce; text-decoration:none; }
a:hover { color:#005fce; text-decoration:underline; }
a:visited { color:#004aa0; text-decoration:none; }

p { padding:0px; margin:0px 0px 20px; }
img { padding:0px; margin:0px 0px 20px; border:none; }
p img, pre img, tt img, li img, h1 img, h2 img { margin-bottom:0px; }

ul { padding:0px; margin:0px 0px 20px 23px; list-style:square; }
ul li { padding:0px; margin:0px 0px 7px 0px; }
ul li ul { padding:5px 0px 0px; margin:0px 0px 7px 23px; }
ul li ol li { list-style:decimal; }
ol { padding:0px; margin:0px 0px 20px 0px; list-style:decimal; }
ol li { padding:0px; margin:0px 0px 7px 23px; list-style-type:decimal; }
ol li ol { padding:5px 0px 0px; margin:0px 0px 7px 0px; }
ol li ol li { list-style-type:lower-alpha; }
ol li ul { padding-top:7px; }
ol li ul li { list-style:square; }

.content { font-size:1.2em; line-height:140%; padding: 20px; }

pre, code { font-size:12px; }
tt { font-size: 1.2em; }
pre { margin:0px 0px 20px; }
pre.codeinput { padding:10px; border:1px solid #d3d3d3; background:#f7f7f7; }
pre.codeoutput { padding:10px 11px; margin:0px 0px 20px; color:#4c4c4c; }
pre.error { color:red; }

@media print { pre.codeinput, pre.codeoutput { word-wrap:break-word; width:100%; } }

span.keyword { color:#0000FF }
span.comment { color:#228B22 }
span.string { color:#A020F0 }
span.untermstring { color:#B20000 }
span.syscmd { color:#B28C00 }
span.typesection { color:#A0522D }

.footer { width:auto; padding:10px 0px; margin:25px 0px 0px; border-top:1px dotted #878787; font-size:0.8em; line-height:140%; font-style:italic; color:#878787; text-align:left; float:none; }
.footer p { margin:0px; }
.footer a { color:#878787; }
.footer a:hover { color:#878787; text-decoration:underline; }
.footer a:visited { color:#878787; }

table th { padding:7px 5px; text-align:left; vertical-align:middle; border: 1px solid #d6d4d4; font-weight:bold; }
table td { padding:7px 5px; text-align:left; vertical-align:top; border:1px solid #d6d4d4; }





</style></head><body><div class="content"><pre class="codeinput"><span class="keyword">function</span> [UserVar,f,lambda,HEmatrix,HErhs]=HelmholtzEquation(UserVar,CtrlVar,MUA,a,b,c,d,RHS)
</pre><p>Solves the in-homogeneous Helmholtz equation with variable coefficients in two dimentions:</p><p><img src="HelmholtzEquation_eq01471109207232534691.png" alt="$$ a(x,y) f(x,y) - \nabla \cdot (b(x,y) \nabla f(x,y)) = c(x,y) - \nabla \cdot \nabla d(x,y) $$" style="width:385px;height:15px;"></p><p>Also possible to specify the right-hand-side directly through</p><pre class="language-matlab">RHS
</pre><p>If RHS is given as input, c and d are not used.</p><p>The in-homogeneous Helmholtz equation with variable coefficents in two spatial dimentions is</p><p><img src="HelmholtzEquation_eq09497764449799483310.png" alt="$$ a(x,y) f(x,y) - \nabla \cdot (b(x,y) \nabla f(x,y)) = c(x,y) - \nabla \cdot \nabla d(x,y) $$" style="width:385px;height:15px;"></p><p>which we can also write as</p><p><img src="HelmholtzEquation_eq01950210437778155504.png" alt="$$ a ( f - \tilde{f} ) - \nabla \cdot (b \nabla (f-\tilde{f})) = 0 $$" style="width:199px;height:18px;"></p><p>where <img src="HelmholtzEquation_eq06444168545079738416.png" alt="$c= -a \tilde{f}$" style="width:54px;height:17px;"> and <img src="HelmholtzEquation_eq05397933670334748401.png" alt="$d= -b \tilde{f}$" style="width:54px;height:17px;">, and <img src="HelmholtzEquation_eq02187778057748804630.png" alt="$\tilde{f}$" style="width:9px;height:17px;"> is a given function</p><p>Examples:</p><p>Smooth a given field over a FE mesh:</p><pre> load('PIG-TWG-RestartFile.mat') ; CtrlVar=CtrlVarInRestartFile;
L=1e3 ; % Smoothing length scale
[UserVar,SmoothedField]=HelmholtzEquation([],CtrlVar,MUA,1,L^2,F.B,0);</pre><pre> figure(1) ; PlotMeshScalarVariable(CtrlVarInRestartFile,MUA,SmoothedField) ; title(' Smoothed field') ; xlabel('x (km)') ; ylabel('y (km)')
figure(2) ; PlotMeshScalarVariable(CtrlVarInRestartFile,MUA,F.B) ; title(' Original field') ; xlabel('x (km)') ; ylabel('y (km)')
figure(3) ; PlotMeshScalarVariable(CtrlVarInRestartFile,MUA,SmoothedField-F.B) ; title(' Smoothed-Original') ; xlabel('x (km)') ; ylabel('y (km)')</pre><pre class="codeinput">narginchk(7,8)


[UserVar,HEmatrix,HErhs]=HelmholtzEquationAssembly(UserVar,CtrlVar,MUA,a,b,c,d);
L=[] ; Lrhs=[] ; lambda=[]; f=[] ;

<span class="keyword">if</span> nargin==8 &amp;&amp; ~isempty(RHS)
HErhs=RHS;
<span class="keyword">end</span>

[f,lambda]=solveKApe(HEmatrix,L,HErhs,Lrhs,f,lambda,CtrlVar);
f=full(f);

<span class="comment">%</span>
<span class="comment">% MLC=BCs2MLC(CtrlVar,MUA,BCsTracer);</span>
<span class="comment">% L=MLC.hL ; Lrhs=MLC.hRhs ; lambda=Lrhs*0;</span>
<span class="comment">% [c1,lambda]=solveKApe(kv,L,rh,Lrhs,c0,lambda,CtrlVar);</span>
<span class="comment">% c1=full(c1);</span>
</pre><pre class="codeinput"><span class="keyword">end</span>
</pre><p class="footer"><br><a href="https://www.mathworks.com/products/matlab/">Published with MATLAB&reg; R2021b</a><br></p></div><!--
-->
<title>HelmholtzEquation</title>
<meta name="generator" content="MATLAB 23.2">
<link rel="schema.DC" href="http://purl.org/dc/elements/1.1/">
<meta name="DC.date" content="2023-10-23">
<meta name="DC.source" content="HelmholtzEquation.m">
<style type="text/css">
html,body,div,span,applet,object,iframe,h1,h2,h3,h4,h5,h6,p,blockquote,pre,a,abbr,acronym,address,big,cite,code,del,dfn,em,font,img,ins,kbd,q,s,samp,small,strike,strong,tt,var,b,u,i,center,dl,dt,dd,ol,ul,li,fieldset,form,label,legend,table,caption,tbody,tfoot,thead,tr,th,td{margin:0;padding:0;border:0;outline:0;font-size:100%;vertical-align:baseline;background:transparent}body{line-height:1}ol,ul{list-style:none}blockquote,q{quotes:none}blockquote:before,blockquote:after,q:before,q:after{content:'';content:none}:focus{outine:0}ins{text-decoration:none}del{text-decoration:line-through}table{border-collapse:collapse;border-spacing:0}

html { min-height:100%; margin-bottom:1px; }
html body { height:100%; margin:0px; font-family:Arial, Helvetica, sans-serif; font-size:10px; color:#000; line-height:140%; background:#fff none; overflow-y:scroll; }
html body td { vertical-align:top; text-align:left; }

h1 { padding:0px; margin:0px 0px 25px; font-family:Arial, Helvetica, sans-serif; font-size:1.5em; color:#d55000; line-height:100%; font-weight:normal; }
h2 { padding:0px; margin:0px 0px 8px; font-family:Arial, Helvetica, sans-serif; font-size:1.2em; color:#000; font-weight:bold; line-height:140%; border-bottom:1px solid #d6d4d4; display:block; }
h3 { padding:0px; margin:0px 0px 5px; font-family:Arial, Helvetica, sans-serif; font-size:1.1em; color:#000; font-weight:bold; line-height:140%; }

a { color:#005fce; text-decoration:none; }
a:hover { color:#005fce; text-decoration:underline; }
a:visited { color:#004aa0; text-decoration:none; }

p { padding:0px; margin:0px 0px 20px; }
img { padding:0px; margin:0px 0px 20px; border:none; }
p img, pre img, tt img, li img, h1 img, h2 img { margin-bottom:0px; }

ul { padding:0px; margin:0px 0px 20px 23px; list-style:square; }
ul li { padding:0px; margin:0px 0px 7px 0px; }
ul li ul { padding:5px 0px 0px; margin:0px 0px 7px 23px; }
ul li ol li { list-style:decimal; }
ol { padding:0px; margin:0px 0px 20px 0px; list-style:decimal; }
ol li { padding:0px; margin:0px 0px 7px 23px; list-style-type:decimal; }
ol li ol { padding:5px 0px 0px; margin:0px 0px 7px 0px; }
ol li ol li { list-style-type:lower-alpha; }
ol li ul { padding-top:7px; }
ol li ul li { list-style:square; }

.content { font-size:1.2em; line-height:140%; padding: 20px; }

pre, code { font-size:12px; }
tt { font-size: 1.2em; }
pre { margin:0px 0px 20px; }
pre.codeinput { padding:10px; border:1px solid #d3d3d3; background:#f7f7f7; }
pre.codeoutput { padding:10px 11px; margin:0px 0px 20px; color:#4c4c4c; }
pre.error { color:red; }

@media print { pre.codeinput, pre.codeoutput { word-wrap:break-word; width:100%; } }

span.keyword { color:#0000FF }
span.comment { color:#228B22 }
span.string { color:#A020F0 }
span.untermstring { color:#B20000 }
span.syscmd { color:#B28C00 }
span.typesection { color:#A0522D }

.footer { width:auto; padding:10px 0px; margin:25px 0px 0px; border-top:1px dotted #878787; font-size:0.8em; line-height:140%; font-style:italic; color:#878787; text-align:left; float:none; }
.footer p { margin:0px; }
.footer a { color:#878787; }
.footer a:hover { color:#878787; text-decoration:underline; }
.footer a:visited { color:#878787; }

table th { padding:7px 5px; text-align:left; vertical-align:middle; border: 1px solid #d6d4d4; font-weight:bold; }
table td { padding:7px 5px; text-align:left; vertical-align:top; border:1px solid #d6d4d4; }





</style>
</head>
<body>
<div class="content">
<pre class="codeinput">
<span class="keyword">function</span> [UserVar,f,lambda,HEmatrix,HErhs]=HelmholtzEquation(UserVar,CtrlVar,MUA,a,b,c,d,RHS)
</pre>
<p>Solves the in-homogeneous Helmholtz equation with variable coefficients in two dimensions:</p>
<p>
<img src="HelmholtzEquation_eq01471109207232534691.png" alt="$$ a(x,y) f(x,y) - \nabla \cdot (b(x,y) \nabla f(x,y)) = c(x,y) - \nabla \cdot \nabla d(x,y) $$" style="width:462px;height:18px;"></p>
<p>Also possible to specify the right-hand-side directly through</p>
<pre class="language-matlab">RHS
</pre>
<p>If RHS is given as input, c and d are not used.</p>
<p>The in-homogeneous Helmholtz equation with variable coefficents in two spatial dimensions is</p>
<p>
<img src="HelmholtzEquation_eq09497764449799483310.png" alt="$$ a(x,y) f(x,y) - \nabla \cdot (b(x,y) \nabla f(x,y)) = c(x,y) - \nabla \cdot \nabla d(x,y) $$" style="width:462px;height:18px;"></p>
<p>which we can also write as</p>
<p>
<img src="HelmholtzEquation_eq01950210437778155504.png" alt="$$ a ( f - \tilde{f} ) - \nabla \cdot (b \nabla (f-\tilde{f})) = 0 $$" style="width:238px;height:22px;"></p>
<p>where <img src="HelmholtzEquation_eq06444168545079738416.png" alt="$c= -a \tilde{f}$" style="width:65px;height:21px;"> and <img src="HelmholtzEquation_eq05397933670334748401.png" alt="$d= -b \tilde{f}$" style="width:65px;height:21px;">, and <img src="HelmholtzEquation_eq02187778057748804630.png" alt="$\tilde{f}$" style="width:10px;height:21px;"> is a given function</p>
<p>Examples:</p>
<p>Smooth a given field over a FE mesh:</p>
<pre> load('PIG-TWG-RestartFile.mat') ; CtrlVar=CtrlVarInRestartFile;
L=1e3 ; % Smoothing length scale
[UserVar,SmoothedField]=HelmholtzEquation([],CtrlVar,MUA,1,L^2,F.B,0);</pre>
<pre> figure(1) ; PlotMeshScalarVariable(CtrlVarInRestartFile,MUA,SmoothedField) ; title(' Smoothed field') ; xlabel('x (km)') ; ylabel('y (km)')
figure(2) ; PlotMeshScalarVariable(CtrlVarInRestartFile,MUA,F.B) ; title(' Original field') ; xlabel('x (km)') ; ylabel('y (km)')
figure(3) ; PlotMeshScalarVariable(CtrlVarInRestartFile,MUA,SmoothedField-F.B) ; title(' Smoothed-Original') ; xlabel('x (km)') ; ylabel('y (km)')</pre>
<pre class="codeinput">narginchk(7,8)


[UserVar,HEmatrix,HErhs]=HelmholtzEquationAssembly(UserVar,CtrlVar,MUA,a,b,c,d);
L=[] ; Lrhs=[] ; lambda=[]; f=[] ;

<span class="keyword">if</span> nargin==8 &amp;&amp; ~isempty(RHS)
HErhs=RHS;
<span class="keyword">end</span>

[f,lambda]=solveKApe(HEmatrix,L,HErhs,Lrhs,f,lambda,CtrlVar);
f=full(f);

<span class="comment">%</span>
<span class="comment">% MLC=BCs2MLC(CtrlVar,MUA,BCsTracer);</span>
<span class="comment">% L=MLC.hL ; Lrhs=MLC.hRhs ; lambda=Lrhs*0;</span>
<span class="comment">% [c1,lambda]=solveKApe(kv,L,rh,Lrhs,c0,lambda,CtrlVar);</span>
<span class="comment">% c1=full(c1);</span>
</pre>
<pre class="codeoutput error">Error using HelmholtzEquation
Not enough input arguments.
</pre>
<pre class="codeinput">
<span class="keyword">end</span>
</pre>
<p class="footer">
<br>
<a href="https://www.mathworks.com/products/matlab/">Published with MATLAB&reg; R2023b</a>
<br>
</p>
</div>
<!--
##### SOURCE BEGIN #####
function [UserVar,f,lambda,HEmatrix,HErhs]=HelmholtzEquation(UserVar,CtrlVar,MUA,a,b,c,d,RHS)
%%
% Solves the in-homogeneous Helmholtz equation with variable coefficients in two dimentions:
% Solves the in-homogeneous Helmholtz equation with variable coefficients in two dimensions:
%
% $$ a(x,y) f(x,y) - \nabla \cdot (b(x,y) \nabla f(x,y)) = c(x,y) - \nabla \cdot \nabla d(x,y) $$
%
Expand All @@ -108,7 +148,7 @@
% If RHS is given as input, c and d are not used.
%
% The in-homogeneous Helmholtz equation with variable coefficents in two spatial
% dimentions is
% dimensions is
%
% $$ a(x,y) f(x,y) - \nabla \cdot (b(x,y) \nabla f(x,y)) = c(x,y) - \nabla \cdot \nabla d(x,y) $$
%
Expand Down Expand Up @@ -159,4 +199,6 @@
##### SOURCE END #####
--></body></html>
-->
</body>
</html>

0 comments on commit 6530a60

Please sign in to comment.