Pore3dCouple

Step 1Step 2Step 3

%step1: packing the elements
clear;
fs.randSeed(1);%random model seed, 1,2,3...
B=obj_Box();%declare a box object
B.name='pore3dCoupleTest';
%--------------initial model------------
B.GPUstatus='auto';%program will test the CPU and GPU speed, and choose the quicker one
B.ballR=0.02/2;
B.sampleW=0.2;%width, length, height, average radius
B.sampleL=0.2;
B.sampleH=0.22;
%B.setType('TriaxialCompression');
B.setType('Fluid');%
B.buildInitialModel();%B.show();
d=B.d;
%you may change the size distribution of elements here, e.g. d.mo.aR=d.aR*0.95;

%---------- gravity sedimentation
B.gravitySediment(1);%you may use B.gravitySediment(10); to increase sedimentation time (10)
%B.compactSample(2);%input is compaction time
mfs.reduceGravity(B.d,10);
d.balance('Standard',2);

%------------return and save result--------------
d.status.dispEnergy();%display the energy of the model
d.show('-aR');
d.mo.setGPU('off');
d.clearData(1);%clear dependent data
d.recordCalHour('Step1Finish');
save(['TempModel/' B.name '1.mat'],'B','d');
save(['TempModel/' B.name '1R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();
%d.show('-Id');

%please run user_Pore3dTest1.m first
%set the material of the model
clear
load('TempModel/pore3dCoupleTest1.mat');
B.setUIoutput();%set output of message
d=B.d;
d.calculateData();
d.mo.setGPU('off'); 
d.getModel();%get xyz from d.mo
d.resetStatus();

%d.show('aR');return

%----------set material of model
matTxt=load('Mats\RockHydro.txt');
Mats{1,1}=material('RockHydro',matTxt,B.ballR);
Mats{1,1}.Id=1;
d.Mats=Mats;
%----------end set material of model
d.showB=2;

%---------assign material to layers and balance the model
d.mo.setShear('on');
%d.removeFixId('Z',[d.GROUP.lefPlaten;d.GROUP.rigPlaten]);
%d.addFixId('Z',[d.GROUP.lefPlaten(end-1),d.GROUP.lefPlaten(end),d.GROUP.rigPlaten(end-1),d.GROUP.rigPlaten(end)]);
d.setGroupMat('sample','RockHydro');
d.groupMat2Model({'sample'});
platenIds=[d.GROUP.lefPlaten;d.GROUP.rigPlaten];
d.mo.aKN(platenIds)=d.mo.aKN(platenIds)*0.1;
d.mo.aMUp(platenIds)=0;
d.mo.mGZ(:)=0;

%---------start computing
d.mo.setGPU('auto');
d.breakGroup();
d.connectGroup('sample');
d.removePrestress(0.1);
d.balance('Standard',2);
d.connectGroup('sample');
d.removePrestress(0.1);
d.balance('Standard',2);

d.show();
%---------end assign material to layers and balance the model1. 

%---------save the data
d.mo.setGPU('off'); 
d.clearData(1);
d.recordCalHour('Step2Finish');
save(['TempModel/' B.name '2.mat'],'B','d');
save(['TempModel/' B.name '2R' num2str(B.ballR) '-distri' num2str(B.distriRate)  'aNum' num2str(d.aNum) '.mat']);
d.calculateData();

warning off;
clear;
load('TempModel\pore3dCoupleTest2.mat');
B.setUIoutput();
d=B.d;
d.calculateData();
d.mo.setGPU('off');
d.getModel();
d.resetStatus();
%-----------------------
d.mo.mVis=d.mo.mVis./d.vRate*0.1;
d.mo.SET.isBalanceData=1;

p=build2pore(d);
p.setInitialPores();
%p.setPlaten('fix')

p.fKFlow(:)=1e-8;
%p.setTimeStep();
p.dT=d.mo.dT;
%tInterval=floor(p.dT/d.mo.dT);
%p.dT=d.mo.dT*tInterval;
%p.SET.tInterval=tInterval;
%------------------------
cXYZ=[median(p.X),median(p.Y),median(p.Z)];
[~,pId]=min(sum(([p.pX,p.pY,p.pZ]-cXYZ).^2,2));
refP=p.Fluids.refP;
cPre=refP+4e8;

p.setPressure(pId,cPre);
p.SET.totalMass0=sum(p.pMass);
pSign=p.pSign;

p.SET.pVolOld=p.pVol;
p.SET.killPoreId=p.SET.pVolOld<10*median(p.SET.pVolOld);

p.setCouple('on');
d.showFilter('SlideY',0.48,1);
% h=line(0,NaN);
totalCircle=20;stepNum=50*4;
d.tic(totalCircle);
save(['data\step\loopNum',num2str(0),'.mat'],'d');
for ii=1:totalCircle
    assert(~any(isnan(p.pPressure)),'boomshakalaka!!!')
    p.SET.pPressure=p.pPressure-refP;
     p.show('SETpPressure');
    view([0,-1,0]);colorbar();
    drawnow;
    for jj=1:stepNum
    p.setPressure(pId,cPre);
    %disp(sum(p.pMass)-p.SET.totalMass0);
    %disp(max(max(sqrt(d.mo.nFnX.^2+d.mo.nFnY.^2+d.mo.nFnZ.^2)./d.mo.nKNe))./median(d.aR));
    p.balance(2);
    if p.isCouple
        p.fluid2Ball();
        p.d.balance();%p.ball2Fluid();
    end
%     k=(max(abs(sum(p.poreFlowMass,2)./p.pMass)));
%     [k,ki]=max(abs(log10(p.SET.pVolOld./p.pVol)));
%     ki;
%     h.XData(end+1)=h.XData(end)+1;
%     h.YData(end+1)=k;
%     drawnow limitrate 
    %p.show('SETpPressure');
    %view([0,-1,0]);
    %drawnow
    end    
    d.toc();
    save(['data\step\loopNum',num2str(ii),'.mat'],'d');
end
d.showFilter('SlideY',0.48,1);
figure;p.show('SETpPressure');view([0,-1,0]);
d.show('Displacement')