testMimo

Test of function mimo (multiple input, multiple output).

Contents

Syntax

testMimo

Description

testMimo compares by mimo.html computed scattered fields (at receivers positions and on region of interest ) with reference solutions (referenceData2D.html and referenceData3D.html). This comparison is done for several grid sizes.

This test should run in environment runtests, see runtests.html via

runtests('m')

Input Arguments

Further parameters are desribed in the code in the section "General settings", "Discretization parameters", and "Experimental set-up".

Output Arguments

Figures, see runtests.html.

More About

See Also

Code

Run this file standalone... (not recommended)

If you call testMimo 'standalone' you can use the following inside if

if 0
    init;
    seti = struct; % old: setInput;
    seti.dim = 2;
    seti.model = 'helmholtz2D';
end

General settings

seti.measType = 'farField'; % (can be choosen in 2D and 3D)
%seti.measType = 'nearField'; % (can be choosen only in 2D!)

useReferenceData3D = 1;

HMode = 0;      % if 1: using helmholtzHMode instead of helmholtz (not available in public version).
useCoarse = 0;  % if 1: using coarse grid yes or no (logical) (not available in public version and not recommended).

Setting: check consistency

if ~exist('seti', 'var') || ~isfield(seti,'dim')
    seti.dim = 2;
    if HMode
        seti.model = ['helmholtzHMode' int2str(seti.dim) 'D'];
    else
        seti.model = ['helmholtz' int2str(seti.dim) 'D'];
    end
end
seti.fileSuffix = sprintf('_testMimo%dD',seti.dim);

Discretization parameters

% Note: Mexp = -1 -> M = floor(1/seti.nCD) = 0 -> twogrid is not used
if seti.dim == 3
    longtest = 1;
    if longtest
        NMat = sort([2.^(6:8) 2.^(6:8)-1]);
    else
        NMat = sort([2.^(6:6) 2.^(6:6)-1]);
    end
elseif seti.dim == 2
    if HMode % HMode not supported in public version
        NMat = sort([2.^(7:10) 2.^(7:10)-1]);
    else
        NMat = sort([2.^(5:9) 2.^(5:9)-1]);
    end
end

if useCoarse % useCoarse is not supported in public version.
    Mexps = [-1, 5/6, 2/3, 1/2];
else
    Mexps = 0; % Something with length 1.
    % mCD = 0 will be set later, so the value Mexps will not be used
    % but the length is used in a for-loop.
end

Experimental set-up

Changeable parameters

Fixed parameters (for this test)

seti.k = 2;

% helmholtz*D
seti.qBall = 0.8;
seti.rBall = 0.015;
seti.rCD = 0.2;

% helmholtzHMode*D (not available in public version)
%seti.qBall = 3.0; % contrast of ball
%seti.rBall = 0.5; % radius of ball
%seti.rCD = 2.5;

seti.measNb = 2;
seti.radSrc = 4;
seti.radMeas = 4;

% Fixed parameters:
seti.contrast = ['referenceBall' int2str(seti.dim) 'D'];
seti.incType = 'planeWave';
seti.incNb = 1;

setGeomSim

disp('# setGeomSim')
seti = setGeomSim(seti); % set geometry and simulation
%seti = setGeomSim(seti,1); % to debug with output...

testing mimo

fprintf('testing mimo for dim = %i, model = %s\n', seti.dim, seti.model);
fprintf('\n        |         | rel. error  | rel. error  | grid. gen.  |        ');
fprintf('\nPoints  | Points  | of u_s      | of u_s      | +ref.u gen. | comp.  ');
fprintf('\nin CD(N)| in CD(M)| on ROI      | at RX       | time        | time.  ');
fprintf('\n------- + ------- + ----------- + ----------- + ----------- + --------\n');

