Commit ef182913 authored by DeAn Wei's avatar DeAn Wei
Browse files

1. add openmp for HandleExtendedDislocation_MD.

2. add SpecifyEquation for specifing function in Tecplot file.
parent 6e9fd4b5
Loading
Loading
Loading
Loading
+1 −1
Original line number Diff line number Diff line
@@ -13,7 +13,7 @@ Set(CMAKE_CXX_COMPILER "g++")
set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} ${OpenMP_C_FLAGS}")  
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")  

#FIND_PACKAGE( OpenMP REQUIRED)  
FIND_PACKAGE( OpenMP REQUIRED)  
if(OPENMP_FOUND)  
    set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${OpenMP_EXE_LINKER_FLAGS}")  
    message(STATUS "OpenMP_C_FLAGS = " ${OpenMP_C_FLAGS})  
+21 −1
Original line number Diff line number Diff line
@@ -15,6 +15,25 @@
#include <vector>
#include <list>

#ifdef _OPENMP

#include <omp.h>

#define INIT_LOCK(a)    omp_init_lock((a))
#define LOCK(a)         omp_set_lock((a))
#define UNLOCK(a)       omp_unset_lock((a))
#define DESTROY_LOCK(a) omp_destroy_lock((a))

#else /* _OPENMP not defined */

#define INIT_LOCK(a)
#define LOCK(a)
#define UNLOCK(a)
#define DESTROY_LOCK(a)

#endif /* end ifdef _OPENMP */


#define real8 double
using namespace std;

