/[escript]/branches/diaplayground/ripley/src/Brick.cpp
ViewVC logotype

Diff of /branches/diaplayground/ripley/src/Brick.cpp

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 3913 by caltinay, Tue Jun 19 06:21:58 2012 UTC revision 3943 by caltinay, Thu Aug 9 04:43:24 2012 UTC
# Line 32  namespace ripley { Line 32  namespace ripley {
32  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,  Brick::Brick(int n0, int n1, int n2, double x0, double y0, double z0,
33               double x1, double y1, double z1, int d0, int d1, int d2) :               double x1, double y1, double z1, int d0, int d1, int d2) :
34      RipleyDomain(3),      RipleyDomain(3),
     m_gNE0(n0),  
     m_gNE1(n1),  
     m_gNE2(n2),  
35      m_x0(x0),      m_x0(x0),
36      m_y0(y0),      m_y0(y0),
37      m_z0(z0),      m_z0(z0),
38      m_l0(x1-x0),      m_l0(x1-x0),
39      m_l1(y1-y0),      m_l1(y1-y0),
40      m_l2(z1-z0),      m_l2(z1-z0)
     m_NX(d0),  
     m_NY(d1),  
     m_NZ(d2)  
41  {  {
42        // ignore subdivision parameters for serial run
43        if (m_mpiInfo->size == 1) {
44            d0=1;
45            d1=1;
46            d2=1;
47        }
48    
49        bool warn=false;
50        // if number of subdivisions is non-positive, try to subdivide by the same
51        // ratio as the number of elements
52        if (d0<=0 && d1<=0 && d2<=0) {
53            warn=true;
54            d0=(int)(pow(m_mpiInfo->size*(n0+1)*(n0+1)/((float)(n1+1)*(n2+1)), 1./3));
55            d1=(int)(d0*n1/(float)n0);
56            d2=m_mpiInfo->size/(d0*d1);
57            if (d0*d1*d2 != m_mpiInfo->size) {
58                // ratios not the same so leave "smallest" side undivided and try
59                // dividing 2 sides only
60                if (n0>=n1) {
61                    if (n1>=n2) {
62                        d0=d1=0;
63                        d2=1;
64                    } else {
65                        d0=d2=0;
66                        d1=1;
67                    }
68                } else {
69                    if (n0>=n2) {
70                        d0=d1=0;
71                        d2=1;
72                    } else {
73                        d0=1;
74                        d1=d2=0;
75                    }
76                }
77            }
78        }
79        if (d0<=0 && d1<=0) {
80            warn=true;
81            d0=int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n1+1)));
82            d1=m_mpiInfo->size/d0;
83            if (d0*d1*d2 != m_mpiInfo->size) {
84                // ratios not the same so subdivide side with more elements only
85                if (n0>n1) {
86                    d1=1;
87                } else {
88                    d0=1;
89                }
90            }
91        } else if (d0<=0 && d2<=0) {
92            warn=true;
93            d0=int(sqrt(m_mpiInfo->size*(n0+1)/(float)(n2+1)));
94            d2=m_mpiInfo->size/d0;
95            if (d0*d1*d2 != m_mpiInfo->size) {
96                // ratios not the same so subdivide side with more elements only
97                if (n0>n2) {
98                    d2=1;
99                } else {
100                    d0=1;
101                }
102            }
103        } else if (d1<=0 && d2<=0) {
104            warn=true;
105            d1=int(sqrt(m_mpiInfo->size*(n1+1)/(float)(n2+1)));
106            d2=m_mpiInfo->size/d1;
107            if (d0*d1*d2 != m_mpiInfo->size) {
108                // ratios not the same so subdivide side with more elements only
109                if (n1>n2) {
110                    d2=1;
111                } else {
112                    d1=1;
113                }
114            }
115        }
116        if (d0<=0) {
117            // d1,d2 are preset, determine d0
118            d0=m_mpiInfo->size/(d1*d2);
119        } else if (d1<=0) {
120            // d0,d2 are preset, determine d1
121            d1=m_mpiInfo->size/(d0*d2);
122        } else if (d2<=0) {
123            // d0,d1 are preset, determine d2
124            d2=m_mpiInfo->size/(d0*d1);
125        }
126    
127        m_NX=d0;
128        m_NY=d1;
129        m_NZ=d2;
130    
131      // ensure number of subdivisions is valid and nodes can be distributed      // ensure number of subdivisions is valid and nodes can be distributed
132      // among number of ranks      // among number of ranks
133      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)      if (m_NX*m_NY*m_NZ != m_mpiInfo->size)
134          throw RipleyException("Invalid number of spatial subdivisions");          throw RipleyException("Invalid number of spatial subdivisions");
135    
136      if ((n0+1)%m_NX > 0 || (n1+1)%m_NY > 0 || (n2+1)%m_NZ > 0)      if (warn) {
137          throw RipleyException("Number of elements+1 must be separable into number of ranks in each dimension");          cout << "Warning: Automatic domain subdivision (d0=" << d0 << ", d1="
138                << d1 << ", d2=" << d2 << "). This may not be optimal!" << endl;
139        }
140    
141        if ((n0+1)%m_NX > 0) {
142            double Dx=m_l0/n0;
143            n0=(int)round((float)(n0+1)/d0+0.5)*d0-1;
144            m_l0=Dx*n0;
145            cout << "Warning: Adjusted number of elements and length. N0="
146                << n0 << ", l0=" << m_l0 << endl;
147        }
148        if ((n1+1)%m_NY > 0) {
149            double Dy=m_l1/n1;
150            n1=(int)round((float)(n1+1)/d1+0.5)*d1-1;
151            m_l1=Dy*n1;
152            cout << "Warning: Adjusted number of elements and length. N1="
153                << n1 << ", l1=" << m_l1 << endl;
154        }
155        if ((n2+1)%m_NZ > 0) {
156            double Dz=m_l2/n2;
157            n2=(int)round((float)(n2+1)/d2+0.5)*d2-1;
158            m_l2=Dz*n2;
159            cout << "Warning: Adjusted number of elements and length. N2="
160                << n2 << ", l2=" << m_l2 << endl;
161        }
162    
163        m_gNE0=n0;
164        m_gNE1=n1;
165        m_gNE2=n2;
166    
167      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2) || (m_NZ > 1 && (n2+1)/m_NZ<2))      if ((m_NX > 1 && (n0+1)/m_NX<2) || (m_NY > 1 && (n1+1)/m_NY<2) || (m_NZ > 1 && (n2+1)/m_NZ<2))
168          throw RipleyException("Too few elements for the number of ranks");          throw RipleyException("Too few elements for the number of ranks");

Legend:
Removed from v.3913  
changed lines
  Added in v.3943

  ViewVC Help
Powered by ViewVC 1.1.26