Mercurial > repos > public > sbplib
changeset 743:f4595f14d696 feature/utux2D
Change schemes to work for special coefficients.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Mon, 21 May 2018 23:11:34 -0700 |
| parents | 07f8311374c6 |
| children | 0be9b4d6737b |
| files | +scheme/LaplaceCurvilinear.m +scheme/Schrodinger2d.m +scheme/Utux2D.m |
| diffstat | 3 files changed, 29 insertions(+), 12 deletions(-) [+] |
line wrap: on
line diff
diff -r 07f8311374c6 -r f4595f14d696 +scheme/LaplaceCurvilinear.m --- a/+scheme/LaplaceCurvilinear.m Wed Mar 14 21:01:34 2018 -0700 +++ b/+scheme/LaplaceCurvilinear.m Mon May 21 23:11:34 2018 -0700 @@ -57,6 +57,10 @@ end % assert(isa(g, 'grid.Curvilinear')) + if isa(a, 'function_handle') + a = grid.evalOn(g, a); + a = spdiag(a); + end m = g.size(); m_u = m(1);
diff -r 07f8311374c6 -r f4595f14d696 +scheme/Schrodinger2d.m --- a/+scheme/Schrodinger2d.m Wed Mar 14 21:01:34 2018 -0700 +++ b/+scheme/Schrodinger2d.m Mon May 21 23:11:34 2018 -0700 @@ -46,6 +46,10 @@ dim = 2; assert(isa(g, 'grid.Cartesian')) + if isa(a, 'function_handle') + a = grid.evalOn(g, a); + a = spdiag(a); + end m = g.size(); m_tot = g.N(); @@ -178,7 +182,7 @@ Hi = obj.Hi; H_gamma = obj.H_boundary{j}; - a = obj.a; + a = e{j}'*obj.a*e{j}; switch type
diff -r 07f8311374c6 -r f4595f14d696 +scheme/Utux2D.m --- a/+scheme/Utux2D.m Wed Mar 14 21:01:34 2018 -0700 +++ b/+scheme/Utux2D.m Mon May 21 23:11:34 2018 -0700 @@ -7,6 +7,7 @@ v0 % Initial data a % Wave speed a = [a1, a2]; + % Can either be a constant vector or a cell array of function handles. H % Discrete norm H_x, H_y % Norms in the x and y directions @@ -45,7 +46,15 @@ default_arg('coupling_type','upwind'); default_arg('a',1/sqrt(2)*[1, 1]); default_arg('opSet',@sbp.D2Standard); + assert(isa(g, 'grid.Cartesian')) + if iscell(a) + a1 = grid.evalOn(g, a{1}); + a2 = grid.evalOn(g, a{2}); + a = {spdiag(a1), spdiag(a2)}; + else + a = {a(1), a(2)}; + end m = g.size(); m_x = m(1); @@ -96,7 +105,7 @@ obj.coupling_type = coupling_type; obj.interpolation_type = interpolation_type; obj.interpolation_damping = interpolation_damping; - obj.D = -(a(1)*obj.Dx + a(2)*obj.Dy); + obj.D = -(a{1}*obj.Dx + a{2}*obj.Dy); end % Closure functions return the opertors applied to the own domain to close the boundary @@ -112,11 +121,11 @@ sigma = -1; % Scalar penalty parameter switch boundary case {'w','W','west','West'} - tau = sigma*obj.a(1)*obj.e_w*obj.H_y; + tau = sigma*obj.a{1}*obj.e_w*obj.H_y; closure = obj.Hi*tau*obj.e_w'; case {'s','S','south','South'} - tau = sigma*obj.a(2)*obj.e_s*obj.H_x; + tau = sigma*obj.a{2}*obj.e_s*obj.H_x; closure = obj.Hi*tau*obj.e_s'; end penalty = -obj.Hi*tau; @@ -223,35 +232,35 @@ switch boundary case {'w','W','west','West'} - tau = sigma_ds*obj.a(1)*obj.e_w*obj.H_y; + tau = sigma_ds*obj.a{1}*obj.e_w*obj.H_y; closure = obj.Hi*tau*obj.e_w'; penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; - beta = int_damp_ds*obj.a(1)... + beta = int_damp_ds*obj.a{1}... *obj.e_w*obj.H_y; closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_w'; case {'e','E','east','East'} - tau = sigma_us*obj.a(1)*obj.e_e*obj.H_y; + tau = sigma_us*obj.a{1}*obj.e_e*obj.H_y; closure = obj.Hi*tau*obj.e_e'; penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; - beta = int_damp_us*obj.a(1)... + beta = int_damp_us*obj.a{1}... *obj.e_e*obj.H_y; closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_e'; case {'s','S','south','South'} - tau = sigma_ds*obj.a(2)*obj.e_s*obj.H_x; + tau = sigma_ds*obj.a{2}*obj.e_s*obj.H_x; closure = obj.Hi*tau*obj.e_s'; penalty = -obj.Hi*tau*I_neighbour2local_ds*e_neighbour'; - beta = int_damp_ds*obj.a(2)... + beta = int_damp_ds*obj.a{2}... *obj.e_s*obj.H_x; closure = closure + obj.Hi*beta*(I_back_forth_ds - I)*obj.e_s'; case {'n','N','north','North'} - tau = sigma_us*obj.a(2)*obj.e_n*obj.H_x; + tau = sigma_us*obj.a{2}*obj.e_n*obj.H_x; closure = obj.Hi*tau*obj.e_n'; penalty = -obj.Hi*tau*I_neighbour2local_us*e_neighbour'; - beta = int_damp_us*obj.a(2)... + beta = int_damp_us*obj.a{2}... *obj.e_n*obj.H_x; closure = closure + obj.Hi*beta*(I_back_forth_us - I)*obj.e_n'; end