% FEMLAB Model M-file
% Generated by FEMLAB 3.1 (FEMLAB 3.1.0.157, $Date: 2004/11/12 07:39:54 $)

flclear fem

% % COMSOL version
% clear vrsn
% vrsn.name = 'COMSOL 3.2';
% vrsn.ext = '';
% vrsn.major = 0;
% vrsn.build = 222;
% vrsn.rcs = '$Name:  $';
% vrsn.date = '$Date: 2005/09/01 18:02:30 $';
% fem.version = vrsn;

% Constants
fem.const.Kappa=1e-6;   %heat diffusivity
fem.const.rho=3300;     %reference density
fem.const.Cp=1250;      %Heat capacity
fem.const.k=fem.const.Kappa*fem.const.rho*fem.const.Cp;
                        %Heat conductivity
fem.const.Rp=70e3;      %Radius of plume
fem.const.g=9.81;       %Acceleration of gravity
fem.const.eta=1e20;     %Mantle viscosity
fem.const.alpha=3e-5;   %Thermal expansion coefficient
fem.const.Tmax=100;     %maximum temperature excess in plume
fem.const.yr=365.24*24*3600; %secons per year
fem.const.tmax=20e6*fem.const.yr;%duration of model
D=700e3;                %depth of model

% Geometry: two squares next to each other
% Only one is really needed but this makes the definition of the square
% mesh easier
g1=rect2(D,D,'base','corner','pos',[0,-D]);
g2=rect2(D,D,'base','corner','pos',[-D,-D]);

% Geometry
clear s
s.objs={g1, g2};
s.name={'R1','R2'};
s.tags={'g1','g2'};

fem.draw=struct('s',s);
fem.geom=geomcsg(fem);

%% Initialize mesh
% pt1 is a vector that specified the location of square mesh node point
% There are several possibilities here, from coarser to finer
% un-comment next for low resolution
% ptl=[0,2,4,6,8,10,12,15,18,22,26,30,36,42,50,60,70,80,90,100];
% un-comment next for medium resolution
% ptl=[0,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,18,20,...
%     22,24,26,28,30,33,36,39,42,46,50,55,60,...
%     65,70,75,80,85,90,95,100];
% un-comment next for high resolution
ptl=[0,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,...
       32,   34,   36,   38,   40,...
       42,   44,   46,   48,   50,...
       52,      55,      58,...
    61,      64,      67,      70,...
          73,      76,      79,...
       82,      85,      88,...
    91,      96,              100];

% In some boundaries, the finer mesh is at high boundary index values
ptl2=100-fliplr(ptl);
% Making mapped quad mesh, for this particular strategy
fem.mesh=meshmap(fem, ...
    'edgelem',{1,ptl2,2,ptl2,3,ptl2,4,ptl2,6,ptl,5,ptl,7,ptl2});
% If you prefer triangles, use this instead
% fem.mesh=meshinit(fem,...
%     'hmaxsub',[1;0.2],...
%     'Hgrad',1.07,...
%     'hmaxedg',[4;0.02]);
% Visualize the mesh
figure(2);
meshplot(fem); drawnow


%% Define applications
% (Default values are not included)
% Application mode 1: Navier Stokes
clear appl
appl.mode.class = 'FlNavierStokes';
appl.gporder = {4,2};
appl.cporder = {2,1};
appl.assignsuffix = '_ns';
clear bnd
% section 3
bnd.u0 = {0,0,0};
bnd.type = {'uv','out','neutral'};
bnd.ind = [2,2,1,3,2,1,2];
appl.bnd = bnd;
clear equ
equ.gporder = {{1;1;2}};
equ.cporder = {{1;1;2}};
% Buoyancy force, must be linked to temperature excess
equ.F_y = 'T*alpha*rho*g';
equ.eta='eta';
equ.ind = [1, 1];
appl.equ = equ;
fem.appl{1} = appl;

% Application mode 2: Convection/Conduction
clear appl
appl.mode.class = 'FlConvCond';
appl.assignsuffix = '_cc';
clear prop
prop.weakconstr=struct('value',{'off'},'dim',{{'lm4'}});
appl.prop = prop;
clear bnd
% T-base is the Gaussian temperature profile at the bottom
bnd.T0 = {0,'Tbase',0,0,0};
bnd.type = {'q0','T','T','qc','cont'};
bnd.ind = [4,2,3,5,2,3,4];
%bnd.ind = [1,2,3,4];
appl.bnd = bnd;
clear equ
equ.k = 'k';
equ.rho = 'rho';
equ.C = 'Cp';
equ.u = 'u';
equ.v = 'v';
equ.ind = [1, 1];
appl.equ = equ;
fem.appl{2} = appl;
fem.sdim = {'x','y'};
fem.border = 1;

%% Boundary expressions 
clear bnd
bnd.ind = [1,2,1,1,2,1,1];
bnd.dim = {'u','v','p','T'};
bnd.expr = {'Tbase',{'','Tmax*exp(-(x/Rp)^2)'}};
fem.bnd = bnd;

%% Multiphysics
fem=multiphysics(fem);

% Extend mesh
fem.xmesh=meshextend(fem);

%% Solve problem
% Choose a series of time index for which the model is output. Here 100
% regularly-spaced times frames within the choosed time interval. The
% computation time step is selected automatically for default
nt=41;
t_all=linspace(0,fem.const.tmax,nt);
fem.sol=femtime(fem, ...
    'solcomp',{'u','T','p','v'}, ...
    'outcomp',{'u','T','p','v'}, ...
    'tlist',t_all, ...
    'tout','tlist');

% Save current fem structure for restart purposes
fem0=fem;

%break

it=nt;
%% make movie
% Every time frame, background color is temperature, arrows for
% displacement, and time frame in Ma as plot title. Movie saves as .avi
% uncomment relevent lines (end in %movie)

M1=avifile('plume.avi'); %movie
for it=1:nt %movie
    figure(1); clf;
    hold on;
    axis equal;
    postplot(fem, ...
        'tridata',{'max(T,0)','cont','internal'}, ...
        'trimap','hot(1024)', ...
        'tridlim',[0,100],...
        'arrowdata',{'u','v'}, ...
        'arrowxspacing',19, ...
        'arrowyspacing',19, ...
        'arrowtype','arrow', ...
        'arrowstyle','proportional', ...
        'arrowscale',2,...
        'arrowcolor',[0.2,0.1,0.0],...
        'solnum',it, ...
        'title',sprintf('time: %g Ma',t_all(it)/(1e6*fem.const.yr)), ...
        'refine',5, ...
        'axis',[-0.1,1.1,-1.1,0.1,-1,1]);
    set(1,'position',[50,50,700,300]);
    set(gca,'box','on');axis equal; axis tight;
    axis([-1,1,-1,0]*D)
    drawnow; %movie
    M1=addframe(M1,getframe(1)); %movie
end %movie
M1=close(M1); %movie