Mercurial > repos > public > sbplib
changeset 681:7368affc8f78 feature/poroelastic
Add D2 variable periodic for second order.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Wed, 07 Feb 2018 15:42:50 -0800 |
| parents | cd1a76c38565 |
| children | 50e77b15d841 |
| files | +sbp/+implementations/d2_variable_periodic_2.m +sbp/D2VariablePeriodic.m spdiagVariable.m spdiagsVariablePeriodic.m stripeMatrixPeriodic.m |
| diffstat | 5 files changed, 193 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/+implementations/d2_variable_periodic_2.m Wed Feb 07 15:42:50 2018 -0800 @@ -0,0 +1,51 @@ +function [H, HI, D1, D2, e_l, e_r, d1_l, d1_r] = d2_variable_periodic_2(m,h) + % m = number of unique grid points, i.e. h = L/m; + + BP = 1; + if(m<2*BP) + error(['Operator requires at least ' num2str(2*BP) ' grid points']); + end + + % Norm + Hv = ones(m,1); + Hv = h*Hv; + H = spdiag(Hv, 0); + HI = spdiag(1./Hv, 0); + + + % Dummy boundary operators + e_l = sparse(m,1); + e_r = rot90(e_l, 2); + + d1_l = sparse(m,1); + d1_r = -rot90(d1_l, 2); + + % D1 operator + diags = -1:1; + stencil = [-1/2 0 1/2]; + D1 = stripeMatrixPeriodic(stencil, diags, m); + D1 = D1/h; + + scheme_width = 3; + scheme_radius = (scheme_width-1)/2; + + r = 1:m; + offset = scheme_width; + r = r + offset; + + function D2 = D2_fun(c) + c = [c(end-scheme_width+1:end); c; c(1:scheme_width) ]; + + Mm1 = -c(r-1)/2 - c(r)/2; + M0 = c(r-1)/2 + c(r) + c(r+1)/2; + Mp1 = -c(r)/2 - c(r+1)/2; + + vals = [Mm1,M0,Mp1]; + diags = -scheme_radius : scheme_radius; + M = spdiagsVariablePeriodic(vals,diags); + + M=M/h; + D2=HI*(-M ); + end + D2 = @D2_fun; +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+sbp/D2VariablePeriodic.m Wed Feb 07 15:42:50 2018 -0800 @@ -0,0 +1,75 @@ +classdef D2VariablePeriodic < sbp.OpSet + properties + D1 % SBP operator approximating first derivative + H % Norm matrix + HI % H^-1 + Q % Skew-symmetric matrix + e_l % Left boundary operator + e_r % Right boundary operator + D2 % SBP operator for second derivative + M % Norm matrix, second derivative + d1_l % Left boundary first derivative + d1_r % Right boundary first derivative + m % Number of grid points. + h % Step size + x % grid + borrowing % Struct with borrowing limits for different norm matrices + end + + methods + function obj = D2VariablePeriodic(m,lim,order) + + x_l = lim{1}; + x_r = lim{2}; + L = x_r-x_l; + obj.h = L/m; + x = linspace(x_l,x_r,m+1)'; + obj.x = x(1:end-1); + + switch order + + case 6 + error('Not impl') + + [obj.H, obj.HI, obj.D1, obj.D2, ... + ~, obj.e_l, obj.e_r, ~, ~, ~, ~, ~,... + obj.d1_l, obj.d1_r] = ... + sbp.implementations.d4_variable_periodic_6(m, obj.h); + obj.borrowing.M.d1 = 0.1878; + obj.borrowing.R.delta_D = 0.3696; + % Borrowing e^T*D1 - d1 from R + + case 4 + error('Not impl') + + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_4(m,obj.h); + obj.borrowing.M.d1 = 0.2505765857; + + obj.borrowing.R.delta_D = 0.577587500088313; + % Borrowing e^T*D1 - d1 from R + case 2 + [obj.H, obj.HI, obj.D1, obj.D2, obj.e_l,... + obj.e_r, obj.d1_l, obj.d1_r] = ... + sbp.implementations.d2_variable_periodic_2(m,obj.h); + obj.borrowing.M.d1 = 0.3636363636; + % Borrowing const taken from Virta 2014 + + obj.borrowing.R.delta_D = 1.000000538455350; + % Borrowing e^T*D1 - d1 from R + + otherwise + error('Invalid operator order %d.',order); + end + obj.borrowing.H11 = obj.H(1,1)/obj.h; % First element in H/h, + obj.m = m; + obj.M = []; + end + function str = string(obj) + str = [class(obj) '_' num2str(obj.order)]; + end + end + + +end
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/spdiagVariable.m Wed Feb 07 15:42:50 2018 -0800 @@ -0,0 +1,17 @@ +function A = spdiagVariable(a,i) + default_arg('i',0); + + if isrow(a) + a = a'; + end + + n = length(a)+abs(i); + + if i > 0 + a = [sparse(i,1); a]; + elseif i < 0 + a = [a; sparse(abs(i),1)]; + end + + A = spdiags(a,i,n,n); +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/spdiagsVariablePeriodic.m Wed Feb 07 15:42:50 2018 -0800 @@ -0,0 +1,42 @@ +function A = spdiagsVariablePeriodic(vals,diags) + % Creates an m x m periodic discretization matrix. + % vals - m x ndiags matrix of values + % diags - 1 x ndiags vector of the 'center diagonals' that vals end up on + % vals that are not on main diagonal are going to spill over to + % off-diagonal corners. + + default_arg('diags',0); + + [m, ~] = size(vals); + + A = sparse(m,m); + + for i = 1:length(diags) + + d = diags(i); + a = vals(:,i); + + % Sub-diagonals + if d < 0 + a_bulk = a(1+abs(d):end); + a_corner = a(1:1+abs(d)-1); + corner_diag = m-abs(d); + A = A + spdiagVariable(a_bulk, d); + A = A + spdiagVariable(a_corner, corner_diag); + + % Super-diagonals + elseif d > 0 + a_bulk = a(1:end-d); + a_corner = a(end-d+1,end); + corner_diag = -m + d; + A = A + spdiagVariable(a_bulk, d); + A = A + spdiagVariable(a_corner, corner_diag); + + % Main diagonal + else + A = A + spdiagVariable(a, 0); + end + + end + +end \ No newline at end of file
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/stripeMatrixPeriodic.m Wed Feb 07 15:42:50 2018 -0800 @@ -0,0 +1,8 @@ +% Creates a periodic discretization matrix of size n x n +% with the values of val on the diagonals diag. +% A = stripeMatrix(val,diags,n) +function A = stripeMatrixPeriodic(val,diags,n) + + D = ones(n,1)*val; + A = spdiagsVariablePeriodic(D,diags); +end \ No newline at end of file