Mercurial > repos > public > sbplib

comparison +sbp/+implementations/intOpAWW_orders_8to8_ratio2to1.m @ 608:c923fe6197ff feature/interpolation

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Add inerpolation operator classes and implementations for MC and AWW.
author Martin Almquist <malmquist@stanford.edu>
date Sat, 14 Oct 2017 22:32:25 -0700
parents
children
comparison
equal deleted inserted replaced
607:0546de4b31a2 608:c923fe6197ff
1 function [IC2F,IF2C,Hc,Hf] = IntOp_orders_8to8_ratio2to1(mc,hc,ACC)
2
3 % ACC is a string.
4 % ACC = 'C2F' creates IC2F with one order of accuracy higher than IF2C.
5 % ACC = 'F2C' creates IF2C with one order of accuracy higher than IC2F.
6 ratio = 2;
7 mf = ratio*mc-1;
8 hf = hc/ratio;
9
10 switch ACC
11 case 'F2C'
12 [stencil_F2C,BC_F2C,HcU,HfU] = ...
13 sbp.implementations.intOpAWW_orders_8to8_ratio_2to1_accC2F4_accF2C5;
14 case 'C2F'
15 [stencil_F2C,BC_F2C,HcU,HfU] = ...
16 sbp.implementations.intOpAWW_orders_8to8_ratio_2to1_accC2F5_accF2C4;
17 end
18
19 stencil_width = length(stencil_F2C);
20 stencil_hw = (stencil_width-1)/2;
21 [BC_rows,BC_cols] = size(BC_F2C);
22
23 %%% Norm matrices %%%
24 Hc = speye(mc,mc);
25 HcUm = length(HcU);
26 Hc(1:HcUm,1:HcUm) = spdiags(HcU',0,HcUm,HcUm);
27 Hc(mc-HcUm+1:mc,mc-HcUm+1:mc) = spdiags(rot90(HcU',2),0,HcUm,HcUm);
28 Hc = Hc*hc;
29
30 Hf = speye(mf,mf);
31 HfUm = length(HfU);
32 Hf(1:HfUm,1:HfUm) = spdiags(HfU',0,HfUm,HfUm);
33 Hf(mf-length(HfU)+1:mf,mf-length(HfU)+1:mf) = spdiags(rot90(HfU',2),0,HfUm,HfUm);
34 Hf = Hf*hf;
35 %%%%%%%%%%%%%%%%%%%%%%
36
37 %%% Create IF2C from stencil and BC
38 IF2C = sparse(mc,mf);
39 for i = BC_rows+1 : mc-BC_rows
40 IF2C(i,ratio*i-1+(-stencil_hw:stencil_hw)) = stencil_F2C; %#ok<SPRIX>
41 end
42 IF2C(1:BC_rows,1:BC_cols) = BC_F2C;
43 IF2C(end-BC_rows+1:end,end-BC_cols+1:end) = rot90(BC_F2C,2);
44 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45
46 %%% Create IC2F using symmetry condition %%%%
47 IC2F = Hf\IF2C.'*Hc;