Minor updates to field coefficient bases and OutVel

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features
 
 project(
   EXP
-  VERSION "7.8.4"
+  VERSION "7.8.5"
   HOMEPAGE_URL https://github.com/EXP-code/EXP
   LANGUAGES C CXX Fortran)
 

expui/FieldBasis.cc

@@ -33,6 +33,7 @@ namespace BasisClasses
     "ascl",
     "delta",
     "lmax",
+    "mmax",
     "nmax",
     "model"
   };
@@ -54,7 +55,7 @@ namespace BasisClasses
 
     // Allocate coefficient storage
     //
-    nfld = p.size() + 2;
+    nfld = p.size() + 1;  // Add the density to the phase-space return
     allocateStore();
 
     // Okay to register
@@ -67,7 +68,7 @@ namespace BasisClasses
   //
   void FieldBasis::configure()
   {
-    nfld      = 2;		// Weight and density fields by default
+    nfld      = 1;		// Density field by default
     lmax      = 4;
     nmax      = 10;
     rmin      = 1.0e-4;
@@ -137,7 +138,8 @@ namespace BasisClasses
 
       // Compute interpolation functionoid
       //
-      double rmin = r.front(), rmax = r.back();
+      rmin = r.front();
+      rmax = r.back();
 
       interp = std::make_shared<Linear1d>(r, d);
       densfunc = [this](double r)
@@ -163,12 +165,15 @@ namespace BasisClasses
 
     // Generate the orthogonal function instance
     //
-    ortho = std::make_shared<OrthoFunction>(nmax, densfunc, rmin, rmax, rmapping, dof);
+    ortho = std::make_shared<OrthoFunction>
+	      (nmax-1, densfunc, rmin, rmax, rmapping, dof);
+    //        ^
+    //        |
+    //        +--- This is the polynomial order, not the rank
 
     // Initialize fieldlabels
     //
     fieldLabels.clear();
-    fieldLabels.push_back("weight");
     fieldLabels.push_back("density");
 
     // Debug
@@ -204,9 +209,11 @@ namespace BasisClasses
 
   void FieldBasis::initialize()
   {
-    // Remove matched keys
+    // Check for unmatched keys
     //
-    // for (auto v : valid_keys) current_keys.erase(v);
+    auto unmatched = YamlCheck(conf, valid_keys);
+    if (unmatched.size())
+      throw YamlConfigError("Basis::FieldBasis", "parameter", unmatched, __FILE__, __LINE__);
 
     // Assign values from YAML
     //
@@ -215,6 +222,7 @@ namespace BasisClasses
       if (conf["model"    ])    model     = conf["model"    ].as<std::string>();
       if (conf["nfld"     ])    nfld      = conf["nfld"     ].as<int>();
       if (conf["lmax"     ])    lmax      = conf["lmax"     ].as<int>();
+      if (conf["mmax"     ])    lmax      = conf["mmax"     ].as<int>();
       if (conf["nmax"     ])    nmax      = conf["nmax"     ].as<int>();
       if (conf["dof"      ])    dof       = conf["dof"      ].as<int>();
       if (conf["rmin"     ])    rmin      = conf["rmin"     ].as<double>();
@@ -283,7 +291,7 @@ namespace BasisClasses
       // Sanity test dimensions
       if (nfld!=p->nfld || lmax!=p->mmax || nmax!=p->nmax) {
 	std::ostringstream serr;
-	serr << "FieldBasis::set_coefs: dimension error! "
+	serr << "FieldBasis::set_coefs: dimension error for dof=2! "
 	     << " nfld [" << nfld << "!= " << p->nfld << "]"
 	     << " mmax [" << lmax << "!= " << p->mmax << "]"
 	     << " nmax [" << nmax << "!= " << p->nmax << "]";
@@ -298,7 +306,7 @@ namespace BasisClasses
       // Sanity test dimensions
       if (nfld!=p->nfld || lmax!=p->lmax || nmax!=p->nmax) {
 	std::ostringstream serr;
-	serr << "FieldBasis::set_coefs: dimension error! "
+	serr << "FieldBasis::set_coefs: dimension error for dof=3! "
 	     << " nfld [" << nfld << "!= " << p->nfld << "]"
 	     << " lmax [" << lmax << "!= " << p->lmax << "]"
 	     << " nmax [" << nmax << "!= " << p->nmax << "]";
@@ -312,8 +320,7 @@ namespace BasisClasses
 			      double x, double y, double z,
 			      double u, double v, double w)
   {
-    constexpr std::complex<double> I(0, 1);
-    constexpr double fac0 = 0.25*M_2_SQRTPI;
+    constexpr double fac2 = 0.5*M_2_SQRTPI/M_SQRT2; // 1/sqrt(2*pi)=0.3989422804
 
     int tid = omp_get_thread_num();
     PS3 pos{x, y, z}, vel{u, v, w};
@@ -337,45 +344,42 @@ namespace BasisClasses
 
       auto p = (*ortho)(R);
 
-      (*coefs[tid])(0, 0, 0) += mass*p(0)*fac0;
-
       for (int m=0; m<=lmax; m++) {
 
-	std::complex<double> P = std::exp(-I*(phi*m));
+	auto P = std::exp(std::complex<double>(0, -phi*m))*fac2;
 
 	for (int n=0; n<nmax; n++) {
 
-	  (*coefs[tid])(1, m, n) += mass*P*p(n);
+	  (*coefs[tid])(0, m, n) += mass*P*p(n);
 
 	  for (int k=0; k<vec.size(); k++)
-	    (*coefs[tid])(k+2, m, n) += mass*P*p(n)*vec[k];
+	    (*coefs[tid])(k+1, m, n) += mass*P*p(n)*vec[k];
 	}
       }	 
 
     } else {
 
       auto p = (*ortho)(r);
 
-      (*coefs[tid])(0, 0, 0) += mass*p(0);
-
       for (int l=0, lm=0; l<=lmax; l++) {
 
 	double s = 1.0;
 
 	for (int m=0; m<=l; m++, lm++) {
 
 	  // Spherical harmonic value
-	  std::complex<double> P =
-	    std::exp(-I*(phi*m))*Ylm(l, m)*plgndr(l, m, cth) * s;
-
-	  s *= -1.0;		// Flip sign for next evaluation
+	  auto P = std::exp(std::complex<double>(0, -phi*m)) *
+	    Ylm(l, m)*plgndr(l, m, cth) * s;
+
+	  // Flip sign for next evaluation
+	  s *= -1.0;
 
 	  for (int n=0; n<nmax; n++) {
 
-	    (*coefs[tid])(1, lm, n) += mass*P*p(n);
+	    (*coefs[tid])(0, lm, n) += mass*P*p(n);
 
 	    for (int k=0; k<vec.size(); k++)
-	      (*coefs[tid])(k+2, lm, n) += mass*P*p(n)*vec[k];
+	      (*coefs[tid])(k+1, lm, n) += mass*P*p(n)*vec[k];
 	  }
 	}
       }
@@ -407,8 +411,10 @@ namespace BasisClasses
     //
     if (not coefret) {
       if (dof==2) coefret = std::make_shared<CoefClasses::CylFldStruct>();
-      else        coefret = std::make_shared<CoefClasses::SphFldStruct  >();
-				// Add the configuration data
+      else        coefret = std::make_shared<CoefClasses::SphFldStruct>();
+
+      // Add the configuration data
+      //
       std::ostringstream sout; sout << node;
       coefret->buf = sout.str();
     }
@@ -443,18 +449,21 @@ namespace BasisClasses
   std::vector<double>
   FieldBasis::sph_eval(double r, double costh, double phi)
   {
-    constexpr std::complex<double> I(0, 1);
+    constexpr double fac2 = 0.5*M_2_SQRTPI/M_SQRT2; // 1/sqrt(2 pi)
 
     if (dof==2) {
       std::vector<double> ret(nfld, 0);
-      auto p = (*ortho)(r*sqrt(fabs(1.0 - costh*costh)));
+      auto p = (*ortho)(r);
       for (int i=0; i<nfld; i++) {
 	for (int m=0; m<=lmax; m++) {
-	  double fac = m>0 ? 2.0 : 1.0;
+	  auto P = std::exp(std::complex<double>(0, -phi*m));
+
+	  std::complex<double> sum = 0.0;
 	  for (int n=0; n<nmax; n++) {
-	    ret[i] += fac *
-	      std::real((*coefs[0])(i, m, n)*exp(I*(phi*m)))*p[n];
+	    sum += (*coefs[0])(i, m, n)*P * p(n) * fac2;
 	  }
+
+	  ret[i] += std::real(sum);
 	}
       }
       return ret;
@@ -478,9 +487,11 @@ namespace BasisClasses
       for (int i=0; i<nfld; i++) {
 	for (int l=0, L=0; l<=lmax; l++) {
 	  for (int m=0; m<=l; m++, L++) {
+	    auto P = std::exp(std::complex<double>(0, -phi*m));
+
 	    double fac = m>0 ? 2.0 : 1.0;
 	    for (int n=0; n<nmax; n++) {
-	      ret[i] += std::real((*coefs[0])(i, L, n)*exp(I*(phi*m)))*p[n]*
+	      ret[i] += std::real((*coefs[0])(i, L, n)*P)*p(n)*
 		Ylm(l, m)*plgndr(l, m, costh) * fac;
 	    }
 	  }
@@ -751,7 +762,6 @@ namespace BasisClasses
 
     // Field labels (force field labels added below)
     //
-    fieldLabels.push_back("weight");
     fieldLabels.push_back("density");
 
     if (coordinates == Coord::Cylindrical) {

exputil/OrthoFunction.cc

@@ -2,8 +2,8 @@
 #include <OrthoFunction.H>
 
 OrthoFunction::OrthoFunction
-(int nmax, DensFunc W, double rmin, double rmax, double scl, int dof, bool seg) :
-  nmax(nmax), W(W), rmin(rmin), rmax(rmax), scale(scl), dof(dof), segment(seg)
+(int nmax, DensFunc W, double rmin, double rmax, double scl, int dof) :
+  nmax(nmax), W(W), rmin(rmin), rmax(rmax), scale(scl), dof(dof)
 {
   // Quadrature knots/weights
   lq = std::make_shared<LegeQuad>(knots);
@@ -27,14 +27,14 @@ void OrthoFunction::generate()
   norm.resize(nmax+1);	
 
   // Initial values
-  beta[0] = norm[0] = scalar_prod(0, 0);
-  alph[0] = scalar_prod(0, 1)/norm[0];
+  beta(0) = norm(0) = scalar_prod(0, 0);
+  alph(0) = scalar_prod(0, 1)/norm(0);
 
   // Remaining values
   for (int i=1; i<=nmax; i++) {
-    norm[i] = scalar_prod(i, 0);
-    alph[i] = scalar_prod(i, 1)/norm[i];
-    beta[i] = norm[i]/norm[i-1];
+    norm(i) = scalar_prod(i, 0);
+    alph(i) = scalar_prod(i, 1)/norm(i);
+    beta(i) = norm(i)/norm(i-1);
   }
 }
 
@@ -45,10 +45,10 @@ double OrthoFunction::scalar_prod(int n, int m)
   for (int i=0; i<knots; i++) {
     double x = xmin + dx*lq->knot(i);
     double r = x_to_r(x);
-    double y = 2.0*r/scale; if (segment) y = x;
-    Eigen::VectorXd p = poly_eval(y, n) * W(r);
+    Eigen::VectorXd p = poly_eval(r, n) * W(r);
+
     ans += dx*lq->weight(i)*d_r_to_x(x)*pow(r, dof-1) *
-      p(n) * p(n) * pow(y, m);
+      p(n) * p(n) * pow(r, m);
   }
 
   return ans;
@@ -58,11 +58,11 @@ double OrthoFunction::scalar_prod(int n, int m)
 Eigen::VectorXd OrthoFunction::poly_eval(double t, int n)
 {
   Eigen::VectorXd ret(n+1);
-  ret[0] = 1.0;
+  ret(0) = 1.0;
   if (n) {
-    ret[1] = (t - alph[0])*ret[0];
+    ret(1) = (t - alph(0))*ret(0);
     for (int j=1; j<n; j++) {
-      ret[j+1] = (t - alph[j])*ret[j] - beta[j]*ret[j-1];
+      ret(j+1) = (t - alph(j))*ret(j) - beta(j)*ret(j-1);
     }
   }
 
@@ -81,15 +81,14 @@ Eigen::MatrixXd OrthoFunction::testOrtho()
     double x = xmin + dx*lq->knot(i);
     double r = x_to_r(x);
     double f = dx*lq->weight(i)*d_r_to_x(x)*pow(r, dof-1);
-    double y = 2.0*r/scale; if (segment) y = x;
 
     // Evaluate the unnormalized polynomial
-    Eigen::VectorXd p = poly_eval(y, nmax) * W(r);
+    Eigen::VectorXd p = poly_eval(r, nmax) * W(r);
 
     // Compute scalar products with the normalizations
     for (int n1=0; n1<=nmax; n1++) {
       for (int n2=0; n2<=nmax; n2++) {
-	ret(n1, n2) += f * p(n1) * p(n2) / sqrt(norm[n1]*norm[n2]);
+	ret(n1, n2) += f * p(n1) * p(n2) / sqrt(norm(n1)*norm(n2));
       }
     }
   }
@@ -111,13 +110,12 @@ void OrthoFunction::dumpOrtho(const std::string& filename)
       //
       double x = xmin + dx*i/(number-1);
       double r = x_to_r(x);
-      double y = 2.0*r/scale; if (segment) y = x;
 
       // Evaluate polynomial and apply the normalization
       //
-      Eigen::VectorXd p = poly_eval(y, nmax) * W(r);
+      Eigen::VectorXd p = poly_eval(r, nmax) * W(r);
       fout << std::setw(16) << r;
-      for (int n=0; n<nmax; n++) fout << std::setw(16) << p(n)/sqrt(norm[n]);
+      for (int n=0; n<=nmax; n++) fout << std::setw(16) << p(n)/sqrt(norm(n));
       fout << std::endl;
     }
   }
@@ -131,13 +129,9 @@ Eigen::VectorXd OrthoFunction::operator()(double r)
   // Weight
   double w = W(r);
 
-  // Choose system;
-  double y = 2.0*r/scale;
-  if (segment) y = r_to_x(r);
-
   // Generate normalized orthogonal functions with weighting
-  auto p = poly_eval(y, nmax);
-  for (int i=0; i<=nmax; i++) p(i) *= w/sqrt(norm[i]);
+  auto p = poly_eval(r, nmax);
+  for (int i=0; i<=nmax; i++) p(i) *= w/sqrt(norm(i));
 
   // The inner product of p(i), p(j) is Kronecker delta(i, j)
 

exputil/mestel.cc

@@ -158,7 +158,7 @@ double TaperedMestelDisk::get_mass(double R)
     for (int i=1; i<=N; i++) {
       double rr  = rmin + dr*i;
       double cur = get_density(rr);
-      cum += 0.5*dr*(lst + cur) * 2.0*M_PI*rr;
+      cum += 0.5*dr*(lst*(rr-dr) + cur*rr) * 2.0*M_PI;
       lst = cur;
 
       vecR(i) = rr;