fwd_solve_2p5d_1st_order

 FWD_SOLVE_2P5D_1ST_ORDER: data= fwd_solve_2p5d_1st_order( img)
 Fwd solver for Andy Adler's EIT code
 Input:
    img       = image struct
 Output:
    data = measurements struct
 Options: (to return internal FEM information)
    img.fwd_solve.get_all_meas = 1 (data.volt = all FEM nodes, but not CEM)
    img.fwd_solve.get_all_nodes= 1 (data.volt = all nodes, including CEM)

    img.fwd_solve_2p5d_1st_order.k = [ a .. b ]
        solve, integrating over the range k = a .. b      (default: [0 Inf])
        - provide a single k to get a point solution
        - solve over a reduced range a=0, b=3 for a faster solution
    img.fwd_solve_2p5d_1st_order.method = 'name'
        perform numerical integration using the selected method (default: 'quadv')
        'trapz' - trapezoidal integration across the listed points k
        'quadv' - adaptive quadrature (vectorized), from a to b
        'integral' - adaptive quadrature (matlab2012+), from a to b

0001 function data =fwd_solve_2p5d_1st_order(fwd_model, img)
0002 % FWD_SOLVE_2P5D_1ST_ORDER: data= fwd_solve_2p5d_1st_order( img)
0003 % Fwd solver for Andy Adler's EIT code
0004 % Input:
0005 %    img       = image struct
0006 % Output:
0007 %    data = measurements struct
0008 % Options: (to return internal FEM information)
0009 %    img.fwd_solve.get_all_meas = 1 (data.volt = all FEM nodes, but not CEM)
0010 %    img.fwd_solve.get_all_nodes= 1 (data.volt = all nodes, including CEM)
0011 %
0012 %    img.fwd_solve_2p5d_1st_order.k = [ a .. b ]
0013 %        solve, integrating over the range k = a .. b      (default: [0 Inf])
0014 %        - provide a single k to get a point solution
0015 %        - solve over a reduced range a=0, b=3 for a faster solution
0016 %    img.fwd_solve_2p5d_1st_order.method = 'name'
0017 %        perform numerical integration using the selected method (default: 'quadv')
0018 %        'trapz' - trapezoidal integration across the listed points k
0019 %        'quadv' - adaptive quadrature (vectorized), from a to b
0020 %        'integral' - adaptive quadrature (matlab2012+), from a to b
0021 
0022 % (C) 2016 A Boyle
0023 % License: GPL version 2 or version 3
0024 
0025 % correct input paralemeters if function was called with only img
0026 if ischar(fwd_model) && strcmp(fwd_model,'UNIT_TEST'); do_unit_test; return; end
0027 
0028 if nargin == 1
0029    img= fwd_model;
0030 end
0031 fwd_model= img.fwd_model;
0032 
0033 img = data_mapper(img);
0034 if ~ismember(img.current_params, supported_params)
0035     error('EIDORS:PhysicsNotSupported', '%s does not support %s', ...
0036     'FWD_SOLVE_2P5D_1ST_ORDER',img.current_params);
0037 end
0038 % all calcs use conductivity
0039 img = convert_img_units(img, 'conductivity');
0040 
0041 pp= fwd_model_parameters( fwd_model, 'skip_VOLUME' );
0042 
0043 k = [0 Inf]; try k = img.fwd_model.fwd_solve_2p5d_1st_order.k; end % default to 8 full waves: k pi h / 2 z = 8 * 2 pi
0044 method = 'quadv'; try method = img.fwd_model.fwd_solve_2p5d_1st_order.method; end % default to 8 full waves: k pi h / 2 z = 8 * 2 pi
0045 
0046 img.fwd_model.system_mat = @system_mat_2p5d_1st_order;
0047 gnd = fwd_model.gnd_node;
0048 img.fwd_model.system_mat_2p5d_1st_order.k = 0;
0049 img.fwd_model.system_mat_2p5d_1st_order.factory = 1;
0050 Sf = system_mat_2p5d_1st_order( img ); % returns a function Sf(k) that builds a system matrix for 'k'
0051 if length(k) == 1 % singleton k
0052    v=2*potential_k(Sf(k), pp, gnd); % integral of a delta function
0053 else
0054    switch method
0055       case 'trapz'
0056          % less accurate: trapz
0057          n = 0;
0058          tol = 1e-8;
0059          k(isinf(k)) = 1/sqrt(tol);
0060          vf = zeros(pp.n_elec*pp.n_stim,length(k)); % voltages under electrodes elec x stim, frequency domain
0061          for ki = k
0062             n = n + 1;
0063             vf(:,n) = potential_k(Sf(ki), pp, gnd);
0064          end
0065          v=2/pi*trapz(k,vf,2);
0066       case 'quadv'
0067          % more accurate: adaptive gaussian quadrature
0068          tol = 1e-8;
0069          kend = min(1/sqrt(tol), max(k)); % don't go too far... k=Inf is a singular matrix, stop adjacent to numeric singularity
0070          % quadv is scheduled to be removed from matlab eventually... but it is
0071          % WAY faster than integral with any tolerance configuration I could identify
0072          v=2/pi*quadv(@(kk) potential_k(Sf(kk), pp, gnd), k(1), kend);
0073       case 'integral'
0074          tol = 1e-8;
0075          kend = min(1/sqrt(tol), max(k)); % don't go too far... k=Inf is a singular matrix, stop adjacent to numeric singularity
0076          % the integral solution is about 10x slower (5.47 seconds vs. 0.60 seconds for UNIT_TEST)
0077          % ... I played with AbsTol and RelTol but wasn't able to affect the outcome
0078          v=2/pi*integral(@(kk) potential_k(Sf(kk), pp, gnd), k(1), kend, 'ArrayVAlued', true);
0079       otherwise
0080          error(['unrecognized method: ' method]);
0081    end
0082 end
0083 % up to now, we've been dealing with voltages at the electrodes,
0084 % now we transform to dipole measurements v=(m+ - m-)
0085 data.meas = pp.v2meas * v;
0086 data.time= NaN; % unknown
0087 data.name= 'solved by fwd_solve_2p5d_1st_order';
0088 try; if img.fwd_solve.get_all_meas == 1
0089    data.volt = v(1:pp.n_node,:); % but not on CEM nodes
0090 end; end
0091 try; if img.fwd_solve.get_all_nodes== 1
0092    data.volt = v;                % all, including CEM nodes
0093 end; end
0094 
0095 % v0: a zero initialized array n x m for n nodes and m stimulus
0096 function v = potential_k(S, pp, gnd)
0097    vf_nodes = zeros(pp.n_node, pp.n_stim); % voltages at all nodes x number of stim, frequency domain
0098    idx = 1:size(S,1); % pp.n_node + #CEM elec
0099    idx( gnd ) = [];
0100    idx2 = find(any(pp.N2E));
0101    vf_nodes(idx,:) = left_divide( S(idx,idx), 1/2*pp.QQ(idx,:));
0102    vf_elecs = pp.N2E(:,idx2) * vf_nodes(idx2,:); % electrode voltages, frequency domain
0103    v = vf_elecs(:); 
0104 
0105 function do_unit_test
0106    ne = 16;
0107    isOctave= exist('OCTAVE_VERSION');
0108    % 2D
0109    imdl2 = mk_geophysics_model('h2a', ne);
0110    imdl2.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0111    img2 = mk_image(imdl2.fwd_model, 1);
0112    % we get a 2D solution for klim=0; z = 1/2; h = 1;
0113    img2.fwd_model.fwd_solve_2p5d_1st_order.k = [0 logspace(-2,0,12) 3]; % default
0114    tic
0115    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0116    v2 = fwd_solve(img2);
0117    fprintf('fwd_solve 2D t= %f s\n', toc);
0118    % 2.5D
0119    tic
0120    img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'trapz';
0121    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0122    v2p5t = fwd_solve_2p5d_1st_order(img2);
0123    fprintf('fwd_solve_2p5d_1st_order 2D t= %f s (%d k''s, trapz)\n', toc, length(img2.fwd_model.fwd_solve_2p5d_1st_order.k));
0124    tic
0125    img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'quadv';
0126    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0127    v2p5q = fwd_solve_2p5d_1st_order(img2);
0128    fprintf('fwd_solve_2p5d_1st_order 2D t= %f s (%d k''s, quadv)\n', toc, length(img2.fwd_model.fwd_solve_2p5d_1st_order.k));
0129    if isOctave
0130       fprintf('fwd_solve_2p5d_1st_order 2D t= --- s (-- k''s, integral) SKIPPED [not supported in Octave]\n');
0131    else
0132       tic
0133       img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'integral';
0134       img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0135       v2p5i = fwd_solve_2p5d_1st_order(img2);
0136       fprintf('fwd_solve_2p5d_1st_order 2D t= %f s (%d k''s, integral)\n', toc, length(img2.fwd_model.fwd_solve_2p5d_1st_order.k));
0137    end
0138    tic
0139    img2.fwd_model.fwd_solve_2p5d_1st_order.k = 0;
0140    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0141    v2p5k0 = fwd_solve_2p5d_1st_order(img2);
0142    fprintf('fwd_solve_2p5d_1st_order 2D t= %f s (%d k''s)\n', toc, length(img2.fwd_model.fwd_solve_2p5d_1st_order.k));
0143    % 3D
0144    imdl3 = mk_geophysics_model('h3a', ne);
0145    imdl3.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0146    img3 = mk_image(imdl3.fwd_model, 1);
0147    tic
0148    img3.fwd_model.nodes(1,:) = img3.fwd_model.nodes(1,:) + rand(1,3)*1e-8; % defeat cache
0149    v3 = fwd_solve(img3);
0150    fprintf('fwd_solve 3D t= %f s\n', toc);
0151    % analytic half-space (dipole)
0152    tic
0153    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0154    vh = fwd_solve_halfspace(img2);
0155    fprintf('fwd_solve_halfspace 3D t= %f s\n', toc);
0156    scale=1e3; uvh=unique(round(vh.meas*scale)/scale);
0157    unit_test_cmp('h2a 2p5d values TEST', uvh, ...
0158    [ 0.0060; 0.0080; 0.0110; 0.0160; 0.0320 ], 1/scale);
0159    tol = norm(vh.meas)*0.025; % 2.5% error tolerance
0160    unit_test_cmp('2D   (h2a)                    vs analytic TEST', norm(vh.meas - v2.meas),     0, -tol);
0161    unit_test_cmp('2.5D (h2a + 2.5D, k=0)        vs 2D       TEST', norm(v2.meas - v2p5k0.meas), 0, tol);
0162    unit_test_cmp('2.5D (h2a + 2.5D trapz,2*tol) vs analytic TEST', norm(vh.meas - v2p5t.meas),   0, 2*tol);
0163    unit_test_cmp('2.5D (h2a + 2.5D quadv**)     vs analytic TEST', norm(vh.meas - v2p5q.meas),   0, tol);
0164    if ~isOctave
0165       unit_test_cmp('2.5D (h2a + 2.5D integral)    vs analytic TEST', norm(vh.meas - v2p5i.meas),   0, tol);
0166    end
0167    unit_test_cmp('3D   (h3b)                    vs analytic TEST', norm(vh.meas - v3.meas),     0, 2*tol);
0168 pause(2);
0169 clf; h=plot([vh.meas, v2p5q.meas, v3.meas, v2.meas],':');
0170      legend('analytic', 'FEM 2.5D', 'FEM 3D', 'FEM 2D', 'Location','Best');
0171      set(gca,'box','off');
0172      set(h,'LineWidth',2);
0173      set(h(1),'Marker','o','MarkerSize',7);
0174      set(h(2),'Marker','square','MarkerSize',7);
0175      set(h(3),'Marker','diamond','MarkerSize',3);
0176      set(h(4),'Marker','+');
0177 pause(2);
0178 clf; h=plot([vh.meas, v2p5q.meas, v3.meas],':');
0179      legend('analytic', 'FEM 2.5D', 'FEM 3D','Location','Best');
0180      set(gca,'box','off');
0181      set(h,'LineWidth',2);
0182      set(h(1),'Marker','o','MarkerSize',7);
0183      set(h(2),'Marker','square','MarkerSize',7);
0184      set(h(3),'Marker','diamond','MarkerSize',3);