Skip to content

Use parfor for parallelism¤

Parallel Computing Toolboxis needed.

1
demo_pushover_parfor_alpha_sy();
Output
Total number of cases: 156 alpha count: 12 sY count: 13 Running serial computation... Serial elapsed time: 12.567896 s Running parallel computation... Starting parallel pool (parpool) using the 'Processes' profile ... Connected to parallel pool with 12 workers. Requested workers: 12 Actual parallel workers used: 12 Parallel pool workers: 12 ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ ============================================================ OpenSeesMatlab v3.8.0.1 OpenSees MEX Interface for MATLAB Copyright (c) 2026, By Yexiang Yan Type 'help OpenSeesMatlab' in MATLAB for documentation. Documentation also available at https://openseesmatlab.readthedocs.io/en/latest/ ============================================================ Parallel elapsed time: 4.485445 s Speedup (serial/parallel): 2.802 CaseID alpha sY maxDisp maxForce success errorMessage ______ _____ __ _______ ________ _______ ____________ 1 0.01 24 2 100.42 true "" 2 0.01 26 2 108.23 true "" 3 0.01 28 2 116.04 true "" 4 0.01 30 2 123.85 true "" 5 0.01 32 2 131.66 true "" 6 0.01 34 2 139.47 true "" 7 0.01 36 2 147.28 true "" 8 0.01 38 2 155.09 true "" 9 0.01 40 2 162.91 true "" 10 0.01 42 2 170.72 true ""
figure_0.png
figure_1.png
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
function [resultsSerial, resultsParallel] = demo_pushover_parfor_alpha_sy()
%DEMO_PUSHOVER_PARFOR_ALPHA_SY
% Parallel pushover study for a 2D three-bar truss model.
%


    %% Parameter sets
    alphaList = 0.01:0.01:0.12;
    sYList    = 24:2:48;


    nAlpha = numel(alphaList);
    nSY    = numel(sYList);
    nCases = nAlpha * nSY;


    fprintf('Total number of cases: %d\n', nCases);
    fprintf('alpha count: %d\n', nAlpha);
    fprintf('sY count: %d\n', nSY);


    %% Generate parameter combinations
    caseTable = zeros(nCases, 3);   % [caseID, alpha, sY]
    k = 0;
    for i = 1:nAlpha
        for j = 1:nSY
            k = k + 1;
            caseTable(k, :) = [k, alphaList(i), sYList(j)];
        end
    end


    %% Preallocate result containers
    resultsSerial   = preallocateResults(nCases);
    resultsParallel = preallocateResults(nCases);


    %% ---------------- Serial run ----------------
    fprintf('\nRunning serial computation...\n');
    tSerial = tic;


    for icase = 1:nCases
        resultsSerial(icase) = runSinglePushoverCase(caseTable(icase, :));
    end


    timeSerial = toc(tSerial);
    fprintf('Serial elapsed time: %.6f s\n', timeSerial);


    %% ---------------- Parallel run ----------------
    if license('test', 'Distrib_Computing_Toolbox')
        fprintf('\nRunning parallel computation...\n');


        pc = parcluster('local');
        nWorkers = 12;
        pc.NumWorkers = nWorkers;
        saveProfile(pc);

        pool = gcp('nocreate');
        if isempty(pool)
            pool = parpool('local', nWorkers);
        elseif pool.NumWorkers ~= nWorkers
            delete(pool);
            pool = parpool('local', nWorkers);
        end


        fprintf('Requested workers: %d\n', 12);
        fprintf('Actual parallel workers used: %d\n', pool.NumWorkers);


        fprintf('Parallel pool workers: %d\n', pool.NumWorkers);


        tParallel = tic;


        parfor icase = 1:nCases
            resultsParallel(icase) = runSinglePushoverCase(caseTable(icase, :));
        end


        timeParallel = toc(tParallel);
        fprintf('Parallel elapsed time: %.6f s\n', timeParallel);


        fprintf('Speedup (serial/parallel): %.3f\n\n\n', timeSerial / timeParallel);
    else
        warning('Parallel Computing Toolbox not found. Parallel run is skipped.');
        timeParallel = NaN;
    end


    %% Show summary table for parallel results if available, otherwise serial
    if all(arrayfun(@(x) ~isempty(x.CaseID), resultsParallel))
        resultsToShow = resultsParallel;
    else
        resultsToShow = resultsSerial;
    end


    T = struct2table(resultsToShow);
    nShow = min(10, height(T));
    disp(T(1:nShow, {'CaseID','alpha','sY','maxDisp','maxForce','success','errorMessage'}))


    %% Plot pushover curves for varying sY at fixed alpha
    alphaFixed = 0.05;
    cmap = parula(nSY);


    figure('Color', 'w', 'Name', sprintf('Pushover Curves - alpha = %.3g', alphaFixed));
    hold on;
    grid on;
    box on;


    for is = 1:nSY
        for icase = 1:nCases
            if resultsToShow(icase).success && ...
               abs(resultsToShow(icase).alpha - alphaFixed) < 1e-12 && ...
               abs(resultsToShow(icase).sY    - sYList(is)) < 1e-12


                plot(resultsToShow(icase).disp, resultsToShow(icase).baseForce, ...
                    'LineWidth', 1.5, ...
                    'Color', cmap(is, :), ...
                    'DisplayName', sprintf('sY = %.3g', resultsToShow(icase).sY));
                break;
            end
        end
    end


    xlabel('Node 4 displacement');
    ylabel('Load factor \times P_x');
    xlim([0 3]);
    title(sprintf('Pushover Curves (fixed \\alpha = %.3g)', alphaFixed));
    legend('Location', 'best');


    %% Plot pushover curves for varying alpha at fixed sY
    sYFixed = 36.0;
    cmap = parula(nAlpha);


    figure('Color', 'w', 'Name', sprintf('Pushover Curves - sY = %.3g', sYFixed));
    hold on;
    grid on;
    box on;


    for ia = 1:nAlpha
        for icase = 1:nCases
            if resultsToShow(icase).success && ...
               abs(resultsToShow(icase).sY    - sYFixed) < 1e-12 && ...
               abs(resultsToShow(icase).alpha - alphaList(ia)) < 1e-12


                plot(resultsToShow(icase).disp, resultsToShow(icase).baseForce, ...
                    'LineWidth', 1.5, ...
                    'Color', cmap(ia, :), ...
                    'DisplayName', sprintf('\\alpha = %.3g', resultsToShow(icase).alpha));
                break;
            end
        end
    end


    xlabel('Node 4 displacement');
    ylabel('Load factor \times P_x');
    xlim([0 3]);
    title(sprintf('Pushover Curves (fixed sY = %.3g)', sYFixed));
    legend('Location', 'best');