@@ -32,7 +51,7 @@ typedef struct {
typedef struct {
        double  burgMag;
        int     seed;
        int     type;
        int     type, nThreads;
        bool    help;

        vector<string>  inpFiles, outFiles;
@@ -67,6 +86,7 @@ typedef enum{
    OPT_TYPE,
    OPT_PRIVATEVALS,
    OPT_AUXFILE,
    OPT_THREADS,
	OPT_MAX
}OPT_t;

+1 −1
Original line number Diff line number Diff line
@@ -28,8 +28,8 @@ typedef struct{

typedef struct {
    int     timestep;
    real8   box[6];

    vector<vector<real8> > box;
    vector<string> bounds;
    vector<string> variables;
    vector<Atom_t> atom;
+6 −2
Original line number Diff line number Diff line
Set (SRCS 
    Main.cpp  ReadData.cpp  Util.cpp Math.cpp HandleExtendedDislocation.cpp
    WriteData.cpp Parse.cpp)
    WriteData.cpp Parse.cpp Custom.cpp)

set(CMAKE_C_FLAGS "${CMAKE_C_FLAGS} -O3 -Wno-unknown-pragmas -DPARALLEL=1 " )
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -O3 " )
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ")

if(OPENMP_FOUND)  
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -openmp")
endif(OPENMP_FOUND)

add_executable(prodata ${SRCS})
add_definitions( -DPARALLEL=1)

+89 −78
Original line number Diff line number Diff line
@@ -127,9 +127,14 @@ void HandleExtendedDislocation_DDD(InArgs_t *inArgs)
        vector<string>().swap(table.variables);

        if(!ReadTecplotNormalData(inArgs->inpFiles[file], table, secLine))continue;

        vector<string>().swap(words);
        words = split(secLine, "\"");

        if(words.size() == 0){
            Fatal("can not find the timesolution in file %s", inArgs->inpFiles[file].c_str());
        }
 
        if(table.data.size() < 2)Fatal("the data size is %d", table.data.size());
 
        colX = GetColIDFromTable(table, "X");
@@ -385,15 +390,15 @@ typedef struct {
}Probe_t;

typedef struct {
    int timestep;
    long int timestep;
    double x,y,z;
    double separation;
    vector<double> vars;
}EDState_t;

bool com(Atom_t p, Atom_t q){
    if(fabs(p.y - q.y) < 0.1){
        if(fabs(p.z - q.z)< 0.1){
    if(fabs(p.y - q.y) < 1.0E-20){
        if(fabs(p.z - q.z)< 1.0E-20){
            return(p.x<q.x);
        }else{
            return(p.z < q.z);
@@ -414,21 +419,22 @@ bool com3(EDState_t p, EDState_t q)

void HandleExtendedDislocation_MD(InArgs_t *inArgs)
{
    int         i, j, file, index, lastIndex, firstNbr, indexVar;
    int         i, j, file, index, firstNbr, indexVar;
    int         nums[2] = {18, 24}, logfile;
    int         stepID;
    double      cutofflen = 2.556, alpha = 0.005, beta = 1.0;
    double      position[3], separation, dis, maxDis;
    double      cutofflen = 2.556, alpha = 0.1, beta = 1.0;
    double      position[3];
    MgData_t    mg;
    LineList_t  list;
    double      dir[3] = {0, 1, 0}, p0[3], dval=5; 
    double      dir[3] = {0, 1, 0}, p0[3], dval=5, effeSepRange[2] = {1.0,40.0}; 
    bool        noP0 = 1;
    string      paraName("para"), dvarName("dvar"), dvar("c_vcna"), p0Name("p0"), dirName("dir");
    string      numsName("nums");
    string      numsName("nums"), ppstring("pp"), separationName("separation");
    int     lastIndex = 0;
    double  dis, minDis;

    vector<Atom_t>::iterator     it; 
    vector<EDState_t>   states;
    EDState_t   state;

    ofstream out;
    out.open(inArgs->outFiles[0].c_str(), ios::out);
@@ -477,76 +483,88 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
            nums[1] = atoi(inArgs->priVars[index].vals[1].c_str());
        }
    }
    printf("The range of effective atoms of a paritial (nums) is [%d,%d]\n", nums[0], nums[1]);
    printf("The range of effective number of a paritial (nums) is [%d,%d]\n", nums[0], nums[1]);

    if((index = GetValID(inArgs->priVars, separationName)) < inArgs->priVars.size()){
        if(inArgs->priVars[index].vals.size() == 2){
            effeSepRange[0] = atof(inArgs->priVars[index].vals[0].c_str());
            effeSepRange[1] = atof(inArgs->priVars[index].vals[1].c_str());
        }
    }
    printf("The range of effective separtion of the extended dislocation (separation) is [%f,%f]\n", effeSepRange[0], effeSepRange[1]);

    if(inArgs->help)return;
    logfile = ReadDataFromMDLogFile(inArgs->auxFiles, list);

    states.resize(inArgs->inpFiles.size());
    for(file=0; file<inArgs->inpFiles.size(); file++){
        states[file].timestep = 1E10;
    }

#pragma omp parallel for private(mg, it, i, j, indexVar, firstNbr, lastIndex, dis, minDis) shared(states) 
    for(file=0; file<inArgs->inpFiles.size(); file++){
        #pragma omp critical
        {
            ReadMGDataFile(inArgs->inpFiles[file], mg);
        }

        for(i=0; i<mg.variables.size(); i++){
            if(dvar == mg.variables[i])break;
        }
        if((i=indexVar) == mg.variables.size()){
        if((indexVar = i) == mg.variables.size()){
            Fatal("There is no %s in the mg file %s", dvar.c_str(), 
                    inArgs->inpFiles[file].c_str());
        }
        
        for(it = mg.atom.end(), i=mg.atom.size()-1; it != mg.atom.begin(); it--, i--){
            if((*it).x < mg.box[0] + (mg.box[1]-mg.box[0])*0.05 ||
               (*it).y < mg.box[2] + (mg.box[3]-mg.box[2])*0.05 ||
               (*it).z < mg.box[4] + (mg.box[5]-mg.box[4])*0.05 ||
               (*it).x > mg.box[0] + (mg.box[1]-mg.box[0])*0.95 ||
               (*it).y > mg.box[2] + (mg.box[3]-mg.box[2])*0.95 ||
               (*it).z > mg.box[4] + (mg.box[5]-mg.box[4])*0.95){
        for(i=0; i<mg.atom.size(); ){
            if(mg.atom[i].vars[indexVar] != dval){
                mg.atom.erase(mg.atom.begin()+i);
                continue;
            }

                if(mg.atom[i].vars[indexVar] == dval){
                    mg.atom.erase(it); 
            if(mg.atom[i].x < mg.box[0][0] + 2 ||
               mg.atom[i].y < mg.box[1][0] + 2 ||
               mg.atom[i].z < mg.box[2][0] + 2 ||
               mg.atom[i].x > mg.box[0][1] - 2 ||
               mg.atom[i].y > mg.box[1][1] - 2 ||
               mg.atom[i].z > mg.box[2][1] - 2){
                mg.atom.erase(mg.atom.begin()+i);
                continue;
            }

            i++;
        }

        if(mg.atom.size() == 0){
            continue;
        }
        
        if(noP0 ==1 && file == 0){
            j = 0;
            for(i=0; i<mg.atom.size(); i++){
                 if(mg.atom[i].vars[indexVar] != dval)continue;
                 j++;
                p0[0] += mg.atom[i].x;
                p0[2] += mg.atom[i].z;
            }
            if(j==0){
                Fatal("The");
            }
            p0[1]  = mg.box[2];
            p0[0] /= ((double)j);
            p0[2] /= ((double)j);
            printf("The initial probe point (p0) is %f %f %f\n", p0[0], p0[1], p0[2]);
            Fatal("The initial probe point (p0) is %f %f %f\n", p0[0], p0[1], p0[2]);
        }

        sort(mg.atom.begin(), mg.atom.end(), com);

        vector<Probe_t> probes;
        Probe_t         probe;
        double          d = 0.0;

        probe.y = p0[1];
        probe.y = mg.atom[0].y;
        probe.z = p0[2];

        lastIndex = 0;  
        while(1){
            vector<Atom_t *>().swap(probe.nbr);

            probe.x += (dir[0]*d*cutofflen*alpha);
            probe.y += (dir[1]*d*cutofflen*alpha);
            probe.z += (dir[2]*d*cutofflen*alpha);
            
            firstNbr = 1;
            maxDis = -1E10;
            minDis = 1E10;
            for(i=lastIndex; i<mg.atom.size(); i++){
                dis = sqrt(pow((mg.atom[i].y-probe.y), 2) +
                      pow((mg.atom[i].z-probe.z), 2) );
                if(dis < cutofflen && mg.atom[i].vars[indexVar] == dval){
                if(dis < minDis){
                    minDis = dis;
                }

                if(dis < cutofflen){
                    if(firstNbr){
                        lastIndex = i;
                        firstNbr = 0;
@@ -558,24 +576,19 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
                probes.push_back(probe);
                vector<Atom_t *>().swap(probe.nbr);
            }
            if(probe.y > mg.atom.back().y ||
               probe.z > mg.atom.back().z)break;
            d++;
        }

#if 0
        printf("Timestep %d, Atoms %d, Bouds %s %s %s, ", 
               mg.timestep, (int)mg.atom.size(), mg.bounds[0].c_str(), 
               mg.bounds[1].c_str(), mg.bounds[2].c_str());
//        printf("Box %e %e %e %e %e %e\n", mg.box[0], mg.box[1],
//               mg.box[2], mg.box[3], mg.box[4], mg.box[5]);
            if(probe.y > mg.atom.back().y)break;

        printf("variables are ");
        for(i=0; i<mg.variables.size(); i++){
            printf("%s(%d) ", mg.variables[i].c_str(), (int)mg.atom.size());
            if(minDis < alpha){
                probe.x += (dir[0]*alpha);
                probe.y += (dir[1]*alpha);
                probe.z += (dir[2]*alpha);
            }else{
                probe.x += (dir[0]*minDis);
                probe.y += (dir[1]*minDis);
                probe.z += (dir[2]*minDis);
            }
        }
        printf("\n");
#endif

        for(i=0; i<probes.size(); i++){
            probes[i].x = 0.0;
            probes[i].y = 0.0;
@@ -592,35 +605,31 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
        }
        sort(probes.begin(), probes.end(), com2);

        separation = 0;
        double separation = 0;

        if(probes.size() > 0){
            separation = fabs(probes[0].y-probes.back().y);
            if(separation > 0.1 && separation < 20){
                state.x = (probes[0].x+probes.back().x)*0.5;
                state.y = (probes[0].y+probes.back().y)*0.5;
                state.z = (probes[0].z+probes.back().z)*0.5;
                state.timestep = mg.timestep;
                state.separation = separation;
            if(separation >  effeSepRange[0] && separation <  effeSepRange[1]){
                states[file].x = (probes[0].x+probes.back().x)*0.5;
                states[file].y = (probes[0].y+probes.back().y)*0.5;
                states[file].z = (probes[0].z+probes.back().z)*0.5;
                states[file].timestep = mg.timestep;
                states[file].separation = separation;

                if(logfile ==1){
                    state.vars.resize(list.variables.size()-1);
                    states[file].vars.resize(list.variables.size()-1);
                    for(i=0; i<list.data[0].size(); i++){
                        if(state.timestep == ((int)list.data[0][i]))break;
                        if(states[file].timestep == ((int)list.data[0][i]))break;
                    }
                    for(j=0; j<state.vars.size(); j++){
                        state.vars[j] = list.data[j+1][i];
                    for(j=0; j<states[file].vars.size(); j++){
                        states[file].vars[j] = list.data[j+1][i];
                    }
                }
                states.push_back(state);
            }
        }
        
    }

    if(states.size() == 0){
        Fatal("No out data can be dumped, please check parameters");
    }

    sort(states.begin(), states.end(), com3);

    out << "variables = timestep, separation, y";
@@ -632,6 +641,8 @@ void HandleExtendedDislocation_MD(InArgs_t *inArgs)
    out << endl;

    for(i=0; i<states.size(); i++){
        if(states[i].timestep > 1E9)break;
        
        out << states[i].timestep << " ";
        out << states[i].separation << " ";
        out << states[i].y  << " ";
Loading