Helmer Cluster show case 2: a 4.65 Magnitude earthquake shakes Great Toronto Area

Earthquakes in GTA? Yes, I am kidding. But what if a 4.65 magnitude nuclear explosion occurred at 3km deep underground secret lab, located close to the junction of Kirby road and Kipling Ave, Northwest of Vaughan (see google).

Here is a terrain map of the location. The black star indicates the hypocentre location.

Using SPECFEM3D with P=5, grid interval= 200m, half duration 0.2 second (frequency ~2~3Hz), the computation on my 18 core cluster took 1.5 hours. Here are some snapshots:
T=1.5 second

T=2.5 second

T=4.0 second

Notice the converted shear wave distribution behind the P-wave front? At some locations, the shear wave energy is amplified by the local topography and can severely damage the local surface structures. Good news: the modeled topography is 10 times amplified from the real topography. So the results here have little indication to the real case.

Coming soon:
Parallel hybrid spectral element and finite element simulation of airborne electromagnetic survey over a heterogeneous earth model with multi-scale lossy objects. Stay tuned.

Sage: a nice open source software alternative to Matlab

There are quit a few alternatives to Matlab: GNU Octave, Scilab, freeMat, and Sage. They all have their advantages and disadvantages. But at least there is one suiting your needs. Currently I am looking at Sage, which seems a “confederation” of multiple GNU math libraries. What interested me most is its web based worksheet allowing multiple users to share and collaborate. Here is a screenshot:

The full content of this “test” worksheet is here
If you are using firefox, accept the warning.

Helmer Cluster show case 1: seismic wave propagation in the presence of a fluid filled fracture

It is a show case for both my 18 core cluster and the performance of SPECFEM2D (link) software package. I remember someone did that already using variable grid finite difference. But for irregular fracture, finite element is an optimum solution.

The model is 30 m by 30 m box with a fracture of 1 cm thick and 10 m long. See figure as follows:

It is simply generated by gmsh software package:

// Gmsh project created on Tue Apr 03 22:42:21 2012
lc=0.5;
lf=0.05;
Point(1) = {0, 0, 0, lc};
Point(2) = {0, 30, 0, lc};
Point(3) = {30, 30, 0, lc};
Point(4) = {30, 0, 0, lc};
Line(1) = {2, 1};
Line(2) = {1, 4};
Line(3) = {4, 3};
Line(4) = {3, 2};
Point(5) = {13.5, 15, 0, lf};
Point(6) = {13.5, 5, 0, lf};
Point(7) = {13.51, 15, 0, lf};
Point(8) = {13.51, 5, 0, lf};
Line(5) = {5, 7};
Line(6) = {5, 6};
Line(7) = {6, 8};
Line(8) = {8, 7};
Line Loop(9) = {6, 7, 8, -5};
Line Loop(10) = {1, 2, 3, 4, -6, -7, -8, 5};
Plane Surface(11) = {9};
Plane Surface(10) = {10};
Physical Line("abs") = {4, 1, 2, 3}; # absorbing boundary
Physical Surface("M1") = {10}; # material 1
Physical Surface("M2") = {11}; # material 2

and exported to SPECFEM2D mesh format using a Matlab script:

clear all

filehead='fracture';
fname=[filehead,'.msh'];
Fmesh=[filehead,'_mesh'];
Fmat=[filehead,'_mat'];
Fnodes=[filehead,'_nodesCoords'];
Fsur_abs=[filehead,'_surface_abs'];
Fsur_free=[filehead,'_surface_free'];

%%
% import gmsh file
fid=fopen(fname,'r');

tline=fgetl(fid);
while ~feof(fid)
    if strcmp(tline, '$PhysicalNames')
        tline=fgetl(fid);
        phyNum=str2num(tline);
        phyNames(phyNum,1)=struct('dimension',[], 'number',[], 'name',[]);
        phyNames(phyNum,1).dimension=[];phyNames(phyNum,1).number=[];phyNames(phyNum,1).name=[];
        for i=1:phyNum
            tline=fgetl(fid);
            C = textscan(tline,'%d %d %s');
            phyNames(i,1).dimension=C{1}(1);
            phyNames(i,1).number=C{2}(1);
            phyNames(i,1).name=C{3}{1}(2:end-1);
        end
    elseif strcmp(tline,'$Nodes')
        tline=fgetl(fid);
        nodeNum=str2num(tline);
        Nodes(nodeNum,1)=struct('coord',[]);
        Nodes(nodeNum,1).coord=[];

        for i=1:nodeNum
            tline=fgetl(fid);
            C = textscan(tline,'%d %f %f %f');
            Nodes(i,1).coord=[C{2} C{3} C{4}];
        end
    elseif strcmp(tline,'$Elements')
        tline=fgetl(fid);
        meshNum=str2num(tline);
        elements(meshNum,1)=struct('idx',[],'type',[],'numTags',[],'tags',[],'mesh',[]);
        elements(meshNum,1).idx=[];elements(meshNum,1).type=[];elements(meshNum,1).numTags=[];
        elements(meshNum,1).tags=[];elements(meshNum,1).mesh=[];

        for i=1:meshNum
            tline=fgetl(fid);
            C = textscan(tline,'%d');
            elements(i).idx=C{1}(1);
            elements(i).type=C{1}(2);
            elements(i).numTags=C{1}(3);
            for j=1:elements(i).numTags
                elements(i).tags(j)=C{1}(j+3);
            end
            if elements(i).type==1 % two node line
                elements(i).mesh=[C{1}(end-1:end)]';
            elseif elements(i).type==8 % three nodes line 2nd order
                elements(i).mesh=[C{1}(end-2:end)]';
            elseif elements(i).type==3 % 4-node quadrangle
                elements(i).mesh=[C{1}(end-3:end)]';
            elseif elements(i).type==10 % 9-node second order quadrangle
                elements(i).mesh=[C{1}(end-8:end)]';
            else
                error('Only 1st and 2nd order quadrangle type element are supported');
            end
        end

    end
    tline=fgetl(fid);

