From 91af7f365aa202617a8b010596b81b7a04bd7ea9 Mon Sep 17 00:00:00 2001 From: JFL Date: Fri, 7 Jun 2024 09:09:13 +0200 Subject: [PATCH] synch --- include/DFT_Factory.h | 11 ++- include/Lattice.h | 169 +++++++++++++++++++++++++++--------------- include/Species.h | 16 +++- src/Density.cpp | 139 +++++++++++++++++++++------------- src/Species.cpp | 47 +++++++----- 5 files changed, 249 insertions(+), 133 deletions(-) diff --git a/include/DFT_Factory.h b/include/DFT_Factory.h index 794cda3..b299c5e 100644 --- a/include/DFT_Factory.h +++ b/include/DFT_Factory.h @@ -61,7 +61,8 @@ class DFT_Factory options_.addOption("FixedBoundary", &fixed_background_); options_.addOption("FixedBackground", &fixed_background_); options_.addOption("Open_System", &fixed_background_); - + options_.addOption("BoundaryWidth",&boundary_width_); + options_.addOption("MaxIterations", &maxIterations_); options_.addOption("Tolerence", &tol_); @@ -167,6 +168,8 @@ class DFT_Factory theDensity_ = new DensityType(dx_, L_, show_graphics_); #endif + theDensity_->set_boundary_width(boundary_width_); + if(eos_correction_ == LJ_JZG_EOS) eos_ = new LJ_JZG(kT_, rcut1_); // need to add no-shift option // if(eos_correction_ == LJ_MECKE_EOS) eos_ = new LJ_Mecke(kT_, rcut1_); // need to add no-shift option @@ -396,7 +399,9 @@ class DFT_Factory bool verbose_ = true; - string potential_name_ = "LJ"; - string log_file_name_ = "log.dat"; + int boundary_width_ = 0; + + string potential_name_ = "LJ"; + string log_file_name_ = "log.dat"; }; #endif // __LUTSKO_DFT_FACTORY__ diff --git a/include/Lattice.h b/include/Lattice.h index 9ff197e..1833661 100644 --- a/include/Lattice.h +++ b/include/Lattice.h @@ -42,6 +42,9 @@ class Lattice Lattice& operator=(const Lattice &a){copy(a);return *this;} + void set_boundary_width(unsigned b) { boundary_width_ = b;} + unsigned get_boundary_width() { return boundary_width_;} + void copy(const Lattice &ref) { Nx_ = ref.Nx_; @@ -57,9 +60,11 @@ class Lattice L_[0] = ref.L_[0]; L_[1] = ref.L_[1]; - L_[2] = ref.L_[2]; + L_[2] = ref.L_[2]; + + boundary_width_ = ref.boundary_width_; } - + void init(double L[]) { @@ -193,25 +198,30 @@ void init(double L[]) } ////////// Characterizing boundary points: - // We map points on the boundary (ix or iy or iz = 0) - // to an array of NxNy+Nx(Nz-1)+(Ny-1)(Nz-1) points + // Boundary points are those within a distance of b of the edge of the cell. + // We map points on the boundary (ix <= b or ix >= Nx-b; or iy <= b or iy >= Ny-b or iz <= b or iz >= Nz - b ) + // to an array of k(NyNz+ (Nx-k)Nz+(Nx-k)(Ny-k)) points where k = 1+2b + // or, equivalently, k(NxNy + NxNy + NyNz -k(Nx+Ny+Nz) +k*k) points. bool is_boundary_point(long pos) const { int ix, iy, iz; cartesian(pos,ix,iy,iz); - return ((ix % Nx_ == 0) || (iy % Ny_ == 0) || (iz % Nz_ == 0)); + + unsigned b = boundary_width_; + + return ((ix%Nx_ <= b || ix%Nx_ >= Nx_-b) || (iy%Ny_ <= b || iy%Ny_ >= Ny_-b) || (iz%Nz_ <= b || iz%Nz_ >= Nz_-b)); } - long get_Nboundary() const { return Nx_*Ny_+Nx_*Nz_+Ny_*Nz_-Nx_-Ny_-Nz_+1;} + long get_Nboundary() const { int k = 1+2*boundary_width_; return k*(Nx_*Ny_+Nx_*Nz_+Ny_*Nz_-k*(Nx_+Ny_+Nz_)+k*k);} // cartesian coordinates to position long boundary_pos(int v[]) const { return boundary_pos(v[0],v[1],v[2]);} long boundary_pos(int ix, int iy, int iz) const { long pos = -1; - + // First get points in cell while(ix < 0) ix += Nx_; while(iy < 0) iy += Ny_; while(iz < 0) iz += Nz_; @@ -220,18 +230,24 @@ void init(double L[]) while(iy > Ny_-1) iy -= Ny_; while(iz > Nz_-1) iz -= Nz_; - if(iz == 0) - pos = ix*Ny_+iy; - else { - pos = Nx_*Ny_; - if(iy == 0) - pos += ix*(Nz_-1)+iz-1; - else if(ix == 0) { - pos += Nx_*(Nz_-1) + (iy-1)*(Nz_-1)+(iz-1); - } else { - throw std::runtime_error("Non-boundary point given to Lattice::boundary_pos"); - } - } + unsigned b = boundary_width_; + unsigned k = 1+2*b; + + if(iz <= b) + pos = iz*Nx_*Ny_ + ix*Ny_+iy; + else if(iz >= Nz_-b) + pos = (b+1)*Nx_*Ny_ + (iz-(Nz_-b))*Nx_*Ny_ + ix*Ny_+iy; + else if(iy <= b) + pos = k*Nx_*Ny_ + (iy*Nx_ + ix)*(Nz_-k)+iz-1-b; + else if(iy >= Ny_-b) + pos = k*Nx_*Ny_ + (b+1)*Nx_*(Nz_-k) + ((iy-(Ny_-b))*Nx_ + ix)*(Nz_-k)+iz-1-b; + else if(ix <= b) + pos = k*Nx_*Ny_ + k*Nx_*(Nz_-k) + (ix*(Ny_-k)+(iy-b-1))*(Nz_-k)+(iz-b-1); + else if(ix >= Nx_-b) + pos = k*Nx_*Ny_ + k*Nx_*(Nz_-k) + (b+1)*(Ny_-k)*(Nz_-k) + ((ix-(Nx_-b))*(Ny_-k)+(iy-b-1))*(Nz_-k)+(iz-b-1); + else + throw std::runtime_error("Non-boundary point given to Lattice::boundary_pos"); + return pos; } @@ -252,55 +268,84 @@ void init(double L[]) void boundary_cartesian(long pos,int &ix, int &iy, int &iz) const { ix = iy = iz = 0; - - if(pos < Nx_*Ny_) + + unsigned b = boundary_width_; + unsigned k = 1+2*b; + + if(pos < (b+1)*Nx_*Ny_) { - ix = int(pos/Ny_); - iy = pos-Ny_*ix; - } else { - pos -= Nx_*Ny_; - if(pos < Nx_*(Nz_-1)) - { - ix = int(pos/(Nz_-1)); - iz = 1+int(pos-ix*(Nz_-1)); - } else { - pos -= Nx_*(Nz_-1); - if(pos < (Ny_-1)*(Nz_-1)) - { - iy = 1+int(pos/(Nz_-1)); - iz = 1+int(pos - (iy-1)*(Nz_-1)); - } else { - throw std::runtime_error("Non-boundary point given to Lattice::boundary_cartesian"); - } - } - } + //pos = iz*Nx_*Ny_ + ix*Ny_+iy; + iz = pos/(Nx_*Ny_); pos -= iz*Nx_*Ny_; + ix = pos/Ny_; pos -= ix*Ny_; + iy = pos; + } else if(pos < k*Nx_*Ny_) { + //pos = (iz+k-Nz_)*Nx_*Ny_ + ix*Ny_+iy; + iz = Nz_-k+pos/(Nx_*Ny_); pos -= (iz-(Nz_-k))*Nx_*Ny_; + ix = pos/Ny_; pos -= ix*Ny_; + iy = pos; + } else if(pos < k*Nx_*Ny_ + (b+1)*Nx_*(Nz_-k)) { + //pos = k*Nx_*Ny_ + (iy*Nx_ + ix)*(Nz_-k)+iz-1-b; + iy = (pos-k*Nx_*Ny_)/(Nx_*(Nz_-k)); pos -= k*Nx_*Ny_ + iy*Nx_*(Nz_-k); + ix = pos/(Nz_-k); pos -= ix*(Nz_-k); + iz = pos+1+b; + } else if(pos < k*Nx_*Ny_ + k*Nx_*(Nz_-k)) { + //pos = k*Nx_*Ny_ + (b+1)*Nx_*(Nz_-k) + ((iy-(Ny_-b))*Nx_ + ix)*(Nz_-k)+iz-1-b; + iy = Ny_-b + (pos - k*Nx_*Ny_ - (b+1)*Nx_*(Nz_-k))/(Nx_*(Nz_-k)); pos -= k*Nx_*Ny_ + (b+1)*Nx_*(Nz_-k) + (iy-(Ny_-b))*Nx_*(Nz_-k); + ix = pos/(Nz_-k); pos -= ix*(Nz_-k); + iz = pos+1+b; + } else if(pos < k*Nx_*Ny_ + k*Nx_*(Nz_-k) + k*(Ny_-k)*(Nz_-k) ) { + //pos = k*Nx_*Ny_ + k*Nx_*(Nz_-k) + (ix*(Ny_-k)+(iy-b-1))*(Nz_-k)+(iz-b-1); + ix = (pos - k*Nx_*Ny_ - k*Nx_*(Nz_-k))/((Ny_-k)*(Nz_-k)); pos -= k*Nx_*Ny_ + k*Nx_*(Nz_-k) + ix*(Ny_-k)*(Nz_-k); + iy = 1+b+pos/(Nz_-k); pos -= (iy-b-1)*(Nz_-k); + iz = pos+1+b; + } else + throw std::runtime_error("Non-boundary point given to Lattice::boundary_cartesian"); } - // Loop over boundary points. The order is - // ix = 0, iy = 0 to Ny-1, iz = 0 to Nz-1 - // iy = 0, ix = 1 to Nx-1, iz = 0 to Nz-1 - // iz = 0, ix = 1 to Nx-1, iy = 1 to Ny-1 - // giving Ny*Nz + (Nx-1)*Nz + (Nx-1)*(Ny-1) = Ny*Nz+Nx*Nz+Nx*Ny-Nz-Nx-Ny+1 points as expected. + + // Loop over boundary points. The order is (with k = 2b+1) + // x-boundary k*Ny*Nz points + // y-boundary k*(Nx-k)*Nz points + // z-boundary k*(Nx-k)*(Ny-k) points + // giving a total of k*(Ny*Nz + Nx*Nz + Nx*Ny - k*Nz - k*Ny-k*Nx+k*k) points + // as expected. bool get_next_boundary_point(int &ix, int &iy, int &iz) const { - if(iz == 0 && ix != 0 && iy != 0) // on the last boundary + int b = boundary_width_; + + // are we in x-boundary ? + if(ix <= b || ix >= Nx_-b) { - if(iy < Ny_-1) {iy++; return true;} - else if(ix < Nx_-1) {iy = 1; ix++; return true;} - return false; // finished + if(iz < Nz_-1) {iz++; return true;} + if(iy < Ny_-1) {iy++; iz = 0; return true;} + if(ix < b) {ix++; iy=iz=0; return true;} + if(ix == b && b>0) {ix = Nx_-b; iy=iz=0; return true;} + if(ix > b && ix < Nx_-1) {ix++; iy=iz=0; return true;} + ix = b+1; iy=iz=0; return true; // go to y boundary } - if(iy == 0 && ix != 0) + + // are we in y-boundary? + // We cannot visit x-boundary points + if(iy <= b || iy >= Ny_-b) { - if(iz < Nz_-1) {iz++; return true;} - else if(ix < Nx_-1) {iz = 0; ix++; return true;} - else {iz = 0; ix = iy = 1; return true;} // start z-boundary - } - if(ix == 0) + if(iz < Nz_-1) {iz++; return true;} + if(ix < Nx_-b-1) {ix++; iz = 0; return true;} + if(iy < b) {iy++; ix = b+1; iz = 0; return true;} + if(iy == b && b>0) {iy=Ny_-b;ix = b+1; iz = 0; return true;} + if(iy > b && iy < Ny_-1) {iy++; ix = b+1; iz = 0; return true;} + ix = iy = b+1; iz = 0; return true; // go to z-boundary + } + // are we on z-boundary? + // we cannot visit x or y boundary points + if(iz <= b || iz >=Nz_-b) { - if(iz < Nz_-1) {iz++; return true;} - else if(iy < Ny_-1) { iy++; iz = 0; return true;} - else {ix = 1; iy = 0; iz = 0; return true;} + if(iy < Ny_-b-1) {iy++; return true;} + if(ix < Nx_-b-1) {ix++; iy = b+1; return true;} + if(iz < b) {iz++; ix = b+1; iy = b+1; return true;} + if(iz == b && b>0) {iz=Nz_-b;ix = b+1; iy = b+1; return true;} + if(iz > b && iz < Nz_-1) {iz++; ix = b+1; iy = b+1; return true;} } + // no where left to go or we were fed a bad starting point throw std::runtime_error("Lattice::get_next_boundary_point called with bad arguments"); } @@ -331,6 +376,8 @@ void init(double L[]) in.read((char*) &l.dz_, sizeof(double)); in.read((char*) &l.L_, 3*sizeof(double)); + + in.read((char*) &l.boundary_width_, sizeof(unsigned)); return in; } @@ -349,6 +396,8 @@ void init(double L[]) ar & dz_; ar & L_; + if(version > 0) + ar & boundary_width_; } protected: @@ -364,6 +413,10 @@ void init(double L[]) double dz_; ///< Lattice spacing in z direction double L_[3]; ///< Dimensions of physical box in hard-sphere units + + unsigned boundary_width_ = 0; }; +BOOST_CLASS_VERSION(Lattice, 1) + #endif diff --git a/include/Species.h b/include/Species.h index 436001e..d0df0eb 100644 --- a/include/Species.h +++ b/include/Species.h @@ -13,6 +13,15 @@ #include "External_Field.h" #include "EOS.h" +/* Constraints: + 1. Fixed Mass : total mass is not allowed to change: value is either given or taken as current mass. + 2. Homogeneous Boundary : all boundary points forced to have same density + 3. Fixed Background : background(i.e. border) points are held constant and never change. +*/ + + + + class Species { public: @@ -29,14 +38,13 @@ class Species void setFixedMass(double m) { set_fixed_mass(m);} void setFixedBackground(bool fixed) { set_open_system(fixed);} void setHomogeneousBoundary(bool val) { set_homogeneous_boundary(val);} - - // These functions should only take care of the fixed mass and it should be the - // role of another function to setup the physical system (e.g. set_closed_system). + + // fixed mass void set_fixed_mass_only() { fixedMass_ = density_->get_mass(); mu_ = 0.0; } void set_fixed_mass() { fixedMass_ = density_->get_mass(); mu_ = 0.0; fixedBackground_ = false; } void set_fixed_mass(double m) { fixedMass_ = m; if(m > 0.0) {mu_ = 0.0; fixedBackground_ = false;} } - // These functions should only take care of the boundary conditions + // boundary conditions void set_homogeneous_boundary(bool val) { homogeneousBoundary_ = val; } void set_fixed_boundary(bool val) { fixedBackground_ = val; } void set_fixed_background(bool val) { set_fixed_boundary(val); } diff --git a/src/Density.cpp b/src/Density.cpp index 0fc876b..1619349 100644 --- a/src/Density.cpp +++ b/src/Density.cpp @@ -268,7 +268,7 @@ void Density::shrink(double distance, double width, double gap_size) double R2 = distance + gap_size; for(int p=0; p boundary_width_ && iz < Nz_-boundary_width_) continue; + pos++; + + long p = boundary_pos(ix,iy,iz); + if(p != pos) + { + cout << "\t" << ix << " " << iy << " " << iz << " " << pos << " " << p << endl; + throw std::runtime_error("p1"); + } + + int jx,jy,jz; jx = jy = jz = -1; + boundary_cartesian(pos,jx,jy,jz); + if(ix != ix || jy != iy || jz != iz) + { + cout << "\t" << ix << " " << iy << " " << iz << " " << pos << " " << jx << " " << jy << " " << jz << endl; + throw std::runtime_error("p2"); + } + } + + // y-boundary + for(iy = 0; iy < Ny_; iy++) + for(ix = 0;ix < Nx_; ix++) + for(iz = boundary_width_+1; iz < Nz_-boundary_width_; iz++) + { + if(iy > boundary_width_ && iy < Ny_-boundary_width_) continue; + pos++; + + long p = boundary_pos(ix,iy,iz); + if(p != pos) + { + cout << "\t" << ix << " " << iy << " " << iz << " " << pos << " " << p << endl; + throw std::runtime_error("p3"); + } + + int jx,jy,jz; jx = jy = jz = -1; + boundary_cartesian(pos,jx,jy,jz); + if(ix != ix || jy != iy || jz != iz) + { + cout << "\t" << ix << " " << iy << " " << iz << " " << pos << " " << jx << " " << jy << " " << jz << endl; + throw std::runtime_error("p4"); + } + } + + // x-boundary + for(ix = 0;ix < Nx_; ix++) + for(iy = boundary_width_+1;iy < Ny_-boundary_width_; iy++) + for(iz = boundary_width_+1; iz < Nz_-boundary_width_; iz++) + { + if(ix > boundary_width_ && ix < Nx_-boundary_width_) continue; + pos++; + + long p = boundary_pos(ix,iy,iz); + if(p != pos) + { + cout << "\t" << ix << " " << iy << " " << iz << " " << pos << " " << p << endl; + throw std::runtime_error("p5"); + } + + int jx,jy,jz; jx = jy = jz = -1; + boundary_cartesian(pos,jx,jy,jz); + if(ix != ix || jy != iy || jz != iz) + { + cout << "\t" << ix << " " << iy << " " << iz << " " << pos << " " << jx << " " << jy << " " << jz << endl; + throw std::runtime_error("p6"); + } + } + cout << "\tTest of conversion to and from boundary position passed" << endl; + + + ix = iy = iz = 0; + int count = 1; + while(1) + { + try {get_next_boundary_point(ix,iy,iz);} //cout << "\t" << ix << " " << iy << " " << iz << " : " << count << endl;} + catch(...) {break;} + count++; + } + cout << "\tFound " << count << " boundary points calling get_next_boundary_point : expected " << get_Nboundary() << endl; + if(count != get_Nboundary()) {cout << endl; throw std::runtime_error("get_next_boundary_point error\n\n");} } diff --git a/src/Species.cpp b/src/Species.cpp index aa2dc02..b909e1f 100644 --- a/src/Species.cpp +++ b/src/Species.cpp @@ -163,15 +163,16 @@ void Species::beginForceCalculation() double Mtot = density_->getNumberAtoms(); double Mboundary = density_->get_ave_background_density() * density_->dV() * density_->get_Nboundary(); double Mtarget = fixedMass_; - - for (long p=0;pNtot();p++) if (!density_->is_boundary_point(p)) - { - double d = density_->get(p); - d *= (Mtarget-Mboundary)/(Mtot-Mboundary); - density_->set(p,d); - } - } - else *density_ *= (fixedMass_/density_->getNumberAtoms()); + + // TBD: cann't we get rid of the if-else? + for (long p=0;pNtot();p++) + if (!density_->is_boundary_point(p)) + { + double d = density_->get(p); + d *= (Mtarget-Mboundary)/(Mtot-Mboundary); + density_->set(p,d); + } + } else *density_ *= (fixedMass_/density_->getNumberAtoms()); mu_ = 0.0; } @@ -181,6 +182,10 @@ void zero_background_forces() { } + +// Makes adjustments to enforce fixed background (forces on background set to zero) +// or homogeneous background (forces on background set to same value) +// and adjusts forces for fixed mass calculation. double Species::endForceCalculation() { @@ -202,8 +207,9 @@ double Species::endForceCalculation() for(long pos = 0; pos < density_->get_Nboundary(); pos++) dF_.set(density_->boundary_pos_2_pos(pos),average_border_force); } - - + + // Next two blocks should really be combined: there is no advantage to separating them. + // Also, it does not seem we need to worry about setting boundary forces to zero (again). if(fixedMass_ > 0.0 && !fixedBackground_) { mu_ = 0.0; @@ -225,14 +231,21 @@ double Species::endForceCalculation() double Mboundary = density_->get_ave_background_density() * density_->dV() * density_->get_Nboundary(); double Mtarget = fixedMass_ - Mboundary; - for(long p=0;pNtot();p++) if (!density_->is_boundary_point(p)) + for(long p=0;pNtot();p++) + if (!density_->is_boundary_point(p)) mu_ += dF_.get(p)*density_->get(p); + mu_ /= Mtarget; - for(long p=0;pNtot();p++) if (!density_->is_boundary_point(p)) - dF_.set(p, dF_.get(p)-mu_*density_->dV()); - - for(long pos = 0; pos < density_->get_Nboundary(); pos++) - dF_.set(density_->boundary_pos_2_pos(pos),0.0); + + for(long p=0;pNtot();p++) + { + if (!density_->is_boundary_point(p)) + dF_.set(p, dF_.get(p)-mu_*density_->dV()); + else dF_.set(p, 0.0); + } + // added "else" above and eliminated following in order to be more efficient + // for(long pos = 0; pos < density_->get_Nboundary(); pos++) + // dF_.set(density_->boundary_pos_2_pos(pos),0.0); }