fwd_solve_2p5d_1st_order

PURPOSE ^

FWD_SOLVE_2P5D_1ST_ORDER: data= fwd_solve_2p5d_1st_order( img)

SYNOPSIS ^

function data =fwd_solve_2p5d_1st_order(fwd_model, img)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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 v2meas = get_v2meas(pp.n_elec, pp.n_stim, img.fwd_model.stimulation);
0086 data.meas = v2meas' * v;
0087 data.time= NaN; % unknown
0088 data.name= 'solved by fwd_solve_2p5d_1st_order';
0089 try; if img.fwd_solve.get_all_meas == 1
0090    data.volt = v(1:pp.n_node,:); % but not on CEM nodes
0091 end; end
0092 try; if img.fwd_solve.get_all_nodes== 1
0093    data.volt = v;                % all, including CEM nodes
0094 end; end
0095 
0096 % v0: a zero initialized array n x m for n nodes and m stimulus
0097 function v = potential_k(S, pp, gnd)
0098    vf_nodes = zeros(pp.n_node, pp.n_stim); % voltages at all nodes x number of stim, frequency domain
0099    idx = 1:size(S,1); % pp.n_node + #CEM elec
0100    idx( gnd ) = [];
0101    idx2 = find(any(pp.N2E));
0102    vf_nodes(idx,:) = left_divide( S(idx,idx), 1/2*pp.QQ(idx,:));
0103    vf_elecs = pp.N2E(:,idx2) * vf_nodes(idx2,:); % electrode voltages, frequency domain
0104    v = vf_elecs(:); 
0105 
0106 function v2meas = get_v2meas(n_elec,n_stim,stim)
0107     v2meas = sparse(n_elec*n_stim,0);
0108     for i=1:n_stim
0109         meas_pat= stim(i).meas_pattern;
0110         n_meas  = size(meas_pat,1);
0111         v2meas((i-1)*n_elec + 1: i*n_elec,end+(1:n_meas)) = meas_pat';
0112     end
0113 
0114 function do_unit_test
0115    ne = 16;
0116    isOctave= exist('OCTAVE_VERSION');
0117    % 2D
0118    imdl2 = mk_geophysics_model('h2a', ne);
0119    imdl2.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0120    img2 = mk_image(imdl2.fwd_model, 1);
0121    % we get a 2D solution for klim=0; z = 1/2; h = 1;
0122    img2.fwd_model.fwd_solve_2p5d_1st_order.k = [0 logspace(-2,0,12) 3]; % default
0123    tic
0124    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0125    v2 = fwd_solve(img2);
0126    fprintf('fwd_solve 2D t= %f s\n', toc);
0127    % 2.5D
0128    tic
0129    img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'trapz';
0130    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0131    v2p5t = fwd_solve_2p5d_1st_order(img2);
0132    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));
0133    tic
0134    img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'quadv';
0135    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0136    v2p5q = fwd_solve_2p5d_1st_order(img2);
0137    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));
0138    if isOctave
0139       fprintf('fwd_solve_2p5d_1st_order 2D t= --- s (-- k''s, integral) SKIPPED [not supported in Octave]\n');
0140    else
0141       tic
0142       img2.fwd_model.fwd_solve_2p5d_1st_order.method = 'integral';
0143       img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0144       v2p5i = fwd_solve_2p5d_1st_order(img2);
0145       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));
0146    end
0147    tic
0148    img2.fwd_model.fwd_solve_2p5d_1st_order.k = 0;
0149    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0150    v2p5k0 = fwd_solve_2p5d_1st_order(img2);
0151    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));
0152    % 3D
0153    imdl3 = mk_geophysics_model('h3a', ne);
0154    imdl3.fwd_model.stimulation = stim_pattern_geophys(ne, 'Wenner');
0155    img3 = mk_image(imdl3.fwd_model, 1);
0156    tic
0157    img3.fwd_model.nodes(1,:) = img3.fwd_model.nodes(1,:) + rand(1,3)*1e-8; % defeat cache
0158    v3 = fwd_solve(img3);
0159    fprintf('fwd_solve 3D t= %f s\n', toc);
0160    % analytic half-space (dipole)
0161    tic
0162    img2.fwd_model.nodes(1,:) = img2.fwd_model.nodes(1,:) + rand(1,2)*1e-8; % defeat cache
0163    vh = fwd_solve_halfspace(img2);
0164    fprintf('fwd_solve_halfspace 3D t= %f s\n', toc);
0165    scale=1e3; uvh=unique(round(vh.meas*scale)/scale);
0166    unit_test_cmp('h2a 2p5d values TEST', uvh, ...
0167    [ 0.0060; 0.0080; 0.0110; 0.0160; 0.0320 ], 1/scale);
0168    tol = norm(vh.meas)*0.025; % 2.5% error tolerance
0169    unit_test_cmp('2D   (h2a)                    vs analytic TEST', norm(vh.meas - v2.meas),     0, -tol);
0170    unit_test_cmp('2.5D (h2a + 2.5D, k=0)        vs 2D       TEST', norm(v2.meas - v2p5k0.meas), 0, tol);
0171    unit_test_cmp('2.5D (h2a + 2.5D trapz,2*tol) vs analytic TEST', norm(vh.meas - v2p5t.meas),   0, 2*tol);
0172    unit_test_cmp('2.5D (h2a + 2.5D quadv**)     vs analytic TEST', norm(vh.meas - v2p5q.meas),   0, tol);
0173    if ~isOctave
0174       unit_test_cmp('2.5D (h2a + 2.5D integral)    vs analytic TEST', norm(vh.meas - v2p5i.meas),   0, tol);
0175    end
0176    unit_test_cmp('3D   (h3b)                    vs analytic TEST', norm(vh.meas - v3.meas),     0, 2*tol);
0177 pause(2);
0178 clf; h=plot([vh.meas, v2p5q.meas, v3.meas, v2.meas],':');
0179      legend('analytic', 'FEM 2.5D', 'FEM 3D', 'FEM 2D', '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);
0185      set(h(4),'Marker','+');
0186 pause(2);
0187 clf; h=plot([vh.meas, v2p5q.meas, v3.meas],':');
0188      legend('analytic', 'FEM 2.5D', 'FEM 3D','Location','Best');
0189      set(gca,'box','off');
0190      set(h,'LineWidth',2);
0191      set(h(1),'Marker','o','MarkerSize',7);
0192      set(h(2),'Marker','square','MarkerSize',7);
0193      set(h(3),'Marker','diamond','MarkerSize',3);

Generated on Wed 21-Jun-2017 09:29:07 by m2html © 2005