end


function results = preallocateResults(nCases)
%PREALLOCATERESULTS Preallocate result struct array.


    results(nCases, 1) = struct( ...
        'CaseID', [], ...
        'alpha', [], ...
        'sY', [], ...
        'disp', [], ...
        'baseForce', [], ...
        'maxDisp', [], ...
        'maxForce', [], ...
        'success', [], ...
        'errorMessage', "" ...
    );
end

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
function out = runSinglePushoverCase(caseRow)
%RUNSINGLEPUSHOVERCASE Run one pushover analysis for one parameter set.


    caseID = caseRow(1);
    alpha  = caseRow(2);
    sY     = caseRow(3);


    out = struct( ...
        'CaseID', caseID, ...
        'alpha', alpha, ...
        'sY', sY, ...
        'disp', [], ...
        'baseForce', [], ...
        'maxDisp', NaN, ...
        'maxForce', NaN, ...
        'success', false, ...
        'errorMessage', "" ...
    );


    try
        %% Create independent OpenSeesMatlab object
        opsMAT = OpenSeesMatlab();
        ops = opsMAT.opensees;


        %% Model parameters
        A      = 4.0;
        E      = 29000.0;
        Nsteps = 1000;
        Px     = 160.0;
        Py     = 0.0;


        %% Build model
        ops.wipe();
        ops.model("basic", "-ndm", 2, "-ndf", 2);


        ops.node(1, 0.0,   0.0);
        ops.node(2, 72.0,  0.0);
        ops.node(3, 168.0, 0.0);
        ops.node(4, 48.0, 144.0);


        ops.fix(1, 1, 1);
        ops.fix(2, 1, 1);
        ops.fix(3, 1, 1);


        % Hardening material
        % H = alpha/(1-alpha) * E
        ops.uniaxialMaterial("Hardening", 1, E, sY, 0.0, alpha / (1 - alpha) * E);


        ops.element("Truss", 1, 1, 4, A, 1);
        ops.element("Truss", 2, 2, 4, A, 1);
        ops.element("Truss", 3, 3, 4, A, 1);


        %% Load pattern
        ops.timeSeries("Linear", 1);
        ops.pattern("Plain", 1, 1);
        ops.load(4, Px, Py);


        %% Analysis settings
        ops.system("ProfileSPD");
        ops.numberer("Plain");
        ops.constraints("Plain");
        ops.integrator("DisplacementControl", 4, 1, 2.0 / Nsteps);
        ops.algorithm("Newton");
        ops.test("NormUnbalance", 1e-8, 10, 0);
        ops.analysis("Static");


        %% Run pushover
        data = zeros(Nsteps + 1, 2);
        data(1, :) = [0.0, 0.0];


        for j = 1:Nsteps
            ok = ops.analyze(1);
            if ok ~= 0
                error("Analysis failed at step %d.", j);
            end


            data(j + 1, 1) = ops.nodeDisp(4, 1);
            data(j + 1, 2) = ops.getLoadFactor(1) * Px;
        end


        %% Store results
        out.disp      = data(:, 1);
        out.baseForce = data(:, 2);
        out.maxDisp   = max(abs(data(:, 1)));
        out.maxForce  = max(abs(data(:, 2)));
        out.success   = true;


    catch ME
        out.success = false;
        out.errorMessage = string(ME.message);
    end
end