model = mphopen('1113.mph'); %打开文件

param_name = ["d1" "a3" "a2" "a1"]; % 参数名称
param_unit = ["um" "um" "um" "um"]; % 参数单位
param_Value = [0.1 0.3 0.5;...  %参数值
                50 75 100;...
                1 3 5;...
                0.1 1 2];

% For循环运行所有参数组
para_num=size(param_Value,2);
T_max=zeros(para_num^4,1);
Stress1=zeros(para_num^4,1);
Stress2=zeros(para_num^4,1);
para=ones(para_num^4,4);
i=1;
for para1=1:para_num % 循环研究参数名
    for para2=1:para_num
        for  para3=1:para_num
            for para4=1:para_num
                %记录输出
                tic;
                fprintf('Start simulation %d: %s=%.1f | %s=%.1f | %s=%.1f | %s=%.1f \n',i,...
                    param_name(1),param_Value(1,para1),...
                    param_name(2),param_Value(2,para2),...
                    param_name(3),param_Value(3,para3),...
                    param_name(4),param_Value(4,para4));
                %设定参数
                model.param.set(param_name(1),num2str(param_Value(1,para1)),param_unit(1));
                model.param.set(param_name(2),num2str(param_Value(2,para2)),param_unit(2));
                model.param.set(param_name(3),num2str(param_Value(3,para3)),param_unit(3));
                model.param.set(param_name(4),num2str(param_Value(4,para4)),param_unit(4));
                %运行计算
                model.sol('sol1').runAll;
                model.sol('sol2').runAll;
                model.sol('sol3').runAll;
                %记录数据
                para(i,:)=[param_Value(1,para1) param_Value(2,para2) param_Value(3,para3) param_Value(4,para4)];
                T_max(i) = mphglobal(model,'dom1','dataset','dset1','unit','degC'); %导出指定数据集的探针数据
                Stress1(i) = mphglobal(model,'dom2','dataset','dset3','unit','MPa');
                Stress2(i) = mphglobal(model,'dom6','dataset','dset4','unit','MPa');
                fprintf('Done simulation %d: T_max = %.2f | Stress1 = %.2f | Stress2 = %.2f\n',i,T_max(i),Stress1(i),Stress2(i));
                i=i+1;
                toc
            end
        end
    end
end
save('data.mat','para','T_max','Stress1','Stress2');