end
fclose(fid);

%%
%export nodes to SPECFEM2D
fid=fopen(Fnodes,'w');
fprintf(fid,'%dn',length(Nodes));
for i=1:length(Nodes)
    fprintf(fid,'%12.3f%12.3fn',Nodes(i).coord(1),Nodes(i).coord(2));
end
fclose(fid);

%%
% export mesh to SPECFEM2D
fid=fopen(Fmesh,'w');
fprintf(fid,'%8d%8dn',0,0);
j3=0;j10=0;
for i=1:length(elements)
    if elements(i).type==3
        fprintf(fid,'%8d%8d%8d%8dn',elements(i).mesh);
        j3=j3+1;
    elseif elements(i).type==10
        fprintf(fid,'%8d%8d%8d%8d%8d%8d%8d%8d%8dn',elements(i).mesh);
        j10=j10+1;
    end
end
frewind(fid);
if j3>0
    fprintf(fid,'%8d%8dn',j3,4);
end
if j10>0
    fprintf(fid,'%8d%8dn',j10,9);
end
fclose(fid);

%%
% export material to SPECFEM2D
fid=fopen(Fmat,'w');
for i=1:length(phyNames)
    allNames(i)=phyNames(i).number;
end
for i=1:length(elements)
    if elements(i).type==3 || elements(i).type==10
        mstr=phyNames(allNames==elements(i).tags(1)).name;
        fprintf(fid,'%8dn',str2num(mstr(2:end)));
    end
end
fclose(fid);

%%
% export boundary conditions to SPECFEM2D
for i=1:length(phyNames)
    allNames(i)=phyNames(i).number;
end

meshMatrix=[];j=1;
for i=1:length(elements)
    if elements(i).type==3 || elements(i).type==10
        meshMatrix=[meshMatrix;elements(i).mesh];
    end
end

fabs=fopen(Fsur_abs,'w');fprintf(fabs,'%8d n',0);
ffree=fopen(Fsur_free,'w');fprintf(ffree,'%8d n',0);
ja=0;jf=0;
for i=1:length(elements)
    if elements(i).type==1 || elements(i).type==8   %take edge elements only
        tmp=elements(i).mesh;
        idx=[];
        for j=1:length(meshMatrix(:,1))
            itag=0;
            for k=1:length(tmp)
                if find(meshMatrix(j,:)==tmp(k))
                    itag=itag+1;
                end
            end
            if length(tmp)==itag
                idx=[idx;j];
            end
        end
        if isempty(idx) || length(idx)>1
            error(['No match or Non edge nodes used in line element, please check element ',num2str(i)]);
        end
        mstr=phyNames(allNames==elements(i).tags(1)).name;
        if strcmp(mstr,'free') %unless specified, all boundaries are absorbing
            jf=jf+1;
            if elements(i).type==1
                fprintf(ffree,'%8d%8d%8d%8dn',idx,2,elements(i).mesh);
            elseif elements(i).type==3
                fprintf(ffree,'%8d%8d%8d%8d%8dn',idx,2,elements(i).mesh);
            else
                error('Other type of edge elements are not implemented.');
            end
        else
            ja=ja+1;
            if elements(i).type==1
                fprintf(fabs,'%8d%8d%8d%8dn',idx,2,elements(i).mesh);
            elseif elements(i).type==3
                fprintf(fabs,'%8d%8d%8d%8d%8dn',idx,2,elements(i).mesh);
            else
                error('Other type of edge elements are not implemented.');
            end
        end
    end
end
frewind(ffree);fprintf(ffree,'%8d n',jf);
frewind(fabs);fprintf(fabs,'%8d n',ja);
fclose(ffree);fclose(fabs);

A full simulation 50000 time steps (1d-8 second each step) is completed in 1 hour. Here are a few snapshots:
Time step: 12000

Time step: 17000

Time step: 30000

Next show case of my cluster: Electromagnetic wave diffusion in the presence of irregular orebody using Galerkin Discontinuous high order finite element method. Stick around and know more about what we can do with a home brewed cluster.