comparison +sbp/+implementations/d2_variable_2.m @ 860:b758d1cf4c8e feature/poroelastic

Add computation of HI*M to D2Variable to make adjoint gradient computation easier.
author Martin Almquist <malmquist@stanford.edu>
date Wed, 24 Oct 2018 16:16:43 -0700
parents ded4156e53e2
children e54c2f54dbfe
comparison
equal deleted inserted replaced
859:4c7532db42cd 860:b758d1cf4c8e
25 25
26 % D1 operator 26 % D1 operator
27 diags = -1:1; 27 diags = -1:1;
28 stencil = [-1/2 0 1/2]; 28 stencil = [-1/2 0 1/2];
29 D1 = stripeMatrix(stencil, diags, m); 29 D1 = stripeMatrix(stencil, diags, m);
30 30
31 D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1; 31 D1(1,1)=-1;D1(1,2)=1;D1(m,m-1)=-1;D1(m,m)=1;
32 D1(m,m-1)=-1;D1(m,m)=1; 32 D1(m,m-1)=-1;D1(m,m)=1;
33 D1=D1/h; 33 D1=D1/h;
34 %Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m'); 34 %Q=H*D1 + 1/2*(e_1*e_1') - 1/2*(e_m*e_m');
35 35
38 38
39 scheme_width = 3; 39 scheme_width = 3;
40 scheme_radius = (scheme_width-1)/2; 40 scheme_radius = (scheme_width-1)/2;
41 r = (1+scheme_radius):(m-scheme_radius); 41 r = (1+scheme_radius):(m-scheme_radius);
42 42
43 function D2 = D2_fun(c) 43 function [D2, B] = D2_fun(c)
44 44
45 Mm1 = -c(r-1)/2 - c(r)/2; 45 Mm1 = -c(r-1)/2 - c(r)/2;
46 M0 = c(r-1)/2 + c(r) + c(r+1)/2; 46 M0 = c(r-1)/2 + c(r) + c(r+1)/2;
47 Mp1 = -c(r)/2 - c(r+1)/2; 47 Mp1 = -c(r)/2 - c(r+1)/2;
48 48
52 M(1:2,1:2)=[c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;]; 52 M(1:2,1:2)=[c(1)/2 + c(2)/2 -c(1)/2 - c(2)/2; -c(1)/2 - c(2)/2 c(1)/2 + c(2) + c(3)/2;];
53 M(m-1:m,m-1:m)=[c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;]; 53 M(m-1:m,m-1:m)=[c(m-2)/2 + c(m-1) + c(m)/2 -c(m-1)/2 - c(m)/2; -c(m-1)/2 - c(m)/2 c(m-1)/2 + c(m)/2;];
54 M=M/h; 54 M=M/h;
55 55
56 D2=HI*(-M-c(1)*e_l*d1_l'+c(m)*e_r*d1_r'); 56 D2=HI*(-M-c(1)*e_l*d1_l'+c(m)*e_r*d1_r');
57 B = HI*M;
57 end 58 end
58 D2 = @D2_fun; 59 D2 = @D2_fun;
60
59 end 61 end