Mercurial > repos > public > sbplib
changeset 741:5a9acf282b34 feature/poroelastic
Add scheme Heat2Dcurvilinear. Neumann and Dirichlet seem to work. Only tested for stretched Cartesian grids though.
| author | Martin Almquist <malmquist@stanford.edu> |
|---|---|
| date | Wed, 09 May 2018 19:29:12 -0700 |
| parents | f4e2a6a2df08 |
| children | 08f3ffe63f48 |
| files | +scheme/Heat2dCurvilinear.m |
| diffstat | 1 files changed, 383 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
diff -r f4e2a6a2df08 -r 5a9acf282b34 +scheme/Heat2dCurvilinear.m --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/+scheme/Heat2dCurvilinear.m Wed May 09 19:29:12 2018 -0700 @@ -0,0 +1,383 @@ +classdef Heat2dCurvilinear < scheme.Scheme + +% Discretizes the Laplacian with variable coefficent, curvilinear, +% in the Heat equation way (i.e., the discretization matrix is not necessarily +% symmetric) +% u_t = div * (kappa * grad u ) +% opSet should be cell array of opSets, one per dimension. This +% is useful if we have periodic BC in one direction. + + properties + m % Number of points in each direction, possibly a vector + h % Grid spacing + + grid + dim + + order % Order of accuracy for the approximation + + % Diagonal matrix for variable coefficients + KAPPA % Variable coefficient + + D % Total operator + D1 % First derivatives + + % Second derivatives + D2_kappa + + H, Hi % Inner products + e_l, e_r + d1_l, d1_r % Normal derivatives at the boundary + alpha % Vector of borrowing constants + + % Boundary inner products + H_boundary_l, H_boundary_r + + % Metric coefficients + b % Cell matrix of size dim x dim + J, Ji + + % Numerical boundary flux operators + flux_l, flux_r + + end + + methods + + function obj = Heat2dCurvilinear(g ,order, kappa_fun, opSet) + default_arg('opSet',{@sbp.D2Variable, @sbp.D2Variable}); + default_arg('kappa_fun', @(x,y) 0*x+1); + dim = 2; + + kappa = grid.evalOn(g, kappa_fun); + m = g.size(); + m_tot = g.N(); + + % 1D operators + ops = cell(dim,1); + for i = 1:dim + ops{i} = opSet{i}(m(i), {0, 1}, order); + end + + I = cell(dim,1); + D1 = cell(dim,1); + D2 = cell(dim,1); + H = cell(dim,1); + Hi = cell(dim,1); + e_l = cell(dim,1); + e_r = cell(dim,1); + d1_l = cell(dim,1); + d1_r = cell(dim,1); + + for i = 1:dim + I{i} = speye(m(i)); + D1{i} = ops{i}.D1; + D2{i} = ops{i}.D2; + H{i} = ops{i}.H; + Hi{i} = ops{i}.HI; + e_l{i} = ops{i}.e_l; + e_r{i} = ops{i}.e_r; + d1_l{i} = ops{i}.d1_l; + d1_r{i} = ops{i}.d1_r; + end + + %====== Assemble full operators ======== + KAPPA = spdiag(kappa); + obj.KAPPA = KAPPA; + + % Allocate + obj.D1 = cell(dim,1); + obj.D2_kappa = cell(dim,1); + obj.e_l = cell(dim,1); + obj.e_r = cell(dim,1); + obj.d1_l = cell(dim,1); + obj.d1_r = cell(dim,1); + + % D1 + obj.D1{1} = kron(D1{1},I{2}); + obj.D1{2} = kron(I{1},D1{2}); + + % -- Metric coefficients ---- + coords = g.points(); + x = coords(:,1); + y = coords(:,2); + + % Use non-periodic difference operators for metric even if opSet is periodic. + xmax = max(ops{1}.x); + ymax = max(ops{2}.x); + opSetMetric{1} = sbp.D2Variable(m(1), {0, xmax}, order); + opSetMetric{2} = sbp.D2Variable(m(2), {0, ymax}, order); + D1Metric{1} = kron(opSetMetric{1}.D1, I{2}); + D1Metric{2} = kron(I{1}, opSetMetric{2}.D1); + + x_xi = D1Metric{1}*x; + x_eta = D1Metric{2}*x; + y_xi = D1Metric{1}*y; + y_eta = D1Metric{2}*y; + + J = x_xi.*y_eta - x_eta.*y_xi; + + b = cell(dim,dim); + b{1,1} = y_eta./J; + b{1,2} = -x_eta./J; + b{2,1} = -y_xi./J; + b{2,2} = x_xi./J; + + % Scale factors for boundary integrals + beta = cell(dim,1); + beta{1} = sqrt(x_eta.^2 + y_eta.^2); + beta{2} = sqrt(x_xi.^2 + y_xi.^2); + + J = spdiag(J); + Ji = inv(J); + for i = 1:dim + beta{i} = spdiag(beta{i}); + for j = 1:dim + b{i,j} = spdiag(b{i,j}); + end + end + obj.J = J; + obj.Ji = Ji; + obj.b = b; + %---------------------------- + + % Boundary operators + obj.e_l{1} = kron(e_l{1},I{2}); + obj.e_l{2} = kron(I{1},e_l{2}); + obj.e_r{1} = kron(e_r{1},I{2}); + obj.e_r{2} = kron(I{1},e_r{2}); + + obj.d1_l{1} = kron(d1_l{1},I{2}); + obj.d1_l{2} = kron(I{1},d1_l{2}); + obj.d1_r{1} = kron(d1_r{1},I{2}); + obj.d1_r{2} = kron(I{1},d1_r{2}); + + % D2 coefficients + kappa_coeff = cell(dim,dim); + for j = 1:dim + obj.D2_kappa{j} = sparse(m_tot,m_tot); + kappa_coeff{j} = sparse(m_tot,1); + for i = 1:dim + kappa_coeff{j} = kappa_coeff{j} + b{i,j}*J*b{i,j}*kappa; + end + end + ind = grid.funcToMatrix(g, 1:m_tot); + + % x-dir + j = 1; + for col = 1:m(2) + D_kappa = D2{1}(kappa_coeff{j}(ind(:,col))); + + p = ind(:,col); + obj.D2_kappa{j}(p,p) = D_kappa; + end + + % y-dir + j = 2; + for row = 1:m(1) + D_kappa = D2{2}(kappa_coeff{j}(ind(row,:))); + + p = ind(row,:); + obj.D2_kappa{j}(p,p) = D_kappa; + end + + % Quadratures + obj.H = kron(H{1},H{2}); + obj.Hi = inv(obj.H); + obj.H_boundary_l = cell(dim,1); + obj.H_boundary_l{1} = obj.e_l{1}'*beta{1}*obj.e_l{1}*H{2}; + obj.H_boundary_l{2} = obj.e_l{2}'*beta{2}*obj.e_l{2}*H{1}; + obj.H_boundary_r = cell(dim,1); + obj.H_boundary_r{1} = obj.e_r{1}'*beta{1}*obj.e_r{1}*H{2}; + obj.H_boundary_r{2} = obj.e_r{2}'*beta{2}*obj.e_r{2}*H{1}; + + %=== Differentiation matrix D (without SAT) === + D2_kappa = obj.D2_kappa; + D1 = obj.D1; + D = sparse(m_tot,m_tot); + + d = @kroneckerDelta; % Kronecker delta + db = @(i,j) 1-d(i,j); % Logical not of Kronecker delta + + % 2nd derivatives + for j = 1:dim + D = D + Ji*D2_kappa{j}; + end + + % Mixed terms + for i = 1:dim + for j = 1:dim + for k = 1:dim + D = D + db(i,j)*Ji*D1{j}*b{i,j}*J*KAPPA*b{i,k}*D1{k}; + end + end + end + obj.D = D; + %=========================================% + + % Normal flux operators for BC. + flux_l = cell(dim,1); + flux_r = cell(dim,1); + + d1_l = obj.d1_l; + d1_r = obj.d1_r; + e_l = obj.e_l; + e_r = obj.e_r; + + % Loop over boundaries + for j = 1:dim + flux_l{j} = sparse(m_tot,m_tot); + flux_r{j} = sparse(m_tot,m_tot); + + % Loop over dummy index + for i = 1:dim + % Loop over dummy index + for k = 1:dim + flux_l{j} = flux_l{j} ... + - beta{j}\b{i,j}*J*KAPPA*b{i,k}*( d(j,k)*e_l{k}*d1_l{k}' + db(j,k)*D1{k} ); + + flux_r{j} = flux_r{j} ... + + beta{j}\b{i,j}*J*KAPPA*b{i,k}*( d(j,k)*e_r{k}*d1_r{k}' + db(j,k)*D1{k} ); + end + + end + end + obj.flux_l = flux_l; + obj.flux_r = flux_r; + + % Misc. + obj.m = m; + obj.h = g.scaling(); + obj.order = order; + obj.grid = g; + obj.dim = dim; + obj.alpha = [ops{1}.borrowing.M.d1, ops{2}.borrowing.M.d1]; + + end + + + % Closure functions return the operators applied to the own domain to close the boundary + % Penalty functions return the operators to force the solution. In the case of an interface it returns the operator applied to the other doamin. + % boundary is a string specifying the boundary e.g. 'l','r' or 'e','w','n','s'. + % type is a string specifying the type of boundary condition. + % data is a function returning the data that should be applied at the boundary. + % neighbour_scheme is an instance of Scheme that should be interfaced to. + % neighbour_boundary is a string specifying which boundary to interface to. + function [closure, penalty] = boundary_condition(obj, boundary, type, symmetric, tuning) + default_arg('type','Neumann'); + default_arg('symmetric', false); + default_arg('tuning',1.2); + + % j is the coordinate direction of the boundary + % nj: outward unit normal component. + % nj = -1 for west, south, bottom boundaries + % nj = 1 for east, north, top boundaries + [j, nj] = obj.get_boundary_number(boundary); + switch nj + case 1 + e = obj.e_r{j}; + flux = obj.flux_r{j}; + H_gamma = obj.H_boundary_r{j}; + case -1 + e = obj.e_l{j}; + flux = obj.flux_l{j}; + H_gamma = obj.H_boundary_l{j}; + end + + Hi = obj.Hi; + Ji = obj.Ji; + KAPPA = obj.KAPPA; + kappa_gamma = e'*KAPPA*e; + h = obj.h(j); + alpha = h*obj.alpha(j); + + switch type + + % Dirichlet boundary condition + case {'D','d','dirichlet','Dirichlet'} + + if ~symmetric + closure = -Ji*Hi*flux'*e*H_gamma*(e' ); + penalty = Ji*Hi*flux'*e*H_gamma; + else + closure = Ji*Hi*flux'*e*H_gamma*(e' )... + -tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma*(e' ) ; + penalty = -Ji*Hi*flux'*e*H_gamma ... + +tuning*2/alpha*Ji*Hi*e*kappa_gamma*H_gamma; + end + + % Normal flux boundary condition + case {'N','n','neumann','Neumann'} + closure = -Ji*Hi*e*H_gamma*(e'*flux ); + penalty = Ji*Hi*e*H_gamma; + + % Unknown boundary condition + otherwise + error('No such boundary condition: type = %s',type); + end + end + + function [closure, penalty] = interface(obj,boundary,neighbour_scheme,neighbour_boundary) + % u denotes the solution in the own domain + % v denotes the solution in the neighbour domain + error('Interface not implemented'); + end + + % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. + function [j, nj] = get_boundary_number(obj, boundary) + + switch boundary + case {'w','W','west','West', 'e', 'E', 'east', 'East'} + j = 1; + case {'s','S','south','South', 'n', 'N', 'north', 'North'} + j = 2; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch boundary + case {'w','W','west','West','s','S','south','South'} + nj = -1; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + nj = 1; + end + end + + % Returns the coordinate number and outward normal component for the boundary specified by the string boundary. + function [return_op] = get_boundary_operator(obj, op, boundary) + + switch boundary + case {'w','W','west','West', 'e', 'E', 'east', 'East'} + j = 1; + case {'s','S','south','South', 'n', 'N', 'north', 'North'} + j = 2; + otherwise + error('No such boundary: boundary = %s',boundary); + end + + switch op + case 'e' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.e_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.e_r{j}; + end + case 'd' + switch boundary + case {'w','W','west','West','s','S','south','South'} + return_op = obj.d1_l{j}; + case {'e', 'E', 'east', 'East','n', 'N', 'north', 'North'} + return_op = obj.d1_r{j}; + end + otherwise + error(['No such operator: operatr = ' op]); + end + + end + + function N = size(obj) + N = prod(obj.m); + end + end +end