for j = 1:length(NMat)
    if j > 1 && length(Mexps) > 1 fprintf('------- + ------- + ----------- + ----------- + ----------- + --------\n');end%#ok<SEPEX>

    tic;
    seti.nCD = NMat(j); % number of collocation points in 1 dimension of CD (computational domain)

    seti = setGeomSim(seti);

    if strcmp(seti.model, 'helmholtz2D')
        [uScattRXRef, uScattROIRef] = referenceData2D(seti);
    elseif strcmp(seti.model, 'helmholtzHMode2D')
        [uScattRXRef, uScattROIRef] = referenceDataHMode2D(seti);
    elseif strcmp(seti.model, 'helmholtz3D')
        % If you want to skip the reference solution
        % and compare with the first computation for this N,
        % then change useReferenceData3D to false (0).
        % Note that the reference solution is expensive.
        if useReferenceData3D
            [uScattRXRef, uScattROIRef] = referenceData3D(seti);
        else
            uScattRXRef = 0;
            uScattROIRef = 0;
        end
    else
        disp('testMimo: Unknown Model!');
        return;
    end
    tref = toc;

    for k = 1:length(Mexps)
        tic;

        if useCoarse
            seti.mCD = floor(seti.nCD^Mexps(k));
            [~,seti] = evalc('setCoarse(seti);');
        else
            seti.mCD = 0;
        end

        E = @(x) x;
        N = seti.nROI;
        %E = @(x) extendROItoCD(x,seti.ROImask);
        %N = seti.nCD;

        tgrid = toc;
        tic;

        [FFqMeas,~,FFqROI] = mimo(seti, seti.qROIexact, 'simo');
        uScattRX = FFqMeas;
        if seti.incNb == 1
            uScattROI = FFqROI(:,1);
        else
            error('testMimo does not support more than one transmitter yet.')
        end

        t = toc;

        % If reference solution is not set use this one
        if uScattROIRef == 0
            uScattRXRef = uScattRX;
            uScattROIRef = uScattROI;
        end

        errorROI = norm( seti.ballMask(seti.ROImask).*(uScattROI-uScattROIRef) )...
            /norm( seti.ballMask(seti.ROImask).*uScattROIRef );
        errorRX = norm( uScattRX-uScattRXRef)/norm(uScattRXRef);

        fprintf(' %4d^%i |  %4d^%i |  %10f |  %10f |  %10f | %f\n',seti.nCD, seti.dim, seti.mCD, seti.dim, errorROI, errorRX, tgrid+tref, t);

        part = @real;
        s = N*ones(1,seti.dim);

        EuScattROI = reshape(part(E(uScattROI)), s);
        EuScattROIRef = reshape(part(E(uScattROIRef)), s);
        Ediff = reshape(part(E(uScattROI-uScattROIRef)), s);

        switch seti.dim
            case 2
                sliceWarn = '';
            case 3
                % Multiple slices
                v = round(linspace(1, N, min(10,N)));

                figure(28); slice(EuScattROI,v,N,1); title('Computed scattered field'); colormap(seti.colormap); colorbar; shading flat;
                figure(29); slice(EuScattROIRef,v,N,1); title('Reference scattered field'); colormap(seti.colormap); colorbar; shading flat;
                figure(30); slice(Ediff,v,N,1); title('Difference'); colormap(seti.colormap); colorbar; shading flat;

                % Display a slice of the data.
                sliceWarn = ' (slice z = 0)';
                N2 = ceil((N-1)/2);
                EuScattROI = EuScattROI(:,:,N2);
                EuScattROIRef = EuScattROIRef(:,:,N2);
                Ediff = Ediff(:,:,N2);
        end

        if seti.dim == 2
            figNo1 = 21; figNo2 = 22; figNo3 = 23;
        elseif seti.dim == 3
            figNo1 = 25; figNo2 = 26; figNo3 = 27;
        end
        figure(figNo1); imagesc(EuScattROI); axis xy; title(['Computed scattered field' sliceWarn]);colormap(seti.colormap);colorbar;
        figure(figNo2); imagesc(EuScattROIRef); axis xy; title(['Reference scattered field' sliceWarn]);colormap(seti.colormap);colorbar;
        figure(figNo3); imagesc(Ediff); axis xy; title(['Difference' sliceWarn]);colormap(seti.colormap);colorbar;
        drawnow;
        savePngFig(figNo1,0,seti);
        savePngFig(figNo2,0,seti);
        savePngFig(figNo3,0,seti);


    end
end
fprintf('\n');