#ifndef _SCTL_BOUNDARY_QUADRATURE_HPP_ #define _SCTL_BOUNDARY_QUADRATURE_HPP_ #include #include #include #include #include #include namespace SCTL_NAMESPACE { template class Basis { public: using ValueType = Real; // class EvalOperator { // public: // }; using EvalOpType = Matrix; static constexpr Long Dim() { return DIM; } static constexpr Long Size() { return pow(ORDER); } static const Matrix& Nodes() { static Matrix nodes_(DIM,Size()); auto nodes_1d = [](Integer i) { return 0.5 - 0.5 * sctl::cos((2*i+1) * const_pi() / (2*ORDER)); }; { // Set nodes_ static std::mutex mutex; static std::atomic first_time(true); if (first_time.load(std::memory_order_relaxed)) { std::lock_guard guard(mutex); if (first_time.load(std::memory_order_relaxed)) { Integer N = 1; for (Integer d = 0; d < DIM; d++) { for (Integer j = 0; j < ORDER; j++) { for (Integer i = 0; i < N; i++) { for (Integer k = 0; k < d; k++) { nodes_[k][j*N+i] = nodes_[k][i]; } nodes_[d][j*N+i] = nodes_1d(j); } } N *= ORDER; } std::atomic_thread_fence(std::memory_order_seq_cst); first_time.store(false); } } } return nodes_; } static const Vector& QuadWts() { static Vector wts(Size()); { // Set nodes_ static std::mutex mutex; static std::atomic first_time(true); if (first_time.load(std::memory_order_relaxed)) { std::lock_guard guard(mutex); if (first_time.load(std::memory_order_relaxed)) { StaticArray wts_1d; { // Set wts_1d Vector x_(ORDER); ChebBasis::template Nodes<1>(ORDER, x_); Vector V_cheb(ORDER * ORDER); { // Set V_cheb Vector I(ORDER*ORDER); I = 0; for (Long i = 0; i < ORDER; i++) I[i*ORDER+i] = 1; ChebBasis::template Approx<1>(ORDER, I, V_cheb); } Matrix M(ORDER, ORDER, V_cheb.begin()); Vector w_sample(ORDER); for (Integer i = 0; i < ORDER; i++) { w_sample[i] = (i % 2 ? 0 : -(ORDER/(ValueType)(i*i-1))); } for (Integer j = 0; j < ORDER; j++) { wts_1d[j] = 0; for (Integer i = 0; i < ORDER; i++) { wts_1d[j] += M[j][i] * w_sample[i] / ORDER; } } } wts[0] = 1; Integer N = 1; for (Integer d = 0; d < DIM; d++) { for (Integer j = 1; j < ORDER; j++) { for (Integer i = 0; i < N; i++) { wts[j*N+i] = wts[i] * wts_1d[j]; } } for (Integer i = 0; i < N; i++) { wts[i] *= wts_1d[0]; } N *= ORDER; } std::atomic_thread_fence(std::memory_order_seq_cst); first_time.store(false); } } } return wts; } static void Grad(Vector& dX, const Vector& X) { static Matrix GradOp[DIM]; static std::mutex mutex; static std::atomic first_time(true); if (first_time.load(std::memory_order_relaxed)) { std::lock_guard guard(mutex); if (first_time.load(std::memory_order_relaxed)) { { // Set GradOp auto nodes = Basis::Nodes(); SCTL_ASSERT(nodes.Dim(1) == ORDER); Matrix M(ORDER, ORDER); for (Integer i = 0; i < ORDER; i++) { // Set M Real x = nodes[0][i]; for (Integer j = 0; j < ORDER; j++) { M[j][i] = 0; for (Integer l = 0; l < ORDER; l++) { if (l != j) { Real M_ = 1; for (Integer k = 0; k < ORDER; k++) { if (k != j && k != l) M_ *= (x - nodes[0][k]); if (k != j) M_ /= (nodes[0][j] - nodes[0][k]); } M[j][i] += M_; } } } } for (Integer d = 0; d < DIM; d++) { GradOp[d].ReInit(Size(), Size()); GradOp[d] = 0; Integer stride0 = sctl::pow(ORDER, d); Integer repeat0 = sctl::pow(ORDER, d); Integer stride1 = sctl::pow(ORDER, d+1); Integer repeat1 = sctl::pow(ORDER, DIM-d-1); for (Integer k1 = 0; k1 < repeat1; k1++) { for (Integer i = 0; i < ORDER; i++) { for (Integer j = 0; j < ORDER; j++) { for (Integer k0 = 0; k0 < repeat0; k0++) { GradOp[d][k1*stride1 + i*stride0 + k0][k1*stride1 + j*stride0 + k0] = M[i][j]; } } } } } } std::atomic_thread_fence(std::memory_order_seq_cst); first_time.store(false); } } if (dX.Dim() != X.Dim()*DIM) dX.ReInit(X.Dim()*DIM); for (Long i = 0; i < X.Dim(); i++) { const Matrix Vi(1, Size(), (Iterator)(ConstIterator)X[i].NodeValues_, false); for (Integer k = 0; k < DIM; k++) { Matrix Vo(1, Size(), dX[i*DIM+k].NodeValues_, false); Matrix::GEMM(Vo, Vi, GradOp[k]); } } } static EvalOpType SetupEval(const Matrix& X) { Long N = X.Dim(1); SCTL_ASSERT(X.Dim(0) == DIM); Matrix M(Size(), N); { // Set M auto nodes = Basis::Nodes(); Integer NN = Basis::Size(); Matrix M_(NN, DIM*N); for (Long i = 0; i < DIM*N; i++) { ValueType x = X[0][i]; for (Integer j = 0; j < NN; j++) { ValueType y = 1; for (Integer k = 0; k < NN; k++) { y *= (j==k ? 1 : (nodes[0][k] - x) / (nodes[0][k] - nodes[0][j])); } M_[j][i] = y; } } if (DIM == 1) { SCTL_ASSERT(M.Dim(0) == M_.Dim(0)); SCTL_ASSERT(M.Dim(1) == M_.Dim(1)); M = M_; } else { Integer NNN = 1; M = 1; for (Integer d = 0; d < DIM; d++) { for (Integer k = 1; k < NN; k++) { for (Integer j = 0; j < NNN; j++) { for (Long i = 0; i < N; i++) { M[k*NNN+j][i] = M[j][i] * M_[k][d*N+i]; } } } { // k = 0 for (Integer j = 0; j < NNN; j++) { for (Long i = 0; i < N; i++) { M[j][i] *= M_[0][d*N+i]; } } } NNN *= NN; } } } return M; } static void Eval(Matrix& Y, const Vector& X, const EvalOpType& M) { Long N0 = X.Dim(); Long N1 = M.Dim(1); SCTL_ASSERT(M.Dim(0) == Size()); if (Y.Dim(0) != N0 || Y.Dim(1) != N1) Y.ReInit(N0, N1); for (Long i = 0; i < N0; i++) { const Matrix X_(1,Size(),(Iterator)(ConstIterator)X[i].NodeValues_,false); Matrix Y_(1,N1,Y[i],false); Matrix::GEMM(Y_,X_,M); } } Basis operator+(Basis X) const { for (Long i = 0; i < Size(); i++) X[i] = (*this)[i] + X[i]; return X; } Basis operator-(Basis X) const { for (Long i = 0; i < Size(); i++) X[i] = (*this)[i] - X[i]; return X; } Basis operator*(Basis X) const { for (Long i = 0; i < Size(); i++) X[i] = (*this)[i] * X[i]; return X; } Basis operator*(Real a) const { Basis X = (*this); for (Long i = 0; i < Size(); i++) X[i] *= a; return X; } Basis operator+(Real a) const { Basis X = (*this); for (Long i = 0; i < Size(); i++) X[i] += a; return X; } Basis& operator+=(const Basis& X) { for (Long i = 0; i < Size(); i++) (*this)[i] += X[i]; return *this; } Basis& operator-=(const Basis& X) { for (Long i = 0; i < Size(); i++) (*this)[i] -= X[i]; return *this; } Basis& operator*=(const Basis& X) { for (Long i = 0; i < Size(); i++) (*this)[i] *= X[i]; return *this; } Basis& operator*=(Real a) { for (Long i = 0; i < Size(); i++) (*this)[i] *= a; return *this; } Basis& operator+=(Real a) { for (Long i = 0; i < Size(); i++) (*this)[i] += a; return *this; } Basis& operator=(Real a) { for (Long i = 0; i < Size(); i++) (*this)[i] = a; return *this; } const ValueType& operator[](Long i) const { SCTL_ASSERT(i < Size()); return NodeValues_[i]; } ValueType& operator[](Long i) { SCTL_ASSERT(i < Size()); return NodeValues_[i]; } private: StaticArray NodeValues_; }; template class ElemList { public: using CoordBasis = Basis; using CoordType = typename CoordBasis::ValueType; static constexpr Integer CoordDim() { return COORD_DIM; } static constexpr Integer ElemDim() { return CoordBasis::Dim(); } ElemList(Long Nelem = 0) { ReInit(Nelem); } void ReInit(Long Nelem = 0) { Nelem_ = Nelem; X_.ReInit(Nelem_ * COORD_DIM); } void ReInit(const Vector& X) { Nelem_ = X.Dim() / COORD_DIM; SCTL_ASSERT(X.Dim() == Nelem_ * COORD_DIM); X_ = X; } Long NElem() const { return Nelem_; } CoordBasis& operator()(Long elem, Integer dim) { SCTL_ASSERT(elem >= 0 && elem < Nelem_); SCTL_ASSERT(dim >= 0 && dim < COORD_DIM); return X_[elem*COORD_DIM+dim]; } const CoordBasis& operator()(Long elem, Integer dim) const { if (!(elem >= 0 && elem < Nelem_)) exit(0); SCTL_ASSERT(elem >= 0 && elem < Nelem_); SCTL_ASSERT(dim >= 0 && dim < COORD_DIM); return X_[elem*COORD_DIM+dim]; } const Vector& ElemVector() const { return X_; } private: static_assert(CoordBasis::Dim() <= CoordDim(), "Basis dimension can not be greater than COORD_DIM."); Vector X_; Long Nelem_; //mutable Vector dX_; }; template class Quadrature { static Real machine_epsilon() { Real eps=1; while(eps*(Real)0.5+(Real)1.0>1.0) eps*=0.5; return eps; } template static void DuffyQuad(Matrix& nodes, Vector& weights, const Vector& coord, Integer order, Real adapt = -1.0) { SCTL_ASSERT(coord.Dim() == DIM); static Real eps = machine_epsilon()*16; Matrix qx; Vector qw; { // Set qx, qw Vector qx0, qw0; ChebBasis::quad_rule(order, qx0, qw0); Integer N = sctl::pow(order); qx.ReInit(DIM,N); qw.ReInit(N); qw[0] = 1; Integer N_ = 1; for (Integer d = 0; d < DIM; d++) { for (Integer j = 0; j < order; j++) { for (Integer i = 0; i < N_; i++) { for (Integer k = 0; k < d; k++) { qx[k][j*N_+i] = qx[k][i]; } qx[d][j*N_+i] = qx0[j]; qw[j*N_+i] = qw[i]; } } for (Integer j = 0; j < order; j++) { for (Integer i = 0; i < N_; i++) { qw[j*N_+i] *= qw0[j]; } } N_ *= order; } } Vector X; { // Set X StaticArray X_; X_[0] = 0; X_[1] = adapt; for (Integer i = 0; i < DIM; i++) { X_[2*i+2] = sctl::fabs(coord[i]); X_[2*i+3] = sctl::fabs(coord[i]-1); } std::sort((Iterator)X_, (Iterator)X_+2*DIM+2); X.PushBack(std::max(0, X_[2*DIM]-1)); for (Integer i = 0; i < 2*DIM+2; i++) { if (X[X.Dim()-1] < X_[i]) { if (X.Dim()) X.PushBack(X_[i]); } } ///////////////////////////////////////////////////////////////////////////////////////////////// Vector r(1); r[0] = X[0]; for (Integer i = 1; i < X.Dim(); i++) { while (r[r.Dim() - 1] > 0.0 && (order*0.5) * r[r.Dim() - 1] < X[i]) r.PushBack((order*0.5) * r[r.Dim() - 1]); // TODO r.PushBack(X[i]); } X = r; ///////////////////////////////////////////////////////////////////////////////////////////////// } Vector nds, wts; for (Integer k = 0; k < X.Dim()-1; k++) { for (Integer dd = 0; dd < 2*DIM; dd++) { Integer d0 = (dd>>1); StaticArray range0, range1; { // Set range0, range1 Integer d1 = (dd%2?1:-1); for (Integer d = 0; d < DIM; d++) { range0[d*2+0] = std::max(0,std::min(1,coord[d] - X[k] )); range0[d*2+1] = std::max(0,std::min(1,coord[d] + X[k] )); range1[d*2+0] = std::max(0,std::min(1,coord[d] - X[k+1])); range1[d*2+1] = std::max(0,std::min(1,coord[d] + X[k+1])); } range0[d0*2+0] = std::max(0,std::min(1,coord[d0] + d1*X[k+0])); range0[d0*2+1] = std::max(0,std::min(1,coord[d0] + d1*X[k+0])); range1[d0*2+0] = std::max(0,std::min(1,coord[d0] + d1*X[k+1])); range1[d0*2+1] = std::max(0,std::min(1,coord[d0] + d1*X[k+1])); } { // if volume(range0, range1) == 0 then continue Real v0 = 1, v1 = 1; for (Integer d = 0; d < DIM; d++) { if (d == d0) { v0 *= sctl::fabs(range0[d*2+0]-range1[d*2+0]); v1 *= sctl::fabs(range0[d*2+0]-range1[d*2+0]); } else { v0 *= range0[d*2+1]-range0[d*2+0]; v1 *= range1[d*2+1]-range1[d*2+0]; } } if (v0 < eps && v1 < eps) continue; } for (Integer i = 0; i < qx.Dim(1); i++) { // Set nds, wts Real w = qw[i]; Real z = qx[d0][i]; for (Integer d = 0; d < DIM; d++) { Real y = qx[d][i]; nds.PushBack((range0[d*2+0]*(1-y) + range0[d*2+1]*y)*(1-z) + (range1[d*2+0]*(1-y) + range1[d*2+1]*y)*z); if (d == d0) { w *= abs(range1[d*2+0] - range0[d*2+0]); } else { w *= (range0[d*2+1] - range0[d*2+0])*(1-z) + (range1[d*2+1] - range1[d*2+0])*z; } } wts.PushBack(w); } } } nodes = Matrix(nds.Dim()/DIM,DIM,nds.begin()).Transpose(); weights = wts; } template static void TensorProductGaussQuad(Matrix& nodes, Vector& weights, Integer order) { Vector coord(DIM); coord = 0; coord[0] = -10; DuffyQuad(nodes, weights, coord, order); } template static void SetupSingular(Matrix& M_singular, const Matrix& trg_nds, const ElemList& elem_lst, const Kernel& kernel, Integer order_singular = 10, Integer order_direct = 10, Real Rqbx = 0) { using CoordBasis = typename ElemList::CoordBasis; using CoordEvalOpType = typename CoordBasis::EvalOpType; using DensityEvalOpType = typename DensityBasis::EvalOpType; constexpr Integer CoordDim = ElemList::CoordDim(); constexpr Integer ElemDim = ElemList::ElemDim(); constexpr Integer KDIM0 = Kernel::SrcDim(); constexpr Integer KDIM1 = Kernel::TrgDim(); const Long Nelem = elem_lst.NElem(); const Integer Ntrg = trg_nds.Dim(1); SCTL_ASSERT(trg_nds.Dim(0) == ElemDim); const Vector& X = elem_lst.ElemVector(); Vector dX; CoordBasis::Grad(dX, X); Vector Xt, Xnt; { // Set Xt, Xnt auto Meval = CoordBasis::SetupEval(trg_nds); eval_basis(Xt, X, CoordDim, trg_nds.Dim(1), Meval); Xnt = Xt; Vector dX_; eval_basis(dX_, dX, 2*CoordDim, trg_nds.Dim(1), Meval); for (Long i = 0; i < Ntrg; i++) { for (Long j = 0; j < Nelem; j++) { auto Xn = Xnt.begin() + (j*Ntrg+i)*CoordDim; auto dX0 = dX_.begin() + (j*Ntrg+i)*2*CoordDim; StaticArray normal; normal[0] = dX0[2]*dX0[5] - dX0[4]*dX0[3]; normal[1] = dX0[4]*dX0[1] - dX0[0]*dX0[5]; normal[2] = dX0[0]*dX0[3] - dX0[2]*dX0[1]; Real Xa = sctl::sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]); Real invXa = 1/Xa; normal[0] *= invXa; normal[1] *= invXa; normal[2] *= invXa; Real sqrt_Xa = sqrt(Xa); Xn[0] = normal[0]*sqrt_Xa*Rqbx; Xn[1] = normal[1]*sqrt_Xa*Rqbx; Xn[2] = normal[2]*sqrt_Xa*Rqbx; } } } SCTL_ASSERT(Xt.Dim() == Nelem * Ntrg * CoordDim); auto& M = M_singular; M.ReInit(Nelem * KDIM0 * DensityBasis::Size(), KDIM1 * Ntrg); #pragma omp parallel for schedule(static) for (Long i = 0; i < Ntrg; i++) { // Set M (singular) Matrix quad_nds; Vector quad_wts; { // Set quad_nds, quad_wts StaticArray trg_node_; for (Integer k = 0; k < ElemDim; k++) { trg_node_[k] = trg_nds[k][i]; } Vector trg_node(ElemDim, trg_node_, false); DuffyQuad(quad_nds, quad_wts, trg_node, order_singular, fabs(Rqbx)); } const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds); Integer Nnds = quad_wts.Dim(); Vector X_, dX_, Xa_, Xn_; { // Set X_, dX_ eval_basis(X_, X, CoordDim, Nnds, CoordEvalOp); eval_basis(dX_, dX, CoordDim * ElemDim, Nnds, CoordEvalOp); } if (CoordDim == 3 && ElemDim == 2) { // Compute Xa_, Xn_ Long N = Nelem*Nnds; Xa_.ReInit(N); Xn_.ReInit(N*CoordDim); for (Long j = 0; j < N; j++) { StaticArray normal; normal[0] = dX_[j*6+2]*dX_[j*6+5] - dX_[j*6+4]*dX_[j*6+3]; normal[1] = dX_[j*6+4]*dX_[j*6+1] - dX_[j*6+0]*dX_[j*6+5]; normal[2] = dX_[j*6+0]*dX_[j*6+3] - dX_[j*6+2]*dX_[j*6+1]; Xa_[j] = sctl::sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]); Real invXa = 1/Xa_[j]; Xn_[j*3+0] = normal[0] * invXa; Xn_[j*3+1] = normal[1] * invXa; Xn_[j*3+2] = normal[2] * invXa; } } DensityEvalOpType DensityEvalOp; if (std::is_same::value) { DensityEvalOp = CoordEvalOp; } else { DensityEvalOp = DensityBasis::SetupEval(quad_nds); } for (Long j = 0; j < Nelem; j++) { Matrix M__(Nnds * KDIM0, KDIM1); if (Rqbx == 0) { // Set kernel matrix M__ const Vector X0_(CoordDim, Xt.begin() + (j * Ntrg + i) * CoordDim, false); const Vector X__(Nnds * CoordDim, X_.begin() + j * Nnds * CoordDim, false); const Vector Xn__(Nnds * CoordDim, Xn_.begin() + j * Nnds * CoordDim, false); kernel.template KernelMatrix(M__, X0_, X__, Xn__); } else { Vector X0_(CoordDim); constexpr Integer qbx_order = 6; StaticArray,qbx_order> M___; for (Integer k = 0; k < qbx_order; k++) { // Set kernel matrix M___ for (Integer kk = 0; kk < CoordDim; kk++) X0_[kk] = Xt[(j * Ntrg + i) * CoordDim + kk] + (k+1) * Xnt[(j * Ntrg + i) * CoordDim + kk]; const Vector X__(Nnds * CoordDim, X_.begin() + j * Nnds * CoordDim, false); const Vector Xn__(Nnds * CoordDim, Xn_.begin() + j * Nnds * CoordDim, false); kernel.template KernelMatrix(M___[k], X0_, X__, Xn__); } for (Long k = 0; k < Nnds * KDIM0 * KDIM1; k++) { M__[0][k] = 0; M__[0][k] += 6*M___[0][0][k]; M__[0][k] += -15*M___[1][0][k]; M__[0][k] += 20*M___[2][0][k]; M__[0][k] += -15*M___[3][0][k]; M__[0][k] += 6*M___[4][0][k]; M__[0][k] += -1*M___[5][0][k]; } } for (Long k0 = 0; k0 < KDIM0; k0++) { for (Long k1 = 0; k1 < KDIM1; k1++) { for (Long l = 0; l < DensityBasis::Size(); l++) { Real M_lk = 0; for (Long n = 0; n < Nnds; n++) { Real quad_wt = Xa_[j * Nnds + n] * quad_wts[n]; M_lk += DensityEvalOp[l][n] * quad_wt * M__[n*KDIM0+k0][k1]; } M[(j * KDIM0 + k0) * DensityBasis::Size() + l][k1 * Ntrg + i] = M_lk; } } } } } { // Set M (subtract direct) Matrix quad_nds; Vector quad_wts; TensorProductGaussQuad(quad_nds, quad_wts, order_direct); const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds); Integer Nnds = quad_wts.Dim(); Vector X_, dX_, Xa_, Xn_; { // Set X_, dX_ eval_basis(X_, X, CoordDim, Nnds, CoordEvalOp); eval_basis(dX_, dX, CoordDim * ElemDim, Nnds, CoordEvalOp); } if (CoordDim == 3 && ElemDim == 2) { // Compute Xa_, Xn_ Long N = Nelem*Nnds; Xa_.ReInit(N); Xn_.ReInit(N*CoordDim); for (Long j = 0; j < N; j++) { StaticArray normal; normal[0] = dX_[j*6+2]*dX_[j*6+5] - dX_[j*6+4]*dX_[j*6+3]; normal[1] = dX_[j*6+4]*dX_[j*6+1] - dX_[j*6+0]*dX_[j*6+5]; normal[2] = dX_[j*6+0]*dX_[j*6+3] - dX_[j*6+2]*dX_[j*6+1]; Xa_[j] = sctl::sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]); Real invXa = 1/Xa_[j]; Xn_[j*3+0] = normal[0] * invXa; Xn_[j*3+1] = normal[1] * invXa; Xn_[j*3+2] = normal[2] * invXa; } } DensityEvalOpType DensityEvalOp; if (std::is_same::value) { DensityEvalOp = CoordEvalOp; } else { DensityEvalOp = DensityBasis::SetupEval(quad_nds); } #pragma omp parallel for schedule(static) for (Long i = 0; i < Ntrg; i++) { // Subtract direct contribution for (Long j = 0; j < Nelem; j++) { Matrix M__(Nnds * KDIM0, KDIM1); { // Set kernel matrix M__ const Vector X0_(CoordDim, (Iterator)Xt.begin() + (j * Ntrg + i) * CoordDim, false); const Vector X__(Nnds * CoordDim, X_.begin() + j * Nnds * CoordDim, false); const Vector Xn__(Nnds * CoordDim, Xn_.begin() + j * Nnds * CoordDim, false); kernel.template KernelMatrix(M__, X0_, X__, Xn__); } for (Long k0 = 0; k0 < KDIM0; k0++) { for (Long k1 = 0; k1 < KDIM1; k1++) { for (Long l = 0; l < DensityBasis::Size(); l++) { Real M_lk = 0; for (Long n = 0; n < Nnds; n++) { Real quad_wt = Xa_[j * Nnds + n] * quad_wts[n]; M_lk += DensityEvalOp[l][n] * quad_wt * M__[n*KDIM0+k0][k1]; } M[(j * KDIM0 + k0) * DensityBasis::Size() + l][k1 * Ntrg + i] -= M_lk; } } } } } } } template static void EvalSingular(Matrix& U, const Vector& density, const Matrix& M, Integer KDIM0_, Integer KDIM1_) { if (M.Dim(0) == 0 || M.Dim(1) == 0) { U.ReInit(0,0); return; } const Long Ntrg = M.Dim(1) / KDIM1_; SCTL_ASSERT(M.Dim(1) == KDIM1_ * Ntrg); const Long Nelem = M.Dim(0) / (KDIM0_ * DensityBasis::Size()); SCTL_ASSERT(M.Dim(0) == Nelem * KDIM0_ * DensityBasis::Size()); const Integer dof = density.Dim() / (Nelem * KDIM0_); SCTL_ASSERT(density.Dim() == Nelem * dof * KDIM0_); if (U.Dim(0) != Nelem * dof * KDIM1_ || U.Dim(1) != Ntrg) { U.ReInit(Nelem * dof * KDIM1_, Ntrg); U = 0; } for (Long j = 0; j < Nelem; j++) { const Matrix M_(KDIM0_ * DensityBasis::Size(), KDIM1_ * Ntrg, (Iterator)M[j * KDIM0_ * DensityBasis::Size()], false); Matrix U_(dof, KDIM1_ * Ntrg, U[j*dof*KDIM1_], false); Matrix F_(dof, KDIM0_ * DensityBasis::Size()); for (Long i = 0; i < dof; i++) { for (Long k = 0; k < KDIM0_; k++) { for (Long l = 0; l < DensityBasis::Size(); l++) { F_[i][k * DensityBasis::Size() + l] = density[(j * dof + i) * KDIM0_ + k][l]; } } } Matrix::GEMM(U_, F_, M_); } } template struct PointData { bool operator<(const PointData& p) const { return mid < p.mid; } Long rank; Long surf_rank; Morton mid; StaticArray coord; Real radius2; }; template struct Pair { Pair() {} Pair(T1 x, T2 y) : first(x), second(y) {} bool operator<(const Pair& p) const { return (first < p.first) || (((first == p.first) && (second < p.second))); } T1 first; T2 second; }; template static void BuildNbrList(Vector>& pair_lst, const Vector& Xt, const Vector& trg_surf, const ElemList& elem_lst, Real distance_factor, Real period_length, const Comm& comm) { using CoordBasis = typename ElemList::CoordBasis; constexpr Integer CoordDim = ElemList::CoordDim(); constexpr Integer ElemDim = ElemList::ElemDim(); using PtData = PointData; const Integer rank = comm.Rank(); Real R0 = 0; StaticArray X0; { // Find bounding box Long N = Xt.Dim() / CoordDim; SCTL_ASSERT(Xt.Dim() == N * CoordDim); SCTL_ASSERT(N); StaticArray Xloc; StaticArray Xglb; for (Integer k = 0; k < CoordDim; k++) { Xloc[0*CoordDim+k] = Xt[k]; Xloc[1*CoordDim+k] = Xt[k]; } for (Long i = 0; i < N; i++) { for (Integer k = 0; k < CoordDim; k++) { Xloc[0*CoordDim+k] = std::min(Xloc[0*CoordDim+k], Xt[i*CoordDim+k]); Xloc[1*CoordDim+k] = std::max(Xloc[1*CoordDim+k], Xt[i*CoordDim+k]); } } comm.Allreduce((ConstIterator)Xloc+0*CoordDim, (Iterator)Xglb+0*CoordDim, CoordDim, Comm::CommOp::MIN); comm.Allreduce((ConstIterator)Xloc+1*CoordDim, (Iterator)Xglb+1*CoordDim, CoordDim, Comm::CommOp::MAX); for (Integer k = 0; k < CoordDim; k++) { R0 = std::max(R0, Xglb[1*CoordDim+k]-Xglb[0*CoordDim+k]); } R0 = R0 * 2.0; for (Integer k = 0; k < CoordDim; k++) { X0[k] = Xglb[k] - R0*0.25; } } if (period_length > 0) { R0 = period_length; } Vector PtSrc, PtTrg; Integer order_upsample = (Integer)(const_pi() / distance_factor + 0.5); { // Set PtSrc const Vector& X_elem_lst = elem_lst.ElemVector(); Vector dX_elem_lst; CoordBasis::Grad(dX_elem_lst, X_elem_lst); Matrix nds; Vector wts; TensorProductGaussQuad(nds, wts, order_upsample); const Long Nnds = nds.Dim(1); Vector X, dX; const auto CoordEvalOp = CoordBasis::SetupEval(nds); eval_basis(X, X_elem_lst, CoordDim, Nnds, CoordEvalOp); eval_basis(dX, dX_elem_lst, CoordDim * ElemDim, Nnds, CoordEvalOp); const Long N = X.Dim() / CoordDim; const Long Nelem = elem_lst.NElem(); SCTL_ASSERT(X.Dim() == N * CoordDim); SCTL_ASSERT(N == Nelem * Nnds); Long rank_offset, surf_rank_offset; { // Set rank_offset, surf_rank_offset comm.Scan(Ptr2ConstItr(&N,1), Ptr2Itr(&rank_offset,1), 1, Comm::CommOp::SUM); comm.Scan(Ptr2ConstItr(&Nelem,1), Ptr2Itr(&surf_rank_offset,1), 1, Comm::CommOp::SUM); surf_rank_offset -= Nelem; rank_offset -= N; } PtSrc.ReInit(N); const Real R0inv = 1.0 / R0; for (Long i = 0; i < N; i++) { // Set coord for (Integer k = 0; k < CoordDim; k++) { PtSrc[i].coord[k] = (X[i*CoordDim+k] - X0[k]) * R0inv; } } if (period_length > 0) { // Wrap-around coord for (Long i = 0; i < N; i++) { auto& x = PtSrc[i].coord; for (Integer k = 0; k < CoordDim; k++) { x[k] -= (Long)(x[k]); } } } for (Long i = 0; i < N; i++) { // Set radius2, mid, rank Integer depth = 0; { // Set radius2, depth Real radius2 = 0; for (Integer k0 = 0; k0 < ElemDim; k0++) { Real R2 = 0; for (Integer k1 = 0; k1 < CoordDim; k1++) { Real dX_ = dX[(i*CoordDim+k1)*ElemDim+k0]; R2 += dX_*dX_; } radius2 = std::max(radius2, R2); } radius2 *= R0inv*R0inv * distance_factor*distance_factor; PtSrc[i].radius2 = radius2; Long Rinv = (Long)(1.0/radius2); while (Rinv > 0) { Rinv = (Rinv>>2); depth++; } } PtSrc[i].mid = Morton((Iterator)PtSrc[i].coord, std::min(Morton::MaxDepth(),depth)); PtSrc[i].rank = rank_offset + i; } for (Long i = 0 ; i < Nelem; i++) { // Set surf_rank for (Long j = 0; j < Nnds; j++) { PtSrc[i*Nnds+j].surf_rank = surf_rank_offset + i; } } Vector PtSrcSorted; comm.HyperQuickSort(PtSrc, PtSrcSorted); PtSrc.Swap(PtSrcSorted); } { // Set PtTrg const Long N = Xt.Dim() / CoordDim; SCTL_ASSERT(Xt.Dim() == N * CoordDim); Long rank_offset; { // Set rank_offset comm.Scan(Ptr2ConstItr(&N,1), Ptr2Itr(&rank_offset,1), 1, Comm::CommOp::SUM); rank_offset -= N; } PtTrg.ReInit(N); const Real R0inv = 1.0 / R0; for (Long i = 0; i < N; i++) { // Set coord for (Integer k = 0; k < CoordDim; k++) { PtTrg[i].coord[k] = (Xt[i*CoordDim+k] - X0[k]) * R0inv; } } if (period_length > 0) { // Wrap-around coord for (Long i = 0; i < N; i++) { auto& x = PtTrg[i].coord; for (Integer k = 0; k < CoordDim; k++) { x[k] -= (Long)(x[k]); } } } for (Long i = 0; i < N; i++) { // Set radius2, mid, rank PtTrg[i].radius2 = 0; PtTrg[i].mid = Morton((Iterator)PtTrg[i].coord); PtTrg[i].rank = rank_offset + i; } if (trg_surf.Dim()) { // Set surf_rank SCTL_ASSERT(trg_surf.Dim() == N); for (Long i = 0; i < N; i++) { PtTrg[i].surf_rank = trg_surf[i]; } } else { for (Long i = 0; i < N; i++) { PtTrg[i].surf_rank = -1; } } Vector PtTrgSorted; comm.HyperQuickSort(PtTrg, PtTrgSorted); PtTrg.Swap(PtTrgSorted); } Tree tree(comm); { // Init tree Vector Xall(PtSrc.Dim()+PtTrg.Dim()); { // Set Xall Xall.ReInit((PtSrc.Dim()+PtTrg.Dim())*CoordDim); Long Nsrc = PtSrc.Dim(); Long Ntrg = PtTrg.Dim(); for (Long i = 0; i < Nsrc; i++) { for (Integer k = 0; k < CoordDim; k++) { Xall[i*CoordDim+k] = PtSrc[i].coord[k]; } } for (Long i = 0; i < Ntrg; i++) { for (Integer k = 0; k < CoordDim; k++) { Xall[(Nsrc+i)*CoordDim+k] = PtTrg[i].coord[k]; } } } tree.UpdateRefinement(Xall, 1000, true, period_length>0); } { // Repartition PtSrc, PtTrg PtData splitter; splitter.mid = tree.GetPartitionMID()[rank]; comm.PartitionS(PtSrc, splitter); comm.PartitionS(PtTrg, splitter); } { // Add tree data PtSrc const auto& node_mid = tree.GetNodeMID(); const Long N = node_mid.Dim(); SCTL_ASSERT(N); Vector dsp(N), cnt(N); for (Long i = 0; i < N; i++) { PtData m0; m0.mid = node_mid[i]; dsp[i] = std::lower_bound(PtSrc.begin(), PtSrc.end(), m0) - PtSrc.begin(); } for (Long i = 0; i < N-1; i++) { cnt[i] = dsp[i+1] - dsp[i]; } cnt[N-1] = PtSrc.Dim() - dsp[N-1]; tree.AddData("PtSrc", PtSrc, cnt); } tree.template Broadcast("PtSrc"); { // Build pair_lst Vector cnt; Vector PtSrc; tree.GetData(PtSrc, cnt, "PtSrc"); const auto& node_mid = tree.GetNodeMID(); const auto& node_attr = tree.GetNodeAttr(); Vector> nbr_mid_tmp; for (Long i = 0; i < node_mid.Dim(); i++) { if (node_attr[i].Leaf && !node_attr[i].Ghost) { Vector> child_mid; node_mid[i].Children(child_mid); for (const auto& trg_mid : child_mid) { Integer d0 = trg_mid.Depth(); Vector Src, Trg; { // Set Trg PtData m0, m1; m0.mid = trg_mid; m1.mid = trg_mid.Next(); Long a = std::lower_bound(PtTrg.begin(), PtTrg.end(), m0) - PtTrg.begin(); Long b = std::lower_bound(PtTrg.begin(), PtTrg.end(), m1) - PtTrg.begin(); Trg.ReInit(b-a, PtTrg.begin()+a, false); if (!Trg.Dim()) continue; } Vector> near_elem(Trg.Dim()); for (Integer d = 0; d <= d0; d++) { trg_mid.NbrList(nbr_mid_tmp, d, period_length>0); for (const auto& src_mid : nbr_mid_tmp) { // Set Src PtData m0, m1; m0.mid = src_mid; m1.mid = (d==d0 ? src_mid.Next() : src_mid.Ancestor(d+1)); Long a = std::lower_bound(PtSrc.begin(), PtSrc.end(), m0) - PtSrc.begin(); Long b = std::lower_bound(PtSrc.begin(), PtSrc.end(), m1) - PtSrc.begin(); Src.ReInit(b-a, PtSrc.begin()+a, false); if (!Src.Dim()) continue; for (Long t = 0; t < Trg.Dim(); t++) { // set near_elem[t] <-- {s : dist(s,t) < radius(s)} for (Long s = 0; s < Src.Dim(); s++) { if (Trg[t].surf_rank != Src[s].surf_rank) { Real R2 = 0; for (Integer k = 0; k < CoordDim; k++) { Real dx = (Src[s].coord[k] - Trg[t].coord[k]); R2 += dx * dx; } if (R2 < Src[s].radius2) { near_elem[t].insert(Src[s].surf_rank); } } } } } } for (Long t = 0; t < Trg.Dim(); t++) { // Set pair_lst for (Long elem_idx : near_elem[t]) { pair_lst.PushBack(Pair(elem_idx,Trg[t].rank)); } } } } } } { // Sort and repartition pair_lst Vector> pair_lst_sorted; comm.HyperQuickSort(pair_lst, pair_lst_sorted); Long surf_rank_offset; const Long Nelem = elem_lst.NElem(); comm.Scan(Ptr2ConstItr(&Nelem,1), Ptr2Itr(&surf_rank_offset,1), 1, Comm::CommOp::SUM); surf_rank_offset -= Nelem; comm.PartitionS(pair_lst_sorted, Pair(surf_rank_offset,0)); pair_lst.Swap(pair_lst_sorted); } } template static void BuildNbrListDeprecated(Vector>& pair_lst, const Vector& Xt, const ElemList& elem_lst, const Matrix& surf_nds, Real distance_factor) { using CoordBasis = typename ElemList::CoordBasis; constexpr Integer CoordDim = ElemList::CoordDim(); constexpr Integer ElemDim = ElemList::ElemDim(); const Long Nelem = elem_lst.NElem(); const Long Ntrg = Xt.Dim() / CoordDim; SCTL_ASSERT(Xt.Dim() == Ntrg * CoordDim); Long Nnds, Nsurf_nds; Vector X_surf, X, dX; Integer order_upsample = (Integer)(const_pi() / distance_factor + 0.5); { // Set X, dX const Vector& X_elem_lst = elem_lst.ElemVector(); Vector dX_elem_lst; CoordBasis::Grad(dX_elem_lst, X_elem_lst); Matrix nds_upsample; Vector wts_upsample; TensorProductGaussQuad(nds_upsample, wts_upsample, order_upsample); Nnds = nds_upsample.Dim(1); const auto CoordEvalOp = CoordBasis::SetupEval(nds_upsample); eval_basis(X, X_elem_lst, CoordDim, nds_upsample.Dim(1), CoordEvalOp); eval_basis(dX, dX_elem_lst, CoordDim * ElemDim, nds_upsample.Dim(1), CoordEvalOp); Nsurf_nds = surf_nds.Dim(1); const auto CoordEvalOp_surf = CoordBasis::SetupEval(surf_nds); eval_basis(X_surf, X_elem_lst, CoordDim, Nsurf_nds, CoordEvalOp_surf); } Real d2 = distance_factor * distance_factor; for (Long i = 0; i < Nelem; i++) { std::set near_pts; std::set self_pts; for (Long j = 0; j < Nnds; j++) { Real R2_max = 0; StaticArray X0; for (Integer k = 0; k < CoordDim; k++) { X0[k] = X[(i*Nnds+j)*CoordDim+k]; } for (Integer k0 = 0; k0 < ElemDim; k0++) { Real R2 = 0; for (Integer k1 = 0; k1 < CoordDim; k1++) { Real dX_ = dX[((i*Nnds+j)*CoordDim+k1)*ElemDim+k0]; R2 += dX_*dX_; } R2_max = std::max(R2_max, R2*d2); } for (Long k = 0; k < Ntrg; k++) { Real R2 = 0; for (Integer l = 0; l < CoordDim; l++) { Real dX = Xt[k*CoordDim+l]- X0[l]; R2 += dX * dX; } if (R2 < R2_max) near_pts.insert(k); } } for (Long j = 0; j < Nsurf_nds; j++) { StaticArray X0; for (Integer k = 0; k < CoordDim; k++) { X0[k] = X_surf[(i*Nsurf_nds+j)*CoordDim+k]; } for (Long k = 0; k < Ntrg; k++) { Real R2 = 0; for (Integer l = 0; l < CoordDim; l++) { Real dX = Xt[k*CoordDim+l]- X0[l]; R2 += dX * dX; } if (R2 == 0) self_pts.insert(k); } } for (Long trg_idx : self_pts) { near_pts.erase(trg_idx); } for (Long trg_idx : near_pts) { pair_lst.PushBack(Pair(i,trg_idx)); } } } template static void SetupNearSingular(Matrix& M_near_singular, Vector>& pair_lst, const Vector& Xt_, const Vector& trg_surf, const ElemList& elem_lst, const Kernel& kernel, Integer order_singular, Integer order_direct, Real period_length, const Comm& comm) { static_assert(std::is_same::value); static_assert(std::is_same::value); static_assert(DensityBasis::Dim() == ElemList::ElemDim()); using CoordBasis = typename ElemList::CoordBasis; using CoordEvalOpType = typename CoordBasis::EvalOpType; using DensityEvalOpType = typename DensityBasis::EvalOpType; constexpr Integer CoordDim = ElemList::CoordDim(); constexpr Integer ElemDim = ElemList::ElemDim(); constexpr Integer KDIM0 = Kernel::SrcDim(); constexpr Integer KDIM1 = Kernel::TrgDim(); const Long Nelem = elem_lst.NElem(); BuildNbrList(pair_lst, Xt_, trg_surf, elem_lst, 2.5/order_direct, period_length, comm); const Long Ninterac = pair_lst.Dim(); Vector Xt; { // Set Xt Integer rank = comm.Rank(); Integer np = comm.Size(); Vector splitter_ranks; { // Set splitter_ranks Vector cnt(np); const Long N = Xt_.Dim() / CoordDim; comm.Allgather(Ptr2ConstItr(&N,1), 1, cnt.begin(), 1); scan(splitter_ranks, cnt); } Vector scatter_index, recv_index, recv_cnt(np), recv_dsp(np); { // Set scatter_index, recv_index, recv_cnt, recv_dsp { // Set scatter_index, recv_index Vector> scatter_pair(pair_lst.Dim()); for (Long i = 0; i < pair_lst.Dim(); i++) { scatter_pair[i] = Pair(pair_lst[i].second,i); } omp_par::merge_sort(scatter_pair.begin(), scatter_pair.end()); recv_index.ReInit(scatter_pair.Dim()); scatter_index.ReInit(scatter_pair.Dim()); for (Long i = 0; i < scatter_index.Dim(); i++) { recv_index[i] = scatter_pair[i].first; scatter_index[i] = scatter_pair[i].second; } } for (Integer i = 0; i < np; i++) { recv_dsp[i] = std::lower_bound(recv_index.begin(), recv_index.end(), splitter_ranks[i]) - recv_index.begin(); } for (Integer i = 0; i < np-1; i++) { recv_cnt[i] = recv_dsp[i+1] - recv_dsp[i]; } recv_cnt[np-1] = recv_index.Dim() - recv_dsp[np-1]; } Vector send_index, send_cnt(np), send_dsp(np); { // Set send_index, send_cnt, send_dsp comm.Alltoall(recv_cnt.begin(), 1, send_cnt.begin(), 1); scan(send_dsp, send_cnt); send_index.ReInit(send_cnt[np-1] + send_dsp[np-1]); comm.Alltoallv(recv_index.begin(), recv_cnt.begin(), recv_dsp.begin(), send_index.begin(), send_cnt.begin(), send_dsp.begin()); } Vector Xt_send(send_index.Dim() * CoordDim); for (Long i = 0; i < send_index.Dim(); i++) { // Set Xt_send Long idx = send_index[i] - splitter_ranks[rank]; for (Integer k = 0; k < CoordDim; k++) { Xt_send[i*CoordDim+k] = Xt_[idx*CoordDim+k]; } } Vector Xt_recv(recv_index.Dim() * CoordDim); { // Set Xt_recv for (Long i = 0; i < np; i++) { send_cnt[i] *= CoordDim; send_dsp[i] *= CoordDim; recv_cnt[i] *= CoordDim; recv_dsp[i] *= CoordDim; } comm.Alltoallv(Xt_send.begin(), send_cnt.begin(), send_dsp.begin(), Xt_recv.begin(), recv_cnt.begin(), recv_dsp.begin()); } Xt.ReInit(scatter_index.Dim() * CoordDim); for (Long i = 0; i < scatter_index.Dim(); i++) { // Set Xt Long idx = scatter_index[i]; for (Integer k = 0; k < CoordDim; k++) { Xt[idx*CoordDim+k] = Xt_recv[i*CoordDim+k]; } } } const Vector& X = elem_lst.ElemVector(); Vector dX; CoordBasis::Grad(dX, X); Long elem_rank_offset; { // Set elem_rank_offset comm.Scan(Ptr2ConstItr(&Nelem,1), Ptr2Itr(&elem_rank_offset,1), 1, Comm::CommOp::SUM); elem_rank_offset -= Nelem; } auto& M = M_near_singular; M.ReInit(Ninterac * KDIM0 * DensityBasis::Size(), KDIM1); #pragma omp parallel for schedule(static) for (Long j = 0; j < Ninterac; j++) { // Set M (near-singular) const Long src_idx = pair_lst[j].first - elem_rank_offset; Real adapt = -1.0; Tensor u0; { // Set u0 (project target point to the surface patch in parameter space) ConstIterator Xt_ = Xt.begin() + j * CoordDim; const auto& nodes = CoordBasis::Nodes(); Long min_idx = -1; Real min_R2 = 1e10; for (Long i = 0; i < CoordBasis::Size(); i++) { Real R2 = 0; for (Integer k = 0; k < CoordDim; k++) { Real dX = X[src_idx * CoordDim + k][i] - Xt_[k]; R2 += dX * dX; } if (R2 < min_R2) { min_R2 = R2; min_idx = i; } } SCTL_ASSERT(min_idx >= 0); for (Integer k = 0; k < ElemDim; k++) { u0(k,0) = nodes[k][min_idx]; } for (Integer i = 0; i < 2; i++) { // iterate Matrix X_, dX_; for (Integer k = 0; k < ElemDim; k++) { u0(k,0) = std::min(1.0, u0(k,0)); u0(k,0) = std::max(0.0, u0(k,0)); } const auto eval_op = CoordBasis::SetupEval(Matrix(ElemDim,1,u0.begin(),false)); CoordBasis::Eval(X_, Vector(CoordDim,(Iterator)X.begin()+src_idx*CoordDim,false),eval_op); CoordBasis::Eval(dX_, Vector(CoordDim*ElemDim,dX.begin()+src_idx*CoordDim*ElemDim,false),eval_op); const Tensor x0((Iterator)Xt_); const Tensor x(X_.begin()); const Tensor x_u(dX_.begin()); auto inv = [](const Tensor& M) { Tensor Minv; Real det_inv = 1.0 / (M(0,0)*M(1,1) - M(1,0)*M(0,1)); Minv(0,0) = M(1,1) * det_inv; Minv(0,1) =-M(0,1) * det_inv; Minv(1,0) =-M(1,0) * det_inv; Minv(1,1) = M(0,0) * det_inv; return Minv; }; auto du = inv(x_u.RotateRight()*x_u) * x_u.RotateRight()*(x0-x); u0 = u0 + du; auto x_u_squared = x_u.RotateRight() * x_u; adapt = sctl::sqrt( ((x0-x).RotateRight()*(x0-x))(0,0) / std::max(x_u_squared(0,0),x_u_squared(1,1)) ); } } Matrix quad_nds; Vector quad_wts; DuffyQuad(quad_nds, quad_wts, Vector(ElemDim,u0.begin(),false), order_singular, adapt); const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds); Integer Nnds = quad_wts.Dim(); Vector X_, dX_, Xa_, Xn_; { // Set X_, dX_ const Vector X__(CoordDim, (Iterator)X.begin() + src_idx * CoordDim, false); const Vector dX__(CoordDim * ElemDim, (Iterator)dX.begin() + src_idx * CoordDim * ElemDim, false); eval_basis(X_, X__, CoordDim, Nnds, CoordEvalOp); eval_basis(dX_, dX__, CoordDim * ElemDim, Nnds, CoordEvalOp); } if (CoordDim == 3 && ElemDim == 2) { // Compute Xa_, Xn_ Xa_.ReInit(Nnds); Xn_.ReInit(Nnds*CoordDim); for (Long j = 0; j < Nnds; j++) { StaticArray normal; normal[0] = dX_[j*6+2]*dX_[j*6+5] - dX_[j*6+4]*dX_[j*6+3]; normal[1] = dX_[j*6+4]*dX_[j*6+1] - dX_[j*6+0]*dX_[j*6+5]; normal[2] = dX_[j*6+0]*dX_[j*6+3] - dX_[j*6+2]*dX_[j*6+1]; Xa_[j] = sctl::sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]); Real invXa = 1/Xa_[j]; Xn_[j*3+0] = normal[0] * invXa; Xn_[j*3+1] = normal[1] * invXa; Xn_[j*3+2] = normal[2] * invXa; } } DensityEvalOpType DensityEvalOp; if (std::is_same::value) { DensityEvalOp = CoordEvalOp; } else { DensityEvalOp = DensityBasis::SetupEval(quad_nds); } Matrix M__(Nnds * KDIM0, KDIM1); { // Set kernel matrix M__ const Vector X0_(CoordDim, (Iterator)Xt.begin() + j * CoordDim, false); kernel.template KernelMatrix(M__, X0_, X_, Xn_); } for (Long k0 = 0; k0 < KDIM0; k0++) { for (Long k1 = 0; k1 < KDIM1; k1++) { for (Long l = 0; l < DensityBasis::Size(); l++) { Real M_lk = 0; for (Long n = 0; n < Nnds; n++) { Real quad_wt = Xa_[n] * quad_wts[n]; M_lk += DensityEvalOp[l][n] * quad_wt * M__[n*KDIM0+k0][k1]; } M[(j * KDIM0 + k0) * DensityBasis::Size() + l][k1] = M_lk; } } } } { // Set M (subtract direct) Matrix quad_nds; Vector quad_wts; TensorProductGaussQuad(quad_nds, quad_wts, order_direct); const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds); Integer Nnds = quad_wts.Dim(); Vector X_, dX_, Xa_, Xn_; { // Set X_, dX_ eval_basis(X_, X, CoordDim, Nnds, CoordEvalOp); eval_basis(dX_, dX, CoordDim * ElemDim, Nnds, CoordEvalOp); } if (CoordDim == 3 && ElemDim == 2) { // Compute Xa_, Xn_ Long N = Nelem*Nnds; Xa_.ReInit(N); Xn_.ReInit(N*CoordDim); for (Long j = 0; j < N; j++) { StaticArray normal; normal[0] = dX_[j*6+2]*dX_[j*6+5] - dX_[j*6+4]*dX_[j*6+3]; normal[1] = dX_[j*6+4]*dX_[j*6+1] - dX_[j*6+0]*dX_[j*6+5]; normal[2] = dX_[j*6+0]*dX_[j*6+3] - dX_[j*6+2]*dX_[j*6+1]; Xa_[j] = sctl::sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]); Real invXa = 1/Xa_[j]; Xn_[j*3+0] = normal[0] * invXa; Xn_[j*3+1] = normal[1] * invXa; Xn_[j*3+2] = normal[2] * invXa; } } DensityEvalOpType DensityEvalOp; if (std::is_same::value) { DensityEvalOp = CoordEvalOp; } else { DensityEvalOp = DensityBasis::SetupEval(quad_nds); } #pragma omp parallel for schedule(static) for (Long j = 0; j < Ninterac; j++) { // Subtract direct contribution const Long src_idx = pair_lst[j].first - elem_rank_offset; Matrix M__(Nnds * KDIM0, KDIM1); { // Set kernel matrix M__ const Vector X0_(CoordDim, (Iterator)Xt.begin() + j * CoordDim, false); Vector X__(Nnds * CoordDim, X_.begin() + src_idx * Nnds * CoordDim, false); Vector Xn__(Nnds * CoordDim, Xn_.begin() + src_idx * Nnds * CoordDim, false); kernel.template KernelMatrix(M__, X0_, X__, Xn__); } for (Long k0 = 0; k0 < KDIM0; k0++) { for (Long k1 = 0; k1 < KDIM1; k1++) { for (Long l = 0; l < DensityBasis::Size(); l++) { Real M_lk = 0; for (Long n = 0; n < Nnds; n++) { Real quad_wt = Xa_[src_idx * Nnds + n] * quad_wts[n]; M_lk += DensityEvalOp[l][n] * quad_wt * M__[n*KDIM0+k0][k1]; } M[(j * KDIM0 + k0) * DensityBasis::Size() + l][k1] -= M_lk; } } } } } } template static void EvalNearSingular(Vector& U, const Vector& density, const Matrix& M, const Vector>& pair_lst, Long Nelem_, Long Ntrg_, Integer KDIM0_, Integer KDIM1_, const Comm& comm) { const Long Ninterac = pair_lst.Dim(); const Integer dof = density.Dim() / Nelem_ / KDIM0_; SCTL_ASSERT(density.Dim() == Nelem_ * dof * KDIM0_); Long elem_rank_offset; { // Set elem_rank_offset comm.Scan(Ptr2ConstItr(&Nelem_,1), Ptr2Itr(&elem_rank_offset,1), 1, Comm::CommOp::SUM); elem_rank_offset -= Nelem_; } Vector U_loc(Ninterac*dof*KDIM1_); for (Long j = 0; j < Ninterac; j++) { const Long src_idx = pair_lst[j].first - elem_rank_offset; const Matrix M_(KDIM0_ * DensityBasis::Size(), KDIM1_, (Iterator)M[j * KDIM0_ * DensityBasis::Size()], false); Matrix U_(dof, KDIM1_, U_loc.begin() + j*dof*KDIM1_, false); Matrix F_(dof, KDIM0_ * DensityBasis::Size()); for (Long i = 0; i < dof; i++) { for (Long k = 0; k < KDIM0_; k++) { for (Long l = 0; l < DensityBasis::Size(); l++) { F_[i][k * DensityBasis::Size() + l] = density[(src_idx * dof + i) * KDIM0_ + k][l]; } } } Matrix::GEMM(U_, F_, M_); } if (U.Dim() != Ntrg_ * dof * KDIM1_) { U.ReInit(Ntrg_ * dof * KDIM1_); U = 0; } { // Set U Integer rank = comm.Rank(); Integer np = comm.Size(); Vector splitter_ranks; { // Set splitter_ranks Vector cnt(np); comm.Allgather(Ptr2ConstItr(&Ntrg_,1), 1, cnt.begin(), 1); scan(splitter_ranks, cnt); } Vector scatter_index, send_index, send_cnt(np), send_dsp(np); { // Set scatter_index, send_index, send_cnt, send_dsp { // Set scatter_index, send_index Vector> scatter_pair(pair_lst.Dim()); for (Long i = 0; i < pair_lst.Dim(); i++) { scatter_pair[i] = Pair(pair_lst[i].second,i); } omp_par::merge_sort(scatter_pair.begin(), scatter_pair.end()); send_index.ReInit(scatter_pair.Dim()); scatter_index.ReInit(scatter_pair.Dim()); for (Long i = 0; i < scatter_index.Dim(); i++) { send_index[i] = scatter_pair[i].first; scatter_index[i] = scatter_pair[i].second; } } for (Integer i = 0; i < np; i++) { send_dsp[i] = std::lower_bound(send_index.begin(), send_index.end(), splitter_ranks[i]) - send_index.begin(); } for (Integer i = 0; i < np-1; i++) { send_cnt[i] = send_dsp[i+1] - send_dsp[i]; } send_cnt[np-1] = send_index.Dim() - send_dsp[np-1]; } Vector recv_index, recv_cnt(np), recv_dsp(np); { // Set recv_index, recv_cnt, recv_dsp comm.Alltoall(send_cnt.begin(), 1, recv_cnt.begin(), 1); scan(recv_dsp, recv_cnt); recv_index.ReInit(recv_cnt[np-1] + recv_dsp[np-1]); comm.Alltoallv(send_index.begin(), send_cnt.begin(), send_dsp.begin(), recv_index.begin(), recv_cnt.begin(), recv_dsp.begin()); } Vector U_send(scatter_index.Dim() * dof * KDIM1_); for (Long i = 0; i < scatter_index.Dim(); i++) { Long idx = scatter_index[i]*dof*KDIM1_; for (Long k = 0; k < dof * KDIM1_; k++) { U_send[i*dof*KDIM1_ + k] = U_loc[idx + k]; } } Vector U_recv(recv_index.Dim() * dof * KDIM1_); { // Set U_recv for (Long i = 0; i < np; i++) { send_cnt[i] *= dof * KDIM1_; send_dsp[i] *= dof * KDIM1_; recv_cnt[i] *= dof * KDIM1_; recv_dsp[i] *= dof * KDIM1_; } comm.Alltoallv(U_send.begin(), send_cnt.begin(), send_dsp.begin(), U_recv.begin(), recv_cnt.begin(), recv_dsp.begin()); } for (Long i = 0; i < recv_index.Dim(); i++) { // Set U Long idx = (recv_index[i] - splitter_ranks[rank]) * dof * KDIM1_; for (Integer k = 0; k < dof * KDIM1_; k++) { U[idx + k] += U_recv[i*dof*KDIM1_ + k]; } } } } template static void Direct(Vector& U, const Vector& Xt, const ElemList& elem_lst, const Vector& density, const Kernel& kernel, Integer order_direct, const Comm& comm) { using CoordBasis = typename ElemList::CoordBasis; using CoordEvalOpType = typename CoordBasis::EvalOpType; using DensityEvalOpType = typename DensityBasis::EvalOpType; constexpr Integer CoordDim = ElemList::CoordDim(); constexpr Integer ElemDim = ElemList::ElemDim(); constexpr Integer KDIM0 = Kernel::SrcDim(); constexpr Integer KDIM1 = Kernel::TrgDim(); const Long Nelem = elem_lst.NElem(); const Integer dof = density.Dim() / Nelem / KDIM0; SCTL_ASSERT(density.Dim() == Nelem * dof * KDIM0); Matrix quad_nds; Vector quad_wts; TensorProductGaussQuad(quad_nds, quad_wts, order_direct); const CoordEvalOpType CoordEvalOp = CoordBasis::SetupEval(quad_nds); Integer Nnds = quad_wts.Dim(); const Vector& X = elem_lst.ElemVector(); Vector dX; CoordBasis::Grad(dX, X); Vector X_, dX_, Xa_, Xn_; eval_basis(X_, X, CoordDim, Nnds, CoordEvalOp); eval_basis(dX_, dX, CoordDim*ElemDim, Nnds, CoordEvalOp); if (CoordDim == 3 && ElemDim == 2) { // Compute Xa_, Xn_ Long N = Nelem*Nnds; Xa_.ReInit(N); Xn_.ReInit(N*CoordDim); for (Long j = 0; j < N; j++) { StaticArray normal; normal[0] = dX_[j*6+2]*dX_[j*6+5] - dX_[j*6+4]*dX_[j*6+3]; normal[1] = dX_[j*6+4]*dX_[j*6+1] - dX_[j*6+0]*dX_[j*6+5]; normal[2] = dX_[j*6+0]*dX_[j*6+3] - dX_[j*6+2]*dX_[j*6+1]; Xa_[j] = sctl::sqrt(normal[0]*normal[0]+normal[1]*normal[1]+normal[2]*normal[2]); Real invXa = 1/Xa_[j]; Xn_[j*3+0] = normal[0] * invXa; Xn_[j*3+1] = normal[1] * invXa; Xn_[j*3+2] = normal[2] * invXa; } } Vector Fa_; { // Set Fa_ Vector F_; if (std::is_same::value) { eval_basis(F_, density, dof * KDIM0, Nnds, CoordEvalOp); } else { const DensityEvalOpType EvalOp = DensityBasis::SetupEval(quad_nds); eval_basis(F_, density, dof * KDIM0, Nnds, EvalOp); } Fa_.ReInit(F_.Dim()); const Integer DensityDOF = dof * KDIM0; SCTL_ASSERT(F_.Dim() == Nelem * Nnds * DensityDOF); for (Long j = 0; j < Nelem; j++) { for (Integer k = 0; k < Nnds; k++) { Long idx = j * Nnds + k; Real quad_wt = Xa_[idx] * quad_wts[k]; for (Integer l = 0; l < DensityDOF; l++) { Fa_[idx * DensityDOF + l] = F_[idx * DensityDOF + l] * quad_wt; } } } } { // Evaluate potential const Long Ntrg = Xt.Dim() / CoordDim; SCTL_ASSERT(Xt.Dim() == Ntrg * CoordDim); if (U.Dim() != Ntrg * dof * KDIM1) { U.ReInit(Ntrg * dof * KDIM1); U = 0; } ParticleFMM::Eval(U, Xt, X_, Xn_, Fa_, kernel, comm); } } public: template void Setup(const ElemList& elem_lst, const Vector& Xt, const Kernel& kernel, Integer order_singular, Integer order_direct, Real period_length, const Comm& comm) { Xt_.ReInit(0); M_singular.ReInit(0,0); M_near_singular.ReInit(0,0); pair_lst.ReInit(0); order_direct_ = order_direct; period_length_ = period_length; comm_ = comm; Profile::Tic("Setup", &comm_); static_assert(std::is_same::value); static_assert(std::is_same::value); static_assert(DensityBasis::Dim() == ElemList::ElemDim()); Xt_ = Xt; M_singular.ReInit(0,0); Profile::Tic("SetupNearSingular", &comm_); SetupNearSingular(M_near_singular, pair_lst, Xt_, Vector(), elem_lst, kernel, order_singular, order_direct_, period_length_, comm_); Profile::Toc(); Profile::Toc(); } template void Setup(const ElemList& elem_lst, const Kernel& kernel, Integer order_singular, Integer order_direct, Real period_length, const Comm& comm, Real Rqbx = 0) { Xt_.ReInit(0); M_singular.ReInit(0,0); M_near_singular.ReInit(0,0); pair_lst.ReInit(0); order_direct_ = order_direct; period_length_ = period_length; comm_ = comm; Profile::Tic("Setup", &comm_); static_assert(std::is_same::value); static_assert(std::is_same::value); static_assert(std::is_same::value); static_assert(PotentialBasis::Dim() == ElemList::ElemDim()); static_assert(DensityBasis::Dim() == ElemList::ElemDim()); Vector trg_surf; { // Set Xt_ using CoordBasis = typename ElemList::CoordBasis; Matrix trg_nds = PotentialBasis::Nodes(); auto Meval = CoordBasis::SetupEval(trg_nds); eval_basis(Xt_, elem_lst.ElemVector(), ElemList::CoordDim(), trg_nds.Dim(1), Meval); { // Set trg_surf const Long Nelem = elem_lst.NElem(); const Long Nnds = trg_nds.Dim(1); Long elem_offset; { // Set elem_offset comm.Scan(Ptr2ConstItr(&Nelem,1), Ptr2Itr(&elem_offset,1), 1, Comm::CommOp::SUM); elem_offset -= Nelem; } trg_surf.ReInit(elem_lst.NElem() * trg_nds.Dim(1)); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnds; j++) { trg_surf[i*Nnds+j] = elem_offset + i; } } } } Profile::Tic("SetupSingular", &comm_); SetupSingular(M_singular, PotentialBasis::Nodes(), elem_lst, kernel, order_singular, order_direct_, Rqbx); Profile::Toc(); Profile::Tic("SetupNearSingular", &comm_); SetupNearSingular(M_near_singular, pair_lst, Xt_, trg_surf, elem_lst, kernel, order_singular, order_direct_, period_length_, comm_); Profile::Toc(); Profile::Toc(); } template void Eval(Vector& U, const ElemList& elements, const Vector& F, const Kernel& kernel) const { Profile::Tic("Eval", &comm_); Matrix U_singular; Vector U_direct, U_near_sing; Profile::Tic("EvalDirect", &comm_); Direct(U_direct, Xt_, elements, F, kernel, order_direct_, comm_); Profile::Toc(); Profile::Tic("EvalSingular", &comm_); EvalSingular(U_singular, F, M_singular, kernel.SrcDim(), kernel.TrgDim()); Profile::Toc(); Profile::Tic("EvalNearSingular", &comm_); EvalNearSingular(U_near_sing, F, M_near_singular, pair_lst, elements.NElem(), Xt_.Dim() / ElemList::CoordDim(), kernel.SrcDim(), kernel.TrgDim(), comm_); SCTL_ASSERT(U_near_sing.Dim() == U_direct.Dim()); Profile::Toc(); const Long dof = U_direct.Dim() / (elements.NElem() * PotentialBasis::Size() * kernel.TrgDim()); SCTL_ASSERT(U_direct .Dim() == elements.NElem() * PotentialBasis::Size() * dof * kernel.TrgDim()); SCTL_ASSERT(U_near_sing.Dim() == elements.NElem() * PotentialBasis::Size() * dof * kernel.TrgDim()); if (U.Dim() != elements.NElem() * dof * kernel.TrgDim()) { U.ReInit(elements.NElem() * dof * kernel.TrgDim()); } for (int i = 0; i < elements.NElem(); i++) { for (int j = 0; j < PotentialBasis::Size(); j++) { for (int k = 0; k < dof*kernel.TrgDim(); k++) { Real& U_ = U[i*dof*kernel.TrgDim()+k][j]; U_ = 0; U_ += U_direct [(i*PotentialBasis::Size()+j)*dof*kernel.TrgDim()+k]; U_ += U_near_sing[(i*PotentialBasis::Size()+j)*dof*kernel.TrgDim()+k]; U_ *= kernel.template ScaleFactor(); } } } if (U_singular.Dim(1)) { SCTL_ASSERT(U_singular.Dim(0) == elements.NElem() * dof * kernel.TrgDim()); SCTL_ASSERT(U_singular.Dim(1) == PotentialBasis::Size()); for (int i = 0; i < elements.NElem(); i++) { for (int j = 0; j < PotentialBasis::Size(); j++) { for (int k = 0; k < dof*kernel.TrgDim(); k++) { U[i*dof*kernel.TrgDim()+k][j] += U_singular[i*dof*kernel.TrgDim()+k][j] * kernel.template ScaleFactor(); } } } } Profile::Toc(); } template void Eval(Vector& U, const ElemList& elements, const Vector& F, const Kernel& kernel) const { Profile::Tic("Eval", &comm_); Matrix U_singular; Vector U_direct, U_near_sing; Profile::Tic("EvalDirect", &comm_); Direct(U_direct, Xt_, elements, F, kernel, order_direct_, comm_); Profile::Toc(); Profile::Tic("EvalSingular", &comm_); EvalSingular(U_singular, F, M_singular, kernel.SrcDim(), kernel.TrgDim()); Profile::Toc(); Profile::Tic("EvalNearSingular", &comm_); EvalNearSingular(U_near_sing, F, M_near_singular, pair_lst, elements.NElem(), Xt_.Dim() / ElemList::CoordDim(), kernel.SrcDim(), kernel.TrgDim(), comm_); SCTL_ASSERT(U_near_sing.Dim() == U_direct.Dim()); Profile::Toc(); Long Nt = Xt_.Dim() / ElemList::CoordDim(); const Long dof = U_direct.Dim() / (Nt * kernel.TrgDim()); SCTL_ASSERT(U_direct.Dim() == Nt * dof * kernel.TrgDim()); if (U.Dim() != U_direct.Dim()) { U.ReInit(U_direct.Dim()); } for (int i = 0; i < U.Dim(); i++) { U[i] = (U_direct[i] + U_near_sing[i]) * kernel.template ScaleFactor(); } if (U_singular.Dim(1)) { SCTL_ASSERT(U_singular.Dim(0) == elements.NElem() * dof * kernel.TrgDim()); const Long Nnodes = U_singular.Dim(1); for (int i = 0; i < elements.NElem(); i++) { for (int j = 0; j < Nnodes; j++) { for (int k = 0; k < dof*kernel.TrgDim(); k++) { Real& U_ = U[(i*Nnodes+j)*dof*kernel.TrgDim()+k]; U_ += U_singular[i*dof*kernel.TrgDim()+k][j] * kernel.template ScaleFactor(); } } } } Profile::Toc(); } template static void test(Integer order_singular = 10, Integer order_direct = 5, const Comm& comm = Comm::World()) { constexpr Integer COORD_DIM = 3; constexpr Integer ELEM_DIM = COORD_DIM-1; using ElemList = ElemList>; using DensityBasis = Basis; using PotentialBasis = Basis; int np = comm.Size(); int rank = comm.Rank(); auto build_torus = [rank,np](ElemList& elements, long Nt, long Np, Real Rmajor, Real Rminor){ auto nodes = ElemList::CoordBasis::Nodes(); auto torus = [](Real theta, Real phi, Real Rmajor, Real Rminor) { Real R = Rmajor + Rminor * cos(phi); Real X = R * cos(theta); Real Y = R * sin(theta); Real Z = Rminor * sin(phi); return std::make_tuple(X,Y,Z); }; long start = Nt*Np*(rank+0)/np; long end = Nt*Np*(rank+1)/np; elements.ReInit(end - start); for (long ii = start; ii < end; ii++) { long i = ii / Np; long j = ii % Np; for (int k = 0; k < ElemList::CoordBasis::Size(); k++) { Real X, Y, Z; Real theta = 2 * const_pi() * (i + nodes[0][k]) / Nt; Real phi = 2 * const_pi() * (j + nodes[1][k]) / Np; std::tie(X,Y,Z) = torus(theta, phi, Rmajor, Rminor); elements(ii-start,0)[k] = X; elements(ii-start,1)[k] = Y; elements(ii-start,2)[k] = Z; } } }; ElemList elements_src, elements_trg; build_torus(elements_src, 28, 16, 2, 1.0); build_torus(elements_trg, 29, 17, 2, 0.99); Vector Xt; Vector U_onsurf, U_offsurf; Vector density_sl, density_dl; { // Set Xt, elements_src, elements_trg, density_sl, density_dl, U Real X0[COORD_DIM] = {3,2,1}; std::function potential = [X0](Real* U, Real* X, Real* Xn) { Real dX[COORD_DIM] = {X[0]-X0[0],X[1]-X0[1],X[2]-X0[2]}; Real Rinv = 1/sqrt(dX[0]*dX[0]+dX[1]*dX[1]+dX[2]*dX[2]); U[0] = Rinv; }; std::function potential_normal_derivative = [X0](Real* U, Real* X, Real* Xn) { Real dX[COORD_DIM] = {X[0]-X0[0],X[1]-X0[1],X[2]-X0[2]}; Real Rinv = 1/sqrt(dX[0]*dX[0]+dX[1]*dX[1]+dX[2]*dX[2]); Real RdotN = dX[0]*Xn[0]+dX[1]*Xn[1]+dX[2]*Xn[2]; U[0] = -RdotN * Rinv*Rinv*Rinv; }; DiscretizeSurfaceFn(density_sl, elements_src, potential_normal_derivative); DiscretizeSurfaceFn(density_dl, elements_src, potential); DiscretizeSurfaceFn(U_onsurf , elements_src, potential); DiscretizeSurfaceFn(U_offsurf , elements_trg, potential); for (long i = 0; i < elements_trg.NElem(); i++) { // Set Xt for (long j = 0; j < PotentialBasis::Size(); j++) { for (int k = 0; k < COORD_DIM; k++) { Xt.PushBack(elements_trg(i,k)[j]); } } } } GenericKernel Laplace_DxU; GenericKernel Laplace_FxU; Profile::Enable(true); if (1) { // Greeen's identity test (Laplace, on-surface) Profile::Tic("OnSurface", &comm); Quadrature quadrature_DxU, quadrature_FxU; quadrature_FxU.Setup(elements_src, Laplace_FxU, order_singular, order_direct, -1.0, comm); quadrature_DxU.Setup(elements_src, Laplace_DxU, order_singular, order_direct, -1.0, comm); Vector U_sl, U_dl; quadrature_FxU.Eval(U_sl, elements_src, density_sl, Laplace_FxU); quadrature_DxU.Eval(U_dl, elements_src, density_dl, Laplace_DxU); Profile::Toc(); Real max_err = 0; Vector err(U_onsurf.Dim()); for (long i = 0; i < U_sl.Dim(); i++) { for (long j = 0; j < PotentialBasis::Size(); j++) { err[i][j] = 0.5*U_onsurf[i][j] - (U_sl[i][j] + U_dl[i][j]); max_err = std::max(max_err, fabs(err[i][j])); } } { // Print error Real glb_err; comm.Allreduce(Ptr2ConstItr(&max_err,1), Ptr2Itr(&glb_err,1), 1, Comm::CommOp::MAX); if (!comm.Rank()) std::cout<<"Error = "< quadrature_DxU, quadrature_FxU; quadrature_FxU.Setup(elements_src, Xt, Laplace_FxU, order_singular, order_direct, -1.0, comm); quadrature_DxU.Setup(elements_src, Xt, Laplace_DxU, order_singular, order_direct, -1.0, comm); Vector U_sl, U_dl; quadrature_FxU.Eval(U_sl, elements_src, density_sl, Laplace_FxU); quadrature_DxU.Eval(U_dl, elements_src, density_dl, Laplace_DxU); Profile::Toc(); Real max_err = 0; Vector err(elements_trg.NElem()); for (long i = 0; i < elements_trg.NElem(); i++) { for (long j = 0; j < PotentialBasis::Size(); j++) { err[i][j] = U_offsurf[i][j] - (U_sl[i*PotentialBasis::Size()+j] + U_dl[i*PotentialBasis::Size()+j]); max_err = std::max(max_err, fabs(err[i][j])); } } { // Print error Real glb_err; comm.Allreduce(Ptr2ConstItr(&max_err,1), Ptr2Itr(&glb_err,1), 1, Comm::CommOp::MAX); if (!comm.Rank()) std::cout<<"Error = "<>; using DensityBasis = Basis; using PotentialBasis = Basis; int np = comm.Size(); int rank = comm.Rank(); auto build_sphere = [rank,np](ElemList& elements, Real X, Real Y, Real Z, Real R){ auto nodes = ElemList::CoordBasis::Nodes(); long start = 2*COORD_DIM*(rank+0)/np; long end = 2*COORD_DIM*(rank+1)/np; elements.ReInit(end - start); for (long ii = start; ii < end; ii++) { long i = ii / 2; long j = ii % 2; for (int k = 0; k < ElemList::CoordBasis::Size(); k++) { Real coord[COORD_DIM]; coord[(i+0)%COORD_DIM] = (j ? -1.0 : 1.0); coord[(i+1)%COORD_DIM] = 2.0 * nodes[j?1:0][k] - 1.0; coord[(i+2)%COORD_DIM] = 2.0 * nodes[j?0:1][k] - 1.0; Real R0 = sqrt(coord[0]*coord[0] + coord[1]*coord[1] + coord[2]*coord[2]); elements(ii-start,0)[k] = X + R * coord[0] / R0; elements(ii-start,1)[k] = Y + R * coord[1] / R0; elements(ii-start,2)[k] = Z + R * coord[2] / R0; } } }; ElemList elements; build_sphere(elements, 0.0, 0.0, 0.0, 1.00); Vector density_sl; { // Set density_sl std::function sigma = [](Real* U, Real* X, Real* Xn) { Real R = sqrt(X[0]*X[0]+X[1]*X[1]+X[2]*X[2]); Real sinp = sqrt(X[1]*X[1] + X[2]*X[2]) / R; Real cosp = -X[0] / R; U[0] = -1.5; U[1] = 0; U[2] = 0; }; DiscretizeSurfaceFn(density_sl, elements, sigma); } GenericKernel Stokes_DxU; GenericKernel Stokes_FxU; Profile::Enable(true); if (1) { Vector U; Quadrature quadrature_FxU; quadrature_FxU.Setup(elements, Stokes_FxU, order_singular, order_direct, -1.0, comm); quadrature_FxU.Eval(U, elements, density_sl, Stokes_FxU); { // Write VTK output VTUData vtu; vtu.AddElems(elements, U, ORDER); vtu.WriteVTK("U", comm); } { // Write VTK output VTUData vtu; vtu.AddElems(elements, density_sl, ORDER); vtu.WriteVTK("sigma", comm); } } Profile::print(&comm); } private: static void scan(Vector& dsp, const Vector& cnt) { dsp.ReInit(cnt.Dim()); if (cnt.Dim()) dsp[0] = 0; omp_par::scan(cnt.begin(), dsp.begin(), cnt.Dim()); } template static void eval_basis(Vector& value, const Vector X, Integer dof, Integer Nnds, const typename Basis::EvalOpType& EvalOp) { Long Nelem = X.Dim() / dof; SCTL_ASSERT(X.Dim() == Nelem * dof); value.ReInit(Nelem*Nnds*dof); Matrix X_(Nelem*dof, Nnds, value.begin(),false); Basis::Eval(X_, X, EvalOp); for (Long j = 0; j < Nelem; j++) { // Rearrange data Matrix X(Nnds, dof, X_[j*dof], false); X = Matrix(dof, Nnds, X_[j*dof], false).Transpose(); } } template static void DiscretizeSurfaceFn(Vector& U, const ElemList& elements, std::function fn) { using CoordBasis = typename ElemList::CoordBasis; const long Nelem = elements.NElem(); U.ReInit(Nelem * FnDim); Matrix X, X_grad; { // Set X, X_grad Vector coord = elements.ElemVector(); Vector coord_grad; CoordBasis::Grad(coord_grad, coord); const auto Meval = CoordBasis::SetupEval(FnBasis::Nodes()); CoordBasis::Eval(X, coord, Meval); CoordBasis::Eval(X_grad, coord_grad, Meval); } for (long i = 0; i < Nelem; i++) { for (long j = 0; j < FnBasis::Size(); j++) { Real X_[CoordDim], Xn[CoordDim], U_[FnDim]; for (long k = 0; k < CoordDim; k++) { X_[k] = X[i*CoordDim+k][j]; } { // Set Xn Real Xu[CoordDim], Xv[CoordDim]; for (long k = 0; k < CoordDim; k++) { Xu[k] = X_grad[(i*CoordDim+k)*2+0][j]; Xv[k] = X_grad[(i*CoordDim+k)*2+1][j]; } Real dA = 0; for (long k = 0; k < CoordDim; k++) { Xn[k] = Xu[(k+1)%CoordDim] * Xv[(k+2)%CoordDim]; Xn[k] -= Xv[(k+1)%CoordDim] * Xu[(k+2)%CoordDim]; dA += Xn[k] * Xn[k]; } dA = sqrt(dA); for (long k = 0; k < CoordDim; k++) { Xn[k] /= dA; } } fn(U_, X_, Xn); for (long k = 0; k < FnDim; k++) { U[i*FnDim+k][j] = U_[k]; } } } } Vector Xt_; Matrix M_singular; Matrix M_near_singular; Vector> pair_lst; Integer order_direct_; Real period_length_; Comm comm_; }; template class Stellarator { private: static constexpr Integer order_singular = 25; static constexpr Integer order_direct = 35; static constexpr Integer COORD_DIM = 3; static constexpr Integer ELEM_DIM = COORD_DIM-1; using ElemBasis = Basis; using ElemLst = ElemList; struct Laplace3D_dUxF { template static constexpr ValueType ScaleFactor() { return 1 / (4 * const_pi()); } template static void Eval(ValueType (&u)[3][1], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) { ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2]; ValueType rinv = (r2>1e-16 ? 1/sqrt(r2) : 0); ValueType rinv3 = rinv * rinv * rinv; u[0][0] = -r[0] * rinv3; u[1][0] = -r[1] * rinv3; u[2][0] = -r[2] * rinv3; } }; struct BiotSavart3D { template static constexpr ValueType ScaleFactor() { return 1 / (4 * const_pi()); } template static void Eval(ValueType (&u)[3][3], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) { ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2]; ValueType rinv = (r2>1e-16 ? 1/sqrt(r2) : 0); ValueType rinv3 = rinv * rinv * rinv; u[0][0] = (0) * rinv3; u[1][0] = r[2] * rinv3; u[2][0] = -r[1] * rinv3; u[0][1] = -r[2] * rinv3; u[1][1] = (0) * rinv3; u[2][1] = r[0] * rinv3; u[0][2] = r[1] * rinv3; u[1][2] = -r[0] * rinv3; u[2][2] = (0) * rinv3; } }; struct BiotSavartGrad3D { template static constexpr ValueType ScaleFactor() { return 1 / (4 * const_pi()); } template static void Eval(ValueType (&u)[3][9], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) { ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2]; ValueType rinv = (r2>1e-16 ? 1/sqrt(r2) : 0); ValueType rinv2 = rinv * rinv; ValueType rinv3 = rinv2 * rinv; ValueType rinv5 = rinv2 * rinv3; u[0][0] = 0; u[1][0] = - 3 * r[2] * r[0] * rinv5; u[2][0] = 3 * r[1] * r[0] * rinv5; u[0][1] = 0; u[1][1] = - 3 * r[2] * r[1] * rinv5; u[2][1] = -(1) * rinv3 + 3 * r[1] * r[1] * rinv5; u[0][2] = 0; u[1][2] = (1) * rinv3 - 3 * r[2] * r[2] * rinv5; u[2][2] = 3 * r[1] * r[2] * rinv5; u[0][3] = 3 * r[2] * r[0] * rinv5; u[1][3] = 0; u[2][3] = (1) * rinv3 - 3 * r[0] * r[0] * rinv5; u[0][4] = 3 * r[2] * r[1] * rinv5; u[1][4] = 0; u[2][4] = - 3 * r[0] * r[1] * rinv5; u[0][5] = -(1) * rinv3 + 3 * r[2] * r[2] * rinv5; u[1][5] = 0; u[2][5] = - 3 * r[0] * r[2] * rinv5; u[0][6] = - 3 * r[1] * r[0] * rinv5; u[1][6] = -(1) * rinv3 + 3 * r[0] * r[0] * rinv5; u[2][6] = 0; u[0][7] = (1) * rinv3 - 3 * r[1] * r[1] * rinv5; u[1][7] = 3 * r[0] * r[1] * rinv5; u[2][7] = 0; u[0][8] = - 3 * r[1] * r[2] * rinv5; u[1][8] = 3 * r[0] * r[2] * rinv5; u[2][8] = 0; } }; struct Laplace3D_dUxD { template static constexpr ValueType ScaleFactor() { return 1 / (4 * const_pi()); } template static void Eval(ValueType (&u)[3][3], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) { ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2]; ValueType rinv = (r2>1e-16 ? 1/sqrt(r2) : 0); ValueType rdotn = r[0]*n[0] + r[1]*n[1] + r[2]*n[2]; ValueType rinv2 = rinv * rinv; ValueType rinv3 = rinv * rinv2; ValueType rinv5 = rinv3 * rinv2; u[0][0] = -1 * rinv3 + 3 * r[0] * r[0] * rinv5; u[0][1] = -0 * rinv3 + 3 * r[0] * r[1] * rinv5; u[0][2] = -0 * rinv3 + 3 * r[0] * r[2] * rinv5; u[1][0] = -0 * rinv3 + 3 * r[1] * r[0] * rinv5; u[1][1] = -1 * rinv3 + 3 * r[1] * r[1] * rinv5; u[1][2] = -0 * rinv3 + 3 * r[1] * r[2] * rinv5; u[2][0] = -0 * rinv3 + 3 * r[2] * r[0] * rinv5; u[2][1] = -0 * rinv3 + 3 * r[2] * r[1] * rinv5; u[2][2] = -1 * rinv3 + 3 * r[2] * r[2] * rinv5; } }; struct Laplace3D_DxdU { template static constexpr ValueType ScaleFactor() { return 1 / (4 * const_pi()); } template static void Eval(ValueType (&u)[1][3], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) { ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2]; ValueType rinv = (r2>1e-16 ? 1/sqrt(r2) : 0); ValueType rdotn = r[0]*n[0] + r[1]*n[1] + r[2]*n[2]; ValueType rinv2 = rinv * rinv; ValueType rinv3 = rinv * rinv2; ValueType rinv5 = rinv3 * rinv2; u[0][0] = -n[0] * rinv3 + 3*rdotn * r[0] * rinv5; u[0][1] = -n[1] * rinv3 + 3*rdotn * r[1] * rinv5; u[0][2] = -n[2] * rinv3 + 3*rdotn * r[2] * rinv5; } }; struct Laplace3D_Fxd2U { template static constexpr ValueType ScaleFactor() { return 1 / (4 * const_pi()); } template static void Eval(ValueType (&u)[1][9], const ValueType (&r)[3], const ValueType (&n)[3], void* ctx_ptr) { ValueType r2 = r[0]*r[0]+r[1]*r[1]+r[2]*r[2]; ValueType rinv = (r2>1e-16 ? 1/sqrt(r2) : 0); ValueType rinv2 = rinv * rinv; ValueType rinv3 = rinv * rinv2; ValueType rinv5 = rinv3 * rinv2; u[0][0+3*0] = -1 * rinv3 + 3 * r[0] * r[0] * rinv5; u[0][1+3*0] = -0 * rinv3 + 3 * r[0] * r[1] * rinv5; u[0][2+3*0] = -0 * rinv3 + 3 * r[0] * r[2] * rinv5; u[0][0+3*1] = -0 * rinv3 + 3 * r[1] * r[0] * rinv5; u[0][1+3*1] = -1 * rinv3 + 3 * r[1] * r[1] * rinv5; u[0][2+3*1] = -0 * rinv3 + 3 * r[1] * r[2] * rinv5; u[0][0+3*2] = -0 * rinv3 + 3 * r[2] * r[0] * rinv5; u[0][1+3*2] = -0 * rinv3 + 3 * r[2] * r[1] * rinv5; u[0][2+3*2] = -1 * rinv3 + 3 * r[2] * r[2] * rinv5; } }; static Real max_norm(const sctl::Vector& x) { Real err = 0; for (const auto& a : x) err = std::max(err, sctl::fabs(a)); return err; } public: static Vector compute_dot_prod(const Vector& A, const Vector& B) { const Long Nelem = A.Dim() / COORD_DIM; const Long Nnodes = ElemBasis::Size(); SCTL_ASSERT(A.Dim() == Nelem * COORD_DIM); SCTL_ASSERT(B.Dim() == Nelem * COORD_DIM); Vector AdotB(Nelem); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { Real a_dot_b = 0; a_dot_b += A[i*COORD_DIM+0][j]*B[i*COORD_DIM+0][j]; a_dot_b += A[i*COORD_DIM+1][j]*B[i*COORD_DIM+1][j]; a_dot_b += A[i*COORD_DIM+2][j]*B[i*COORD_DIM+2][j]; AdotB[i][j] = a_dot_b; } } return AdotB; } static Real compute_inner_prod(const Vector& area_elem, const Vector& A, const Vector& B) { const auto& quad_wts = ElemBasis::QuadWts(); const Long Nnodes = ElemBasis::Size(); const Long Nelem = area_elem.Dim(); const Long dof = B.Dim() / Nelem; Real sum = 0; for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { Real AdotB = 0; for (Long k = 0; k < dof; k++) { AdotB += A[i*dof+k][j] * B[i*dof+k][j]; } sum += AdotB * area_elem[i][j] * quad_wts[j]; } } return sum; } static void compute_harmonic_vector_potentials(Vector& Jt, Vector& Jp, const Stellarator& S) { Comm comm = Comm::World(); Real gmres_tol = 1e-8; Long max_iter = 100; auto cheb2grid = [] (const Vector& X, Long Mt, Long Mp, Long Nt, Long Np) { const Long dof = X.Dim() / (Mt * Mp); SCTL_ASSERT(X.Dim() == Mt * Mp *dof); Vector Xf(dof*Nt*Np); Xf = 0; const Long Nnodes = ElemBasis::Size(); const Matrix& Mnodes = Basis::Nodes(); for (Long t = 0; t < Nt; t++) { for (Long p = 0; p < Np; p++) { Real theta = t / (Real)Nt; Real phi = p / (Real)Np; Long i = (Long)(theta * Mt); Long j = (Long)(phi * Mp); Real x = theta * Mt - i; Real y = phi * Mp - j; Long elem_idx = i * Mp + j; Vector Interp0(ORDER); Vector Interp1(ORDER); { // Set Interp0, Interp1 auto node = [&Mnodes] (Long i) { return Mnodes[0][i]; }; for (Long i = 0; i < ORDER; i++) { Real wt_x = 1, wt_y = 1; for (Long j = 0; j < ORDER; j++) { if (j != i) { wt_x *= (x - node(j)) / (node(i) - node(j)); wt_y *= (y - node(j)) / (node(i) - node(j)); } Interp0[i] = wt_x; Interp1[i] = wt_y; } } } for (Long ii = 0; ii < ORDER; ii++) { for (Long jj = 0; jj < ORDER; jj++) { Long node_idx = jj * ORDER + ii; for (Long k = 0; k < dof; k++) { Xf[(k*Nt+t)*Np+p] += X[elem_idx*dof+k][node_idx] * Interp0[ii] * Interp1[jj]; } } } } } return Xf; }; auto grid2cheb = [] (const Vector& Xf, Long Nt, Long Np, Long Mt, Long Mp) { Long dof = Xf.Dim() / (Nt*Np); SCTL_ASSERT(Xf.Dim() == dof*Nt*Np); Vector X(Mt*Mp*dof); constexpr Integer INTERP_ORDER = 12; for (Long tt = 0; tt < Mt; tt++) { for (Long pp = 0; pp < Mp; pp++) { for (Long t = 0; t < ORDER; t++) { for (Long p = 0; p < ORDER; p++) { Matrix Mnodes = Basis::Nodes(); Real theta = (tt + Mnodes[0][t]) / Mt; Real phi = (pp + Mnodes[0][p]) / Mp; Long i = (Long)(theta * Nt); Long j = (Long)(phi * Np); Real x = theta * Nt - i; Real y = phi * Np - j; Vector Interp0(INTERP_ORDER); Vector Interp1(INTERP_ORDER); { // Set Interp0, Interp1 auto node = [] (Long i) { return (Real)i - (INTERP_ORDER-1)/2; }; for (Long i = 0; i < INTERP_ORDER; i++) { Real wt_x = 1, wt_y = 1; for (Long j = 0; j < INTERP_ORDER; j++) { if (j != i) { wt_x *= (x - node(j)) / (node(i) - node(j)); wt_y *= (y - node(j)) / (node(i) - node(j)); } Interp0[i] = wt_x; Interp1[i] = wt_y; } } } for (Long k = 0; k < dof; k++) { Real X0 = 0; for (Long ii = 0; ii < INTERP_ORDER; ii++) { for (Long jj = 0; jj < INTERP_ORDER; jj++) { Long idx_i = (i + ii-(INTERP_ORDER-1)/2 + Nt) % Nt; Long idx_j = (j + jj-(INTERP_ORDER-1)/2 + Np) % Np; X0 += Interp0[ii] * Interp1[jj] * Xf[(k*Nt+idx_i)*Np+idx_j]; } } Long elem_idx = tt * Mp + pp; Long node_idx = p * ORDER + t; X[elem_idx*dof+k][node_idx] = X0; } } } } } return X; }; Long Nelem = S.NElem(); if (Jp.Dim() != Nelem * COORD_DIM) Jp.ReInit(Nelem * COORD_DIM); if (Jt.Dim() != Nelem * COORD_DIM) Jt.ReInit(Nelem * COORD_DIM); for (Long k = 0; k < S.Nsurf(); k++) { Long Nt = S.NTor(k)*ORDER, Np = S.NPol(k)*ORDER; const auto& X_ = S.GetElemList().ElemVector(); Vector X(S.NTor(k)*S.NPol(k)*COORD_DIM, (Iterator)X_.begin()+S.ElemDsp(k)*COORD_DIM, false); biest::Surface SS(Nt, Np); biest::SurfaceOp surf_op(comm, Nt, Np); SS.Coord() = cheb2grid(X, S.NTor(k), S.NPol(k), Nt, Np); Vector dX, d2X; surf_op.Grad2D(dX, SS.Coord()); surf_op.Grad2D(d2X, dX); sctl::Vector Jt_(COORD_DIM * Nt * Np); sctl::Vector Jp_(COORD_DIM * Nt * Np); { // Set Jt_, Jp_ Vector DivV, InvLapDivV, GradInvLapDivV; for (sctl::Long i = 0; i < Nt*Np; i++) { // Set V for (sctl::Long k =0; k < COORD_DIM; k++) { Jt_[k * Nt*Np + i] = dX[(k*2+0) * Nt*Np + i]; Jp_[k * Nt*Np + i] = dX[(k*2+1) * Nt*Np + i]; } } surf_op.SurfDiv(DivV, dX, Jt_); surf_op.GradInvSurfLap(GradInvLapDivV, dX, d2X, DivV, gmres_tol * max_norm(Jt_) / max_norm(DivV), max_iter, 1.5); Jt_ = Jt_ - GradInvLapDivV; surf_op.SurfDiv(DivV, dX, Jp_); surf_op.GradInvSurfLap(GradInvLapDivV, dX, d2X, DivV, gmres_tol * max_norm(Jp_) / max_norm(DivV), max_iter, 1.5); Jp_ = Jp_ - GradInvLapDivV; } Vector Jt__(S.NTor(k)*S.NPol(k)*COORD_DIM, (Iterator)Jt.begin()+S.ElemDsp(k)*COORD_DIM, false); Vector Jp__(S.NTor(k)*S.NPol(k)*COORD_DIM, (Iterator)Jp.begin()+S.ElemDsp(k)*COORD_DIM, false); Jt__ = grid2cheb(Jt_, Nt, Np, S.NTor(k), S.NPol(k)); Jp__ = grid2cheb(Jp_, Nt, Np, S.NTor(k), S.NPol(k)); } } static void compute_norm_area_elem(const Stellarator& S, Vector& normal, Vector& area_elem){ // Set normal, area_elem const Vector& X = S.GetElemList().ElemVector(); const Long Nelem = X.Dim() / COORD_DIM; const Long Nnodes = ElemBasis::Size(); Vector dX; ElemBasis::Grad(dX, X); area_elem.ReInit(Nelem); normal.ReInit(Nelem * COORD_DIM); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { Tensor x, n; Tensor dx; x(0) = X[i*COORD_DIM+0][j]; x(1) = X[i*COORD_DIM+1][j]; x(2) = X[i*COORD_DIM+2][j]; dx(0,0) = dX[i*COORD_DIM*2+0][j]; dx(0,1) = dX[i*COORD_DIM*2+1][j]; dx(1,0) = dX[i*COORD_DIM*2+2][j]; dx(1,1) = dX[i*COORD_DIM*2+3][j]; dx(2,0) = dX[i*COORD_DIM*2+4][j]; dx(2,1) = dX[i*COORD_DIM*2+5][j]; n(0) = dx(1,0) * dx(2,1) - dx(2,0) * dx(1,1); n(1) = dx(2,0) * dx(0,1) - dx(0,0) * dx(2,1); n(2) = dx(0,0) * dx(1,1) - dx(1,0) * dx(0,1); Real area_elem_ = sqrt(n(0)*n(0) + n(1)*n(1) + n(2)*n(2)); Real ooae = 1 / area_elem_; n(0) *= ooae; n(1) *= ooae; n(2) *= ooae; normal[i*COORD_DIM+0][j] = n(0); normal[i*COORD_DIM+1][j] = n(1); normal[i*COORD_DIM+2][j] = n(2); area_elem[i][j] = area_elem_; } } } static Vector compute_B(const Stellarator& S, const Vector& sigma, Real alpha, Real beta) { const Long Nelem = S.NElem(); Vector B(S.NElem() * COORD_DIM); if (sigma.Dim()) { const Long Nnodes = ElemBasis::Size(); Vector normal, area_elem; compute_norm_area_elem(S, normal, area_elem); if (S.Nsurf() == 2) { Long Nelem0 = S.NTor(0)*S.NPol(0); for (Long i = 0; i < Nelem0*COORD_DIM; i++) { for (Long j = 0; j < Nnodes; j++) { normal[i][j] *= -1.0; } } } EvalQuadrature(B, S.quadrature_FxdU, S, sigma, S.Laplace_FxdU); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { for (Long k = 0; k < COORD_DIM; k++) { B[i*COORD_DIM+k][j] -= 0.5*sigma[i][j]*normal[i*COORD_DIM+k][j]; } } } } else { B = 0; } if (S.Nsurf() >= 1) B += S.Bt0*alpha; if (S.Nsurf() >= 2) B += S.Bp0*beta; return B; } static Vector compute_dB(const Stellarator& S, const Vector& sigma, Real alpha, Real beta) { const Long Nelem = S.NElem(); Vector dB(S.NElem() * COORD_DIM * COORD_DIM); if (sigma.Dim()) { EvalQuadrature(dB, S.quadrature_Fxd2U, S, sigma, S.Laplace_Fxd2U); } else { dB = 0; } if (S.Nsurf() >= 1) dB += S.dBt0*alpha; if (S.Nsurf() >= 2) dB += S.dBp0*beta; return dB; } static void compute_flux(Real& flux_tor, Real& flux_pol, const Stellarator& S, const Vector& B, const Vector& normal) { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); SCTL_ASSERT(B.Dim() == Nelem*COORD_DIM); SCTL_ASSERT(normal.Dim() == Nelem*COORD_DIM); Vector J(Nelem * COORD_DIM); for (Long i = 0; i < Nelem; i++) { // Set J for (Long j = 0; j < Nnodes; j++) { Tensor b, n; b(0) = B[i*COORD_DIM+0][j]; b(1) = B[i*COORD_DIM+1][j]; b(2) = B[i*COORD_DIM+2][j]; n(0) = normal[i*COORD_DIM+0][j]; n(1) = normal[i*COORD_DIM+1][j]; n(2) = normal[i*COORD_DIM+2][j]; J[i*COORD_DIM+0][j] = n(1) * b(2) - n(2) * b(1); J[i*COORD_DIM+1][j] = n(2) * b(0) - n(0) * b(2); J[i*COORD_DIM+2][j] = n(0) * b(1) - n(1) * b(0); } } Vector A; EvalQuadrature(A, S.quadrature_FxU, S, J, S.Laplace_FxU); Vector circ_pol(S.Nsurf()), circ_tor(S.Nsurf()); { // compute circ Vector dX; ElemBasis::Grad(dX, S.GetElemList().ElemVector()); const auto& quad_wts = ElemBasis::QuadWts(); for (Long k = 0; k < S.Nsurf(); k++) { circ_pol[k] = 0; circ_tor[k] = 0; Long Ndsp = S.ElemDsp(k); for (Long i = 0; i < S.NTor(k)*S.NPol(k); i++) { for (Long j = 0; j < Nnodes; j++) { circ_pol[k] += A[(Ndsp+i)*COORD_DIM+0][j] * dX[(Ndsp+i)*COORD_DIM*2+1][j] * quad_wts[j] / S.NTor(k); circ_pol[k] += A[(Ndsp+i)*COORD_DIM+1][j] * dX[(Ndsp+i)*COORD_DIM*2+3][j] * quad_wts[j] / S.NTor(k); circ_pol[k] += A[(Ndsp+i)*COORD_DIM+2][j] * dX[(Ndsp+i)*COORD_DIM*2+5][j] * quad_wts[j] / S.NTor(k); circ_tor[k] += A[(Ndsp+i)*COORD_DIM+0][j] * dX[(Ndsp+i)*COORD_DIM*2+0][j] * quad_wts[j] / S.NPol(k); circ_tor[k] += A[(Ndsp+i)*COORD_DIM+1][j] * dX[(Ndsp+i)*COORD_DIM*2+2][j] * quad_wts[j] / S.NPol(k); circ_tor[k] += A[(Ndsp+i)*COORD_DIM+2][j] * dX[(Ndsp+i)*COORD_DIM*2+4][j] * quad_wts[j] / S.NPol(k); } } } } if (S.Nsurf() == 1) { flux_tor = circ_pol[0]; flux_pol = 0; } else if (S.Nsurf() == 2) { flux_tor = circ_pol[1] - circ_pol[0]; flux_pol = circ_tor[0] - circ_tor[1]; } else { SCTL_ASSERT(false); } } static Vector compute_A(const Stellarator& S, const Vector& x) { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); SCTL_ASSERT(x.Dim() == Nelem*Nnodes+S.Nsurf()); Vector normal, area_elem; compute_norm_area_elem(S, normal, area_elem); if (S.Nsurf() == 2) { Long Nelem0 = S.NTor(0)*S.NPol(0); for (Long i = 0; i < Nelem0*COORD_DIM; i++) { for (Long j = 0; j < Nnodes; j++) { normal[i][j] *= -1.0; } } } Vector sigma(Nelem); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { sigma[i][j] = x[i*Nnodes+j]; } } Real alpha = (S.Nsurf() >= 1 ? x[Nelem*Nnodes + 0] : 0); Real beta = (S.Nsurf() >= 2 ? x[Nelem*Nnodes + 1] : 0); Vector B = compute_B(S, sigma, alpha, beta); Vector BdotN = compute_dot_prod(B, normal); Real flux_tor, flux_pol; compute_flux(flux_tor, flux_pol, S, B, normal); Vector Ax(Nelem*Nnodes+S.Nsurf()); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { Ax[i*Nnodes+j] = BdotN[i][j]; } } if (S.Nsurf() >= 1) Ax[Nelem*Nnodes + 0] = flux_tor; if (S.Nsurf() >= 2) Ax[Nelem*Nnodes + 1] = flux_pol; return Ax; } static void compute_invA(Vector& sigma, Real& alpha, Real& beta, const Stellarator& S, Real flux_tor, Real flux_pol, const Comm& comm) { typename sctl::ParallelSolver::ParallelOp BIOp = [&S](sctl::Vector* Ax, const sctl::Vector& x) { (*Ax) = compute_A(S, x); }; const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector rhs_(Nelem * Nnodes + S.Nsurf()); rhs_ = 0; if (S.Nsurf() >= 1) rhs_[Nelem * Nnodes + 0] = flux_tor; if (S.Nsurf() >= 2) rhs_[Nelem * Nnodes + 1] = flux_pol; Vector x_(Nelem * Nnodes + S.Nsurf()); x_ = 0; ParallelSolver linear_solver(comm, true); linear_solver(&x_, BIOp, rhs_, 1e-6, 100); sigma.ReInit(Nelem); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { sigma[i][j] = x_[i*Nnodes+j]; } } alpha = (S.Nsurf() >= 1 ? x_[Nelem * Nnodes + 0] : 0); beta = (S.Nsurf() >= 2 ? x_[Nelem * Nnodes + 1] : 0); } static Vector compute_Aadj(const Stellarator& S, const Vector& x) { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); SCTL_ASSERT(x.Dim() == Nelem*Nnodes+S.Nsurf()); Vector normal, area_elem; compute_norm_area_elem(S, normal, area_elem); if (S.Nsurf() == 2) { Long Nelem0 = S.NTor(0)*S.NPol(0); for (Long i = 0; i < Nelem0*COORD_DIM; i++) { for (Long j = 0; j < Nnodes; j++) { normal[i][j] *= -1.0; } } } Vector x0(Nelem); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { x0[i][j] = x[i*Nnodes+j]; } } Real x1 = (S.Nsurf() >= 1 ? x[Nelem*Nnodes + 0] : 0); Real x2 = (S.Nsurf() >= 2 ? x[Nelem*Nnodes + 1] : 0); Vector Ax0; Real Ax1, Ax2; { // Set Ax0, Ax1, Ax2 Vector x0_n(Nelem*COORD_DIM); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { x0_n[i*COORD_DIM+0][j] = x0[i][j] * normal[i*COORD_DIM+0][j]; x0_n[i*COORD_DIM+1][j] = x0[i][j] * normal[i*COORD_DIM+1][j]; x0_n[i*COORD_DIM+2][j] = x0[i][j] * normal[i*COORD_DIM+2][j]; } } EvalQuadrature(Ax0, S.quadrature_dUxF, S, x0_n, S.Laplace_dUxF); Ax0 = x0*(-0.5) - Ax0; Ax1 = (S.Nsurf() >= 1 ? compute_inner_prod(area_elem, compute_dot_prod(S.Bt0, normal), x0) : 0); Ax2 = (S.Nsurf() >= 2 ? compute_inner_prod(area_elem, compute_dot_prod(S.Bp0, normal), x0) : 0); } // TODO: precompute A21adj, A22adj auto compute_A21adj = [&S,&normal,&area_elem] (bool toroidal_flux) { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector density(Nelem * COORD_DIM); { // Set density Real scal[2]; if (S.Nsurf() == 1) { SCTL_ASSERT(toroidal_flux == true); scal[0] = 1.0 / S.NTor(0); scal[1] = 0; } else if (S.Nsurf() == 2) { if (toroidal_flux == true) { scal[0] = -1.0 / S.NTor(0); scal[1] = 1.0 / S.NTor(1); } else { scal[0] = 1.0 / S.NPol(0); scal[1] = -1.0 / S.NPol(1); } } else { SCTL_ASSERT(false); } Vector dX; ElemBasis::Grad(dX, S.GetElemList().ElemVector()); for (Long k = 0; k < S.Nsurf(); k++) { for (Long i_ = 0; i_ < S.NTor(k)*S.NPol(k); i_++) { Long i = S.ElemDsp(k) + i_; for (Long j = 0; j < Nnodes; j++) { Real s = scal[k] / area_elem[i][j]; density[i*COORD_DIM+0][j] = dX[i*COORD_DIM*2+0+(toroidal_flux?1:0)][j] * s; density[i*COORD_DIM+1][j] = dX[i*COORD_DIM*2+2+(toroidal_flux?1:0)][j] * s; density[i*COORD_DIM+2][j] = dX[i*COORD_DIM*2+4+(toroidal_flux?1:0)][j] * s; } } } } Vector Gdensity, nxGdensity(Nelem * COORD_DIM), A21adj; EvalQuadrature(Gdensity, S.quadrature_FxU, S, density, S.Laplace_FxU); for (Long i = 0; i < Nelem; i++) { // Set nxGdensity for (Long j = 0; j < Nnodes; j++) { Tensor Gdensity_, n; Gdensity_(0) = Gdensity[i*COORD_DIM+0][j]; Gdensity_(1) = Gdensity[i*COORD_DIM+1][j]; Gdensity_(2) = Gdensity[i*COORD_DIM+2][j]; n(0) = normal[i*COORD_DIM+0][j]; n(1) = normal[i*COORD_DIM+1][j]; n(2) = normal[i*COORD_DIM+2][j]; nxGdensity[i*COORD_DIM+0][j] = n(1) * Gdensity_(2) - n(2) * Gdensity_(1); nxGdensity[i*COORD_DIM+1][j] = n(2) * Gdensity_(0) - n(0) * Gdensity_(2); nxGdensity[i*COORD_DIM+2][j] = n(0) * Gdensity_(1) - n(1) * Gdensity_(0); } } EvalQuadrature(A21adj, S.quadrature_dUxF, S, nxGdensity, S.Laplace_dUxF); return A21adj; }; if (S.Nsurf() >= 1) Ax0 += compute_A21adj( true) * x1; if (S.Nsurf() >= 2) Ax0 += compute_A21adj(false) * x2; if (S.Nsurf() == 1) { // Add flux part of Ax1, Ax2 Real flux_tor, flux_pol; compute_flux(flux_tor, flux_pol, S, S.Bt0, normal); Ax1 += flux_tor * x1; Ax2 += 0; } else if (S.Nsurf() == 2) { Real flux_tor, flux_pol; compute_flux(flux_tor, flux_pol, S, S.Bt0, normal); Ax1 += flux_tor * x1 + flux_pol * x2; compute_flux(flux_tor, flux_pol, S, S.Bp0, normal); Ax2 += flux_tor * x1 + flux_pol * x2; } else { SCTL_ASSERT(false); } Vector Ax(Nelem*Nnodes+S.Nsurf()); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { Ax[i*Nnodes+j] = Ax0[i][j]; } } if (S.Nsurf() >= 1) Ax[Nelem*Nnodes + 0] = Ax1; if (S.Nsurf() >= 2) Ax[Nelem*Nnodes + 1] = Ax2; return Ax; } static Vector compute_invAadj(const Stellarator& S, const Vector& b, const Comm& comm) { typename sctl::ParallelSolver::ParallelOp BIOp = [&S](sctl::Vector* Ax, const sctl::Vector& x) { (*Ax) = compute_Aadj(S,x); }; const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector x(b.Dim()); x = 0; ParallelSolver linear_solver(comm, true); linear_solver(&x, BIOp, b, 1e-6, 100); return x; } static Vector compute_pressure_jump(const Vector>& Svec, const Vector& pressure) { const Long Nnodes = ElemBasis::Size(); const Long Nsurf = Svec.Dim(); SCTL_ASSERT(pressure.Dim() == Nsurf); Vector> total_pressure(Nsurf); for (Long i = 0; i < Nsurf; i++) { // Set total_pressure const Long Nelem = Svec[i].NElem(); const auto& B = Svec[i].B; total_pressure[i].ReInit(Nelem); for (Long j = 0; j < Nelem; j++) { for (Long k = 0; k < Nnodes; k++) { Real B2 = 0; B2 += B[j*COORD_DIM+0][k] * B[j*COORD_DIM+0][k]; B2 += B[j*COORD_DIM+1][k] * B[j*COORD_DIM+1][k]; B2 += B[j*COORD_DIM+2][k] * B[j*COORD_DIM+2][k]; total_pressure[i][j][k] = 0.5 * B2 + pressure[i]; } } } Vector elem_cnt, elem_dsp; for (Long i = 0; i < Nsurf; i++) { if (i == 0) { elem_cnt.PushBack(Svec[i].NTor(0) * Svec[i].NPol(0)); elem_dsp.PushBack(0); } else { elem_cnt.PushBack(Svec[i].NTor(1) * Svec[i].NPol(1)); elem_dsp.PushBack(elem_dsp[i-1] + elem_cnt[i-1]); } } Vector pressure_jump(elem_dsp[Nsurf-1] + elem_cnt[Nsurf-1]); pressure_jump = 0; for (Long i = 0; i < Nsurf-1; i++) { // Set pressure_jump Long Nelem, offset; if (i == 0) { offset = 0; Nelem = Svec[i].NTor(0) * Svec[i].NPol(0); } else { offset = Svec[i].NTor(0) * Svec[i].NPol(0); Nelem = Svec[i].NTor(1) * Svec[i].NPol(1); } for (Long j = 0; j < Nelem; j++) { for (Long k = 0; k < Nnodes; k++) { Real T0 = total_pressure[i][offset+j][k]; Real T1 = (i+1>& Svec, const Vector& pressure) { Vector pressure_jump = compute_pressure_jump(Svec, pressure); const Long Nnodes = ElemBasis::Size(); const Long Nsurf = Svec.Dim(); Long elem_offset = 0; for (Long i = 0; i < Nsurf; i++) { // Allocate Svec[i].gvec.ReInit(Svec[i].NElem()); Svec[i].gvec = 0; } for (Long i = 0; i < Nsurf-1; i++) { // Set gvec Long Nelem, offset; if (i == 0) { offset = 0; Nelem = Svec[i].NTor(0) * Svec[i].NPol(0); } else { offset = Svec[i].NTor(0) * Svec[i].NPol(0); Nelem = Svec[i].NTor(1) * Svec[i].NPol(1); } for (Long j = 0; j < Nelem; j++) { for (Long k = 0; k < Nnodes; k++) { Real jump = pressure_jump[elem_offset+j][k]; Svec[i].gvec[offset+j][k] = 0.5 * jump * jump; if (i+1>& Svec, const Vector& pressure) { Vector pressure_jump = compute_pressure_jump(Svec, pressure); const Long Nnodes = ElemBasis::Size(); const Long Nsurf = Svec.Dim(); Long elem_offset = 0; for (Long i = 0; i < Nsurf; i++) { // Allocate Svec[i].dgdB.ReInit(Svec[i].NElem() * COORD_DIM); Svec[i].dgdB = 0; } for (Long i = 0; i < Nsurf-1; i++) { // Set dgdB Long Nelem, offset; if (i == 0) { offset = 0; Nelem = Svec[i].NTor(0) * Svec[i].NPol(0); } else { offset = Svec[i].NTor(0) * Svec[i].NPol(0); Nelem = Svec[i].NTor(1) * Svec[i].NPol(1); } for (Long j = 0; j < Nelem; j++) { for (Long k = 0; k < Nnodes; k++) { Real jump = pressure_jump[elem_offset+j][k]; Svec[i].dgdB[(offset+j)*COORD_DIM+0][k] = -jump * Svec[i].B[(offset+j)*COORD_DIM+0][k]; Svec[i].dgdB[(offset+j)*COORD_DIM+1][k] = -jump * Svec[i].B[(offset+j)*COORD_DIM+1][k]; Svec[i].dgdB[(offset+j)*COORD_DIM+2][k] = -jump * Svec[i].B[(offset+j)*COORD_DIM+2][k]; if (i+1>& Svec, const Vector& pressure) { Real g = 0; compute_gvec(Svec, pressure); for (Long i = 0; i < Svec.Dim(); i++) { // Set gvec Vector normal, area_elem, wt(Svec[i].NElem()); compute_norm_area_elem(Svec[i], normal, area_elem); wt = 0.5; if (i == Svec.Dim()-1) { Long Nsurf = Svec[i].Nsurf(); Long Nelem = Svec[i].NTor(Nsurf-1) * Svec[i].NPol(Nsurf-1); Long offset = Svec[i].ElemDsp(Nsurf-1); for (Long j = 0; j < Nelem; j++) { wt[offset + j] = 1.0; } } g += compute_inner_prod(area_elem, Svec[i].gvec, wt); } return g; } Stellarator(const Vector& NtNp = Vector()) { NtNp_ = NtNp; Long Nsurf = NtNp_.Dim() / 2; SCTL_ASSERT(Nsurf*2 == NtNp_.Dim()); Long Nelem = 0; elem_dsp.ReInit(Nsurf+1); elem_dsp[0] = 0; for (Long i = 0; i < Nsurf; i++) { Nelem += NtNp_[i*2+0]*NtNp_[i*2+1]; elem_dsp[i+1] = Nelem; } elements.ReInit(Nelem); for (Long i = 0; i < Nsurf; i++) { InitSurf(i, this->Nsurf()); } } Long ElemIdx(Long s, Long t, Long p) { SCTL_ASSERT(0 <= s && s < Nsurf()); SCTL_ASSERT(0 <= t && t < NtNp_[s*2+0]); SCTL_ASSERT(0 <= p && p < NtNp_[s*2+1]); return elem_dsp[s] + t*NtNp_[s*2+1] + p; } ElemBasis& Elem(Long elem, Integer dim) { return elements(elem,dim); } const ElemBasis& Elem(Long elem, Integer dim) const { return elements(elem,dim); } const ElemLst& GetElemList() const { return elements; } Long Nsurf() const { return elem_dsp.Dim()-1; } Long ElemDsp(Long s) const { return elem_dsp[s]; } Long NTor(Long s) const { return NtNp_[s*2+0]; } Long NPol(Long s) const { return NtNp_[s*2+1]; } Long NElem() const { return elements.NElem(); } static Vector compute_gradient(const Stellarator& S_, const Vector& pressure, const Vector& flux_tor_, const Vector& flux_pol_, Real* g_ptr = nullptr) { Comm comm = Comm::World(); Vector> Svec(S_.Nsurf()); for (Long i = 0; i < S_.Nsurf(); i++) { // Set Svec[i] (quadratures, B) const Long elem_dsp = (i==0 ? 0 : S_.ElemDsp(i-1)); const Long Nnodes = ElemBasis::Size(); Stellarator& S = Svec[i]; if (i == 0) { // Init S Vector NtNp; NtNp.PushBack(S_.NTor(i)); NtNp.PushBack(S_.NPol(i)); S = Stellarator(NtNp); } else { Vector NtNp; NtNp.PushBack(S_.NTor(i-1)); NtNp.PushBack(S_.NPol(i-1)); NtNp.PushBack(S_.NTor(i)); NtNp.PushBack(S_.NPol(i)); S = Stellarator(NtNp); } for (Long j = 0; j < S.NElem(); j++) { // Set S coordinates for (Long k = 0; k < Nnodes; k++) { S.Elem(j,0)[k] = S_.Elem(elem_dsp+j,0)[k]; S.Elem(j,1)[k] = S_.Elem(elem_dsp+j,1)[k]; S.Elem(j,2)[k] = S_.Elem(elem_dsp+j,2)[k]; } } SetupQuadrature(S.quadrature_dBS , S, S.BiotSavartGrad, order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER)); SetupQuadrature(S.quadrature_BS , S, S.BiotSavart , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER)); SetupQuadrature(S.quadrature_FxU , S, S.Laplace_FxU , order_singular, order_direct, -1.0, comm); SetupQuadrature(S.quadrature_FxdU , S, S.Laplace_FxdU , order_singular, order_direct, -1.0, comm); SetupQuadrature(S.quadrature_dUxF , S, S.Laplace_dUxF , order_singular, order_direct, -1.0, comm); SetupQuadrature(S.quadrature_dUxD , S, S.Laplace_dUxD , order_singular, order_direct, -1.0, comm, 0.01 * pow<-2,Real>(ORDER)); SetupQuadrature(S.quadrature_Fxd2U, S, S.Laplace_Fxd2U , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER)); { // Set Bt0, Bp0, dBt0, dBp0 Vector Jt, Jp; compute_harmonic_vector_potentials(Jt, Jp, S); EvalQuadrature(S.Bt0 , S.quadrature_BS , S, Jp, S.BiotSavart); EvalQuadrature(S.Bp0 , S.quadrature_BS , S, Jt, S.BiotSavart); EvalQuadrature(S.dBt0, S.quadrature_dBS, S, Jp, S.BiotSavartGrad); EvalQuadrature(S.dBp0, S.quadrature_dBS, S, Jt, S.BiotSavartGrad); } compute_invA(S.sigma, S.alpha, S.beta, S, flux_tor_[i], flux_pol_[i], comm); S.B = compute_B(S, S.sigma, S.alpha, S.beta); if (0) { // Write VTU VTUData vtu; vtu.AddElems(S.GetElemList(), S.sigma, ORDER); vtu.WriteVTK("sigma"+std::to_string(i), comm); } if (0) { // Write VTU VTUData vtu; vtu.AddElems(S.GetElemList(), S.B, ORDER); vtu.WriteVTK("B"+std::to_string(i), comm); } } compute_gvec(Svec, pressure); compute_dgdB(Svec, pressure); if (g_ptr != nullptr) g_ptr[0] = compute_g(Svec, pressure); auto compute_gradient = [&comm] (const Stellarator& S) { const Long Nnodes = ElemBasis::Size(); const Long Nelem = S.NElem(); const auto& sigma = S.sigma; const auto& alpha = S.alpha; const auto& beta = S.beta; const auto& B = S.B; Vector normal, area_elem; compute_norm_area_elem(S, normal, area_elem); if (S.Nsurf() == 2) { Long Nelem0 = S.NTor(0)*S.NPol(0); for (Long i = 0; i < Nelem0*COORD_DIM; i++) { for (Long j = 0; j < Nnodes; j++) { normal[i][j] *= -1.0; } } } auto compute_H = [] (const ElemList& elem_lst, const Vector& normal) { const Long Nnodes = ElemBasis::Size(); const Long Nelem = elem_lst.NElem(); const Vector X = elem_lst.ElemVector(); Vector dX, d2X, H(Nelem); ElemBasis::Grad(dX, X); ElemBasis::Grad(d2X, dX); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { Tensor I, invI, II; for (Long k0 = 0; k0 < 2; k0++) { for (Long k1 = 0; k1 < 2; k1++) { I(k0,k1) = 0; I(k0,k1) += dX[(i*COORD_DIM+0)*2+k0][j] * dX[(i*COORD_DIM+0)*2+k1][j]; I(k0,k1) += dX[(i*COORD_DIM+1)*2+k0][j] * dX[(i*COORD_DIM+1)*2+k1][j]; I(k0,k1) += dX[(i*COORD_DIM+2)*2+k0][j] * dX[(i*COORD_DIM+2)*2+k1][j]; II(k0,k1) = 0; II(k0,k1) += d2X[(i*COORD_DIM+0)*4+k0*2+k1][j] * normal[i*COORD_DIM+0][j]; II(k0,k1) += d2X[(i*COORD_DIM+1)*4+k0*2+k1][j] * normal[i*COORD_DIM+1][j]; II(k0,k1) += d2X[(i*COORD_DIM+2)*4+k0*2+k1][j] * normal[i*COORD_DIM+2][j]; } } { // Set invI Real detI = I(0,0)*I(1,1)-I(0,1)*I(1,0); invI(0,0) = I(1,1) / detI; invI(0,1) = -I(0,1) / detI; invI(1,0) = -I(1,0) / detI; invI(1,1) = I(0,0) / detI; } { // Set H H[i][j] = 0; H[i][j] += -0.5 * II(0,0)*invI(0,0); H[i][j] += -0.5 * II(0,1)*invI(0,1); H[i][j] += -0.5 * II(1,0)*invI(1,0); H[i][j] += -0.5 * II(1,1)*invI(1,1); } } } return H; }; Vector H = compute_H(S.GetElemList(), normal); auto compute_dg_dnu = [&S,&normal,&area_elem,&H]() { // dg_dnu = (B*B) 2H - (2 B) \cdot (n \cdnot nabla) \nabla G[sigma] + (2 B) \alpha dB0_dnu \hat{\theta} + sigma (\nabla D)^T [2 B] + (2H) sigma (\nabla G)^T [2 B] const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); const Vector& gvec = S.gvec; const Vector& v = S.dgdB; const auto& sigma = S.sigma; const auto& alpha = S.alpha; const auto& beta = S.beta; const auto& B = S.B; Vector dg_dnu0(Nelem), dg_dnu1(Nelem), dg_dnu2(Nelem), dg_dnu3(Nelem), dg_dnu4(Nelem); dg_dnu0 = 0; dg_dnu1 = 0; dg_dnu2 = 0; dg_dnu3 = 0; dg_dnu4 = 0; // dg_dnu0 = (B*B) 2H for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dg_dnu0[i][j] = gvec[i][j] * (2.0*H[i][j]) * 0.5; // multiplicative factor 0.5 is there so that this term is not // counted twice from shape derivative of regions on either side // of the domain. } } // dg_dnu1 = (2 B) \cdot (n \cdnot \nabla) B Vector dB = compute_dB(S, sigma, alpha, beta); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dg_dnu1[i][j] = 0; dg_dnu1[i][j] -= dB[i*9+0][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+0][j]; dg_dnu1[i][j] -= dB[i*9+1][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+0][j]; dg_dnu1[i][j] -= dB[i*9+2][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+0][j]; dg_dnu1[i][j] -= dB[i*9+3][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+1][j]; dg_dnu1[i][j] -= dB[i*9+4][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+1][j]; dg_dnu1[i][j] -= dB[i*9+5][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+1][j]; dg_dnu1[i][j] -= dB[i*9+6][j] * normal[i*COORD_DIM+0][j] * v[i*COORD_DIM+2][j]; dg_dnu1[i][j] -= dB[i*9+7][j] * normal[i*COORD_DIM+1][j] * v[i*COORD_DIM+2][j]; dg_dnu1[i][j] -= dB[i*9+8][j] * normal[i*COORD_DIM+2][j] * v[i*COORD_DIM+2][j]; } } // dg_dnu3 = (sigma (\nabla D)^T [2 B] Vector nablaDtv; EvalQuadrature(nablaDtv, S.quadrature_dUxD, S, v, S.Laplace_dUxD); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dg_dnu3[i][j] = 0; dg_dnu3[i][j] += sigma[i][j] * nablaDtv[i*COORD_DIM+0][j]*normal[i*COORD_DIM+0][j]; dg_dnu3[i][j] += sigma[i][j] * nablaDtv[i*COORD_DIM+1][j]*normal[i*COORD_DIM+1][j]; dg_dnu3[i][j] += sigma[i][j] * nablaDtv[i*COORD_DIM+2][j]*normal[i*COORD_DIM+2][j]; } } // dg_dnu4 = (2H) sigma (\nabla G)^T [2 B] EvalQuadrature(dg_dnu4, S.quadrature_dUxF, S, v, S.Laplace_dUxF); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dg_dnu4[i][j] += 0.5 * v[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j]; dg_dnu4[i][j] += 0.5 * v[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j]; dg_dnu4[i][j] += 0.5 * v[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j]; dg_dnu4[i][j] *= 2*H[i][j] * sigma[i][j]; } } return dg_dnu0 + dg_dnu1 + dg_dnu3 - dg_dnu4; }; Vector dg_dnu = compute_dg_dnu(); auto compute_dg_dsigma = [&S,&normal,&area_elem] () { const Long Nnodes = ElemBasis::Size(); const Long Nelem = S.NElem(); const auto& B = S.B; const Vector& dgdB = S.dgdB; auto compute_dg_dsigma = [&S,&B,&dgdB,&normal]() { // dg_dsigma = \int 2 B \cdot (\nabla G + n/2) Vector B_dot_gradG; EvalQuadrature(B_dot_gradG, S.quadrature_dUxF, S, dgdB, S.Laplace_dUxF); return B_dot_gradG * (-1.0) + compute_dot_prod(dgdB,normal) * 0.5; }; auto compute_dg_dalpha = [&S,&B,&dgdB,&area_elem] () { auto dB_dalpha = compute_B(S, Vector(),1,0); return compute_inner_prod(area_elem, dgdB,dB_dalpha); }; auto compute_dg_dbeta = [&S,&B,&dgdB,&area_elem] () { auto dB_dalpha = compute_B(S, Vector(),0,1); return compute_inner_prod(area_elem, dgdB,dB_dalpha); }; Vector dg_dsigma(Nelem*Nnodes+S.Nsurf()); Vector dg_dsigma_ = compute_dg_dsigma(); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dg_dsigma[i*Nnodes+j] = dg_dsigma_[i][j]; } } if (S.Nsurf() >= 1) dg_dsigma[Nelem*Nnodes+0] = compute_dg_dalpha(); if (S.Nsurf() >= 2) dg_dsigma[Nelem*Nnodes+1] = compute_dg_dbeta (); return dg_dsigma; }; Vector dg_dsigma = compute_dg_dsigma(); Vector dg_dsigma_invA = compute_invAadj(S, dg_dsigma, comm); /////////////////////////////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////////////////////////////// auto compute_grad_adj = [&S,&area_elem] (const Vector& V) { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector du_dX(Nelem*COORD_DIM*2); { // Set du_dX Vector dX; ElemBasis::Grad(dX, S.GetElemList().ElemVector()); auto inv2x2 = [](Tensor M) { Tensor Mout; Real oodet = 1 / (M(0,0) * M(1,1) - M(0,1) * M(1,0)); Mout(0,0) = M(1,1) * oodet; Mout(0,1) = -M(0,1) * oodet; Mout(1,0) = -M(1,0) * oodet; Mout(1,1) = M(0,0) * oodet; return Mout; }; for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { Tensor dX_du; dX_du(0,0) = dX[(i*COORD_DIM+0)*2+0][j]; dX_du(1,0) = dX[(i*COORD_DIM+1)*2+0][j]; dX_du(2,0) = dX[(i*COORD_DIM+2)*2+0][j]; dX_du(0,1) = dX[(i*COORD_DIM+0)*2+1][j]; dX_du(1,1) = dX[(i*COORD_DIM+1)*2+1][j]; dX_du(2,1) = dX[(i*COORD_DIM+2)*2+1][j]; Tensor G; // = dX_du.Transpose() * dX_du; G(0,0) = dX_du(0,0) * dX_du(0,0) + dX_du(1,0) * dX_du(1,0) + dX_du(2,0) * dX_du(2,0); G(0,1) = dX_du(0,0) * dX_du(0,1) + dX_du(1,0) * dX_du(1,1) + dX_du(2,0) * dX_du(2,1); G(1,0) = dX_du(0,1) * dX_du(0,0) + dX_du(1,1) * dX_du(1,0) + dX_du(2,1) * dX_du(2,0); G(1,1) = dX_du(0,1) * dX_du(0,1) + dX_du(1,1) * dX_du(1,1) + dX_du(2,1) * dX_du(2,1); Tensor Ginv = inv2x2(G); du_dX[(i*COORD_DIM+0)*2+0][j] = Ginv(0,0) * dX_du(0,0) + Ginv(0,1) * dX_du(0,1); du_dX[(i*COORD_DIM+1)*2+0][j] = Ginv(0,0) * dX_du(1,0) + Ginv(0,1) * dX_du(1,1); du_dX[(i*COORD_DIM+2)*2+0][j] = Ginv(0,0) * dX_du(2,0) + Ginv(0,1) * dX_du(2,1); du_dX[(i*COORD_DIM+0)*2+1][j] = Ginv(1,0) * dX_du(0,0) + Ginv(1,1) * dX_du(0,1); du_dX[(i*COORD_DIM+1)*2+1][j] = Ginv(1,0) * dX_du(1,0) + Ginv(1,1) * dX_du(1,1); du_dX[(i*COORD_DIM+2)*2+1][j] = Ginv(1,0) * dX_du(2,0) + Ginv(1,1) * dX_du(2,1); } } } Vector dudX_V(Nelem*2); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dudX_V[i*2+0][j] = 0; dudX_V[i*2+1][j] = 0; dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+0)*2+0][j] * V[i*COORD_DIM+0][j] * area_elem[i][j]; dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+1)*2+0][j] * V[i*COORD_DIM+1][j] * area_elem[i][j]; dudX_V[i*2+0][j] += du_dX[(i*COORD_DIM+2)*2+0][j] * V[i*COORD_DIM+2][j] * area_elem[i][j]; dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+0)*2+1][j] * V[i*COORD_DIM+0][j] * area_elem[i][j]; dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+1)*2+1][j] * V[i*COORD_DIM+1][j] * area_elem[i][j]; dudX_V[i*2+1][j] += du_dX[(i*COORD_DIM+2)*2+1][j] * V[i*COORD_DIM+2][j] * area_elem[i][j]; } } Vector grad_dudX_V; ElemBasis::Grad(grad_dudX_V, dudX_V); Vector grad_adj_V(Nelem); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { grad_adj_V[i][j] = -(grad_dudX_V[(i*2+0)*2+0][j] + grad_dudX_V[(i*2+1)*2+1][j]) / area_elem[i][j]; } } return grad_adj_V; }; auto compute_u_dAdnu_v_0 = [&S,&normal,&H,&compute_grad_adj] (const Vector& u_, const Vector& v, Real alpha, Real beta) { const Long Nnodes = ElemBasis::Size(); const Long Nelem = S.NElem(); Vector dAdnu0(Nelem), dAdnu1(Nelem), dAdnu2(Nelem), dAdnu3(Nelem); Vector u(Nelem), u_n(Nelem*COORD_DIM); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { u[i][j] = u_[i*Nnodes+j]; u_n[i*COORD_DIM+0][j] = u[i][j] * normal[i*COORD_DIM+0][j]; u_n[i*COORD_DIM+1][j] = u[i][j] * normal[i*COORD_DIM+1][j]; u_n[i*COORD_DIM+2][j] = u[i][j] * normal[i*COORD_DIM+2][j]; } } // dAdnu0 = u B \cdot grad_nu Vector B = compute_B(S, v, alpha, beta); Vector u_B(Nelem*COORD_DIM); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { u_B[i*COORD_DIM+0][j] = u[i][j] * B[i*COORD_DIM+0][j]; u_B[i*COORD_DIM+1][j] = u[i][j] * B[i*COORD_DIM+1][j]; u_B[i*COORD_DIM+2][j] = u[i][j] * B[i*COORD_DIM+2][j]; } } dAdnu0 = compute_grad_adj(u_B)*(-1.0); // dAdnu1 = (u n) \cdot (n \cdnot \nabla) B Vector dB = compute_dB(S, v, alpha, beta); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dAdnu1[i][j] = 0; dAdnu1[i][j] -= dB[i*9+0][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+0][j]; dAdnu1[i][j] -= dB[i*9+1][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+1][j]; dAdnu1[i][j] -= dB[i*9+2][j] * normal[i*COORD_DIM+0][j] * u_n[i*COORD_DIM+2][j]; dAdnu1[i][j] -= dB[i*9+3][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+0][j]; dAdnu1[i][j] -= dB[i*9+4][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+1][j]; dAdnu1[i][j] -= dB[i*9+5][j] * normal[i*COORD_DIM+1][j] * u_n[i*COORD_DIM+2][j]; dAdnu1[i][j] -= dB[i*9+6][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+0][j]; dAdnu1[i][j] -= dB[i*9+7][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+1][j]; dAdnu1[i][j] -= dB[i*9+8][j] * normal[i*COORD_DIM+2][j] * u_n[i*COORD_DIM+2][j]; } } // dAdnu2 = (2H) v (I/2 + \nabla G)^T [u n] EvalQuadrature(dAdnu2, S.quadrature_dUxF, S, u_n, S.Laplace_dUxF); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dAdnu2[i][j] += 0.5 * u_n[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j]; dAdnu2[i][j] += 0.5 * u_n[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j]; dAdnu2[i][j] += 0.5 * u_n[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j]; dAdnu2[i][j] *= -2*H[i][j] * v[i][j]; } } // dAdnu3 = (v n \cdot \nabla D[u] Vector nablaDt_u_n; EvalQuadrature(nablaDt_u_n, S.quadrature_dUxD, S, u_n, S.Laplace_dUxD); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dAdnu3[i][j] = 0; dAdnu3[i][j] += v[i][j] * nablaDt_u_n[i*COORD_DIM+0][j]*normal[i*COORD_DIM+0][j]; dAdnu3[i][j] += v[i][j] * nablaDt_u_n[i*COORD_DIM+1][j]*normal[i*COORD_DIM+1][j]; dAdnu3[i][j] += v[i][j] * nablaDt_u_n[i*COORD_DIM+2][j]*normal[i*COORD_DIM+2][j]; } } return dAdnu0 + dAdnu1 + dAdnu2 + dAdnu3; }; auto compute_u_dAdnu_v_1 = [&S,&area_elem,&normal,&H,&compute_grad_adj] (const Vector& sigma, Real alpha, Real beta, bool toroidal_flux) { const Long Nnodes = ElemBasis::Size(); const Long Nelem = S.NElem(); Vector B = compute_B(S, sigma, alpha, beta); Vector gradB = compute_dB(S, sigma, alpha, beta); auto compute_v = [&S,&area_elem,&toroidal_flux] (const Vector& X) { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Real scal[2]; if (S.Nsurf() == 1) { SCTL_ASSERT(toroidal_flux == true); scal[0] = 1.0 / S.NTor(0); scal[1] = 0; } else if (S.Nsurf() == 2) { if (toroidal_flux == true) { scal[0] = -1.0 / S.NTor(0); scal[1] = 1.0 / S.NTor(1); } else { scal[0] = 1.0 / S.NPol(0); scal[1] = -1.0 / S.NPol(1); } } else { SCTL_ASSERT(false); } Vector v(Nelem * COORD_DIM); Vector dX; ElemBasis::Grad(dX, X); for (Long k = 0; k < S.Nsurf(); k++) { for (Long i_ = 0; i_ < S.NTor(k)*S.NPol(k); i_++) { Long i = S.ElemDsp(k) + i_; for (Long j = 0; j < Nnodes; j++) { Real s = scal[k] / area_elem[i][j]; v[i*COORD_DIM+0][j] = dX[i*COORD_DIM*2+0+(toroidal_flux?1:0)][j] * s; v[i*COORD_DIM+1][j] = dX[i*COORD_DIM*2+2+(toroidal_flux?1:0)][j] * s; v[i*COORD_DIM+2][j] = dX[i*COORD_DIM*2+4+(toroidal_flux?1:0)][j] * s; } } } return v; }; auto compute_AxB = [&S] (const Vector& A, const Vector& B) { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector J(Nelem * COORD_DIM); for (Long i = 0; i < Nelem; i++) { // Set J for (Long j = 0; j < Nnodes; j++) { Tensor a, b; a(0) = A[i*COORD_DIM+0][j]; a(1) = A[i*COORD_DIM+1][j]; a(2) = A[i*COORD_DIM+2][j]; b(0) = B[i*COORD_DIM+0][j]; b(1) = B[i*COORD_DIM+1][j]; b(2) = B[i*COORD_DIM+2][j]; J[i*COORD_DIM+0][j] = a(1) * b(2) - a(2) * b(1); J[i*COORD_DIM+1][j] = a(2) * b(0) - a(0) * b(2); J[i*COORD_DIM+2][j] = a(0) * b(1) - a(1) * b(0); } } return J; }; auto compute_dphi_dnu0 = [&S,&normal,&compute_AxB,&compute_v,&B,compute_grad_adj] () { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector Gv; Vector v = compute_v(S.GetElemList().ElemVector()); EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU); Vector BxGv = compute_AxB(B,Gv); return compute_grad_adj(BxGv)*(-1.0); }; auto compute_dphi_dnu1 = [&S,&normal,&H,&compute_AxB,&compute_v,&B] () { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector Gv; Vector v = compute_v(S.GetElemList().ElemVector()); EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU); Vector BxGv = compute_AxB(B,Gv); Vector n_dot_BxGv = compute_dot_prod(normal,BxGv); Vector dphi_dnu(Nelem); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dphi_dnu[i][j] = n_dot_BxGv[i][j] * 2*H[i][j]; } } return dphi_dnu; }; auto compute_dphi_dnu2 = [&S,&normal,&H,&compute_AxB,&compute_v,&B] () { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector GnxB; Vector nxB = compute_AxB(normal,B); EvalQuadrature(GnxB, S.quadrature_FxU, S, nxB, S.Laplace_FxU); Vector v = compute_v(S.GetElemList().ElemVector()); Vector v_dot_GnxB = compute_dot_prod(v,GnxB); Vector dphi_dnu(Nelem); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dphi_dnu[i][j] = v_dot_GnxB[i][j] * 2*H[i][j]; } } return dphi_dnu; }; auto compute_dphi_dnu3 = [&S,&normal,&area_elem,&H,&compute_AxB,&compute_v,&B] () { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector GnxB; Vector nxB = compute_AxB(normal,B); EvalQuadrature(GnxB, S.quadrature_FxU, S, nxB, S.Laplace_FxU); Vector dGnxB = compute_v(GnxB); Vector v = compute_v(S.GetElemList().ElemVector()); Vector dv_dnu1(Nelem), dv_dnu2(Nelem); { // Set dv_dnu1, dv_dnu2 for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dv_dnu1[i][j] = 0; dv_dnu1[i][j] += -GnxB[i*COORD_DIM+0][j] * v[i*COORD_DIM+0][j] * 2 * H[i][j]; dv_dnu1[i][j] += -GnxB[i*COORD_DIM+1][j] * v[i*COORD_DIM+1][j] * 2 * H[i][j]; dv_dnu1[i][j] += -GnxB[i*COORD_DIM+2][j] * v[i*COORD_DIM+2][j] * 2 * H[i][j]; dv_dnu2[i][j] = 0; dv_dnu2[i][j] += -dGnxB[i*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j]; dv_dnu2[i][j] += -dGnxB[i*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j]; dv_dnu2[i][j] += -dGnxB[i*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j]; } } } return dv_dnu1 + dv_dnu2; }; auto compute_dphi_dnu4 = [&S,&normal,&compute_AxB,&compute_v,&B] () { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector dGnxB; Vector nxB = compute_AxB(normal,B); EvalQuadrature(dGnxB, S.quadrature_FxdU, S, nxB, S.Laplace_FxdU); Vector v = compute_v(S.GetElemList().ElemVector()); Vector dphi_dnu(Nelem); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { Real dphi_dnu_ = 0; dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGnxB[(i*COORD_DIM+0)*COORD_DIM+0][j] * v[i*COORD_DIM+0][j]; dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGnxB[(i*COORD_DIM+0)*COORD_DIM+1][j] * v[i*COORD_DIM+0][j]; dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGnxB[(i*COORD_DIM+0)*COORD_DIM+2][j] * v[i*COORD_DIM+0][j]; dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGnxB[(i*COORD_DIM+1)*COORD_DIM+0][j] * v[i*COORD_DIM+1][j]; dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGnxB[(i*COORD_DIM+1)*COORD_DIM+1][j] * v[i*COORD_DIM+1][j]; dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGnxB[(i*COORD_DIM+1)*COORD_DIM+2][j] * v[i*COORD_DIM+1][j]; dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGnxB[(i*COORD_DIM+2)*COORD_DIM+0][j] * v[i*COORD_DIM+2][j]; dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGnxB[(i*COORD_DIM+2)*COORD_DIM+1][j] * v[i*COORD_DIM+2][j]; dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGnxB[(i*COORD_DIM+2)*COORD_DIM+2][j] * v[i*COORD_DIM+2][j]; dphi_dnu[i][j] = dphi_dnu_; } } return dphi_dnu; }; auto compute_dphi_dnu5 = [&S,&normal,&compute_AxB,&compute_v,&B] () { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector nxB = compute_AxB(normal,B); Vector dGv; Vector v = compute_v(S.GetElemList().ElemVector()); EvalQuadrature(dGv, S.quadrature_FxdU, S, v, S.Laplace_FxdU); Vector dphi_dnu(Nelem); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { Real dphi_dnu_ = 0; dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGv[(i*COORD_DIM+0)*COORD_DIM+0][j] * nxB[i*COORD_DIM+0][j]; dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGv[(i*COORD_DIM+0)*COORD_DIM+1][j] * nxB[i*COORD_DIM+0][j]; dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGv[(i*COORD_DIM+0)*COORD_DIM+2][j] * nxB[i*COORD_DIM+0][j]; dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGv[(i*COORD_DIM+1)*COORD_DIM+0][j] * nxB[i*COORD_DIM+1][j]; dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGv[(i*COORD_DIM+1)*COORD_DIM+1][j] * nxB[i*COORD_DIM+1][j]; dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGv[(i*COORD_DIM+1)*COORD_DIM+2][j] * nxB[i*COORD_DIM+1][j]; dphi_dnu_ += -normal[i*COORD_DIM+0][j] * dGv[(i*COORD_DIM+2)*COORD_DIM+0][j] * nxB[i*COORD_DIM+2][j]; dphi_dnu_ += -normal[i*COORD_DIM+1][j] * dGv[(i*COORD_DIM+2)*COORD_DIM+1][j] * nxB[i*COORD_DIM+2][j]; dphi_dnu_ += -normal[i*COORD_DIM+2][j] * dGv[(i*COORD_DIM+2)*COORD_DIM+2][j] * nxB[i*COORD_DIM+2][j]; dphi_dnu[i][j] = dphi_dnu_; } } return dphi_dnu; }; auto compute_dphi_dnu6 = [&S,&normal,&compute_AxB,&compute_v,&gradB] () { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector Gv; Vector v = compute_v(S.GetElemList().ElemVector()); EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU); Vector nxGv = compute_AxB(Gv,normal); Vector dphi_dnu(Nelem); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { Real dphi_dnu_ = 0; dphi_dnu_ += -nxGv[i*COORD_DIM+0][j] * gradB[(i*COORD_DIM+0)*COORD_DIM+0][j] * normal[i*COORD_DIM+0][j]; dphi_dnu_ += -nxGv[i*COORD_DIM+1][j] * gradB[(i*COORD_DIM+0)*COORD_DIM+1][j] * normal[i*COORD_DIM+0][j]; dphi_dnu_ += -nxGv[i*COORD_DIM+2][j] * gradB[(i*COORD_DIM+0)*COORD_DIM+2][j] * normal[i*COORD_DIM+0][j]; dphi_dnu_ += -nxGv[i*COORD_DIM+0][j] * gradB[(i*COORD_DIM+1)*COORD_DIM+0][j] * normal[i*COORD_DIM+1][j]; dphi_dnu_ += -nxGv[i*COORD_DIM+1][j] * gradB[(i*COORD_DIM+1)*COORD_DIM+1][j] * normal[i*COORD_DIM+1][j]; dphi_dnu_ += -nxGv[i*COORD_DIM+2][j] * gradB[(i*COORD_DIM+1)*COORD_DIM+2][j] * normal[i*COORD_DIM+1][j]; dphi_dnu_ += -nxGv[i*COORD_DIM+0][j] * gradB[(i*COORD_DIM+2)*COORD_DIM+0][j] * normal[i*COORD_DIM+2][j]; dphi_dnu_ += -nxGv[i*COORD_DIM+1][j] * gradB[(i*COORD_DIM+2)*COORD_DIM+1][j] * normal[i*COORD_DIM+2][j]; dphi_dnu_ += -nxGv[i*COORD_DIM+2][j] * gradB[(i*COORD_DIM+2)*COORD_DIM+2][j] * normal[i*COORD_DIM+2][j]; dphi_dnu[i][j] = dphi_dnu_; } } return dphi_dnu; }; auto compute_dphi_dnu7 = [&S,&normal,&H,&compute_AxB,&compute_v,&sigma] () { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector Gv; Vector v = compute_v(S.GetElemList().ElemVector()); EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU); Vector nxGv = compute_AxB(Gv,normal); Vector dphi_dnu(Nelem); EvalQuadrature(dphi_dnu, S.quadrature_dUxF, S, nxGv, S.Laplace_dUxF); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dphi_dnu[i][j] *= -2*H[i][j] * sigma[i][j]; } } return dphi_dnu; }; auto compute_dphi_dnu8 = [&S,&normal,&H,&compute_AxB,&compute_v,&sigma] () { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector Gv; Vector v = compute_v(S.GetElemList().ElemVector()); EvalQuadrature(Gv, S.quadrature_FxU, S, v, S.Laplace_FxU); Vector nxGv = compute_AxB(Gv,normal); Vector dphi_dnu(Nelem); Vector nablaDt_nxGv; EvalQuadrature(nablaDt_nxGv, S.quadrature_dUxD, S, nxGv, S.Laplace_dUxD); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { dphi_dnu[i][j] = 0; dphi_dnu[i][j] += sigma[i][j] * nablaDt_nxGv[i*COORD_DIM+0][j]*normal[i*COORD_DIM+0][j]; dphi_dnu[i][j] += sigma[i][j] * nablaDt_nxGv[i*COORD_DIM+1][j]*normal[i*COORD_DIM+1][j]; dphi_dnu[i][j] += sigma[i][j] * nablaDt_nxGv[i*COORD_DIM+2][j]*normal[i*COORD_DIM+2][j]; } } return dphi_dnu; }; auto dphi_dnu0 = compute_dphi_dnu0(); auto dphi_dnu1 = compute_dphi_dnu1(); auto dphi_dnu2 = compute_dphi_dnu2(); auto dphi_dnu3 = compute_dphi_dnu3(); auto dphi_dnu4 = compute_dphi_dnu4(); auto dphi_dnu5 = compute_dphi_dnu5(); auto dphi_dnu6 = compute_dphi_dnu6(); auto dphi_dnu7 = compute_dphi_dnu7(); auto dphi_dnu8 = compute_dphi_dnu8(); return (dphi_dnu0+dphi_dnu1+dphi_dnu2+dphi_dnu3+dphi_dnu4+dphi_dnu5+dphi_dnu6+dphi_dnu7+dphi_dnu8); }; { // Set dg_dnu -= dg_dsigma invA dA_dnu sigma dg_dnu -= compute_u_dAdnu_v_0(dg_dsigma_invA, sigma, alpha, beta); if (S.Nsurf() >= 1) dg_dnu -= compute_u_dAdnu_v_1(sigma, alpha, beta, true) * dg_dsigma_invA[Nelem*Nnodes+0]; if (S.Nsurf() >= 2) dg_dnu -= compute_u_dAdnu_v_1(sigma, alpha, beta, false) * dg_dsigma_invA[Nelem*Nnodes+1]; } return dg_dnu; }; Vector dgdnu; { // Set dgdnu dgdnu.ReInit(S_.NElem()); dgdnu = 0; for (Long i = 0; i < S_.Nsurf(); i++) { const Long elem_dsp = (i==0 ? 0 : S_.ElemDsp(i-1)); const Long Nnodes = ElemBasis::Size(); auto dgdnu_ = compute_gradient(Svec[i]); if (0) { // Write VTU VTUData vtu; vtu.AddElems(Svec[i].GetElemList(), dgdnu_, ORDER); vtu.WriteVTK("dgdnu-"+std::to_string(i), comm); } for (Long j = 0; j < (i==0?0:Svec[i].NTor(0)*Svec[i].NPol(0)); j++) { for (Long k = 0; k < Nnodes; k++) { dgdnu[elem_dsp+j][k] -= dgdnu_[j][k]; } } for (Long j = (i==0?0:Svec[i].NTor(0)*Svec[i].NPol(0)); j < dgdnu_.Dim(); j++) { for (Long k = 0; k < Nnodes; k++) { dgdnu[elem_dsp+j][k] += dgdnu_[j][k]; } } } } return dgdnu; } static Vector compute_pressure_jump(const Stellarator& S_, const Vector& pressure, const Vector& flux_tor_, const Vector& flux_pol_, Real* g_ptr = nullptr) { Comm comm = Comm::World(); Vector> Svec(S_.Nsurf()); for (Long i = 0; i < S_.Nsurf(); i++) { // Set Svec[i] (quadratures, B) const Long elem_dsp = (i==0 ? 0 : S_.ElemDsp(i-1)); const Long Nnodes = ElemBasis::Size(); Stellarator& S = Svec[i]; if (i == 0) { // Init S Vector NtNp; NtNp.PushBack(S_.NTor(i)); NtNp.PushBack(S_.NPol(i)); S = Stellarator(NtNp); } else { Vector NtNp; NtNp.PushBack(S_.NTor(i-1)); NtNp.PushBack(S_.NPol(i-1)); NtNp.PushBack(S_.NTor(i)); NtNp.PushBack(S_.NPol(i)); S = Stellarator(NtNp); } for (Long j = 0; j < S.NElem(); j++) { // Set S coordinates for (Long k = 0; k < Nnodes; k++) { S.Elem(j,0)[k] = S_.Elem(elem_dsp+j,0)[k]; S.Elem(j,1)[k] = S_.Elem(elem_dsp+j,1)[k]; S.Elem(j,2)[k] = S_.Elem(elem_dsp+j,2)[k]; } } SetupQuadrature(S.quadrature_BS , S, S.BiotSavart , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER)); SetupQuadrature(S.quadrature_FxU , S, S.Laplace_FxU , order_singular, order_direct, -1.0, comm); SetupQuadrature(S.quadrature_FxdU , S, S.Laplace_FxdU , order_singular, order_direct, -1.0, comm); { // Set Bt0, Bp0, dBt0, dBp0 Vector Jt, Jp; compute_harmonic_vector_potentials(Jt, Jp, S); EvalQuadrature(S.Bt0 , S.quadrature_BS , S, Jp, S.BiotSavart); EvalQuadrature(S.Bp0 , S.quadrature_BS , S, Jt, S.BiotSavart); } compute_invA(S.sigma, S.alpha, S.beta, S, flux_tor_[i], flux_pol_[i], comm); S.B = compute_B(S, S.sigma, S.alpha, S.beta); } if (g_ptr != nullptr) g_ptr[0] = compute_g(Svec, pressure); return compute_pressure_jump(Svec, pressure); } static void test() { Comm comm = Comm::World(); Profile::Enable(true); Long Nsurf = 2; Stellarator S; Vector flux_tor(Nsurf), flux_pol(Nsurf), pressure(Nsurf); { // Init S, flux_tor, flux_pol, pressure Vector NtNp; for (Long i = 0; i < Nsurf; i++) { NtNp.PushBack(30); NtNp.PushBack(4); } S = Stellarator(NtNp); flux_tor = 1; flux_pol = 1; pressure = 0; //flux_tor[0] = 1; //0.791881512; //flux_tor[1] = 1; //flux_pol[0] = 0; //flux_pol[1] = 0; //pressure[0] = 0; //pressure[1] = 0; } { // find equilibrium flux surfaces { //auto filter = [](const Stellarator& S, Vector& f) { // auto cheb2grid = [] (const Vector& X, Long Mt, Long Mp, Long Nt, Long Np) { // const Long dof = X.Dim() / (Mt * Mp); // SCTL_ASSERT(X.Dim() == Mt * Mp *dof); // Vector Xf(dof*Nt*Np); Xf = 0; // const Long Nnodes = ElemBasis::Size(); // const Matrix& Mnodes = Basis::Nodes(); // for (Long t = 0; t < Nt; t++) { // for (Long p = 0; p < Np; p++) { // Real theta = t / (Real)Nt; // Real phi = p / (Real)Np; // Long i = (Long)(theta * Mt); // Long j = (Long)(phi * Mp); // Real x = theta * Mt - i; // Real y = phi * Mp - j; // Long elem_idx = i * Mp + j; // Vector Interp0(ORDER); // Vector Interp1(ORDER); // { // Set Interp0, Interp1 // auto node = [&Mnodes] (Long i) { // return Mnodes[0][i]; // }; // for (Long i = 0; i < ORDER; i++) { // Real wt_x = 1, wt_y = 1; // for (Long j = 0; j < ORDER; j++) { // if (j != i) { // wt_x *= (x - node(j)) / (node(i) - node(j)); // wt_y *= (y - node(j)) / (node(i) - node(j)); // } // Interp0[i] = wt_x; // Interp1[i] = wt_y; // } // } // } // for (Long ii = 0; ii < ORDER; ii++) { // for (Long jj = 0; jj < ORDER; jj++) { // Long node_idx = jj * ORDER + ii; // for (Long k = 0; k < dof; k++) { // Xf[(k*Nt+t)*Np+p] += X[elem_idx*dof+k][node_idx] * Interp0[ii] * Interp1[jj]; // } // } // } // } // } // return Xf; // }; // auto grid2cheb = [] (const Vector& Xf, Long Nt, Long Np, Long Mt, Long Mp) { // Long dof = Xf.Dim() / (Nt*Np); // SCTL_ASSERT(Xf.Dim() == dof*Nt*Np); // Vector X(Mt*Mp*dof); // constexpr Integer INTERP_ORDER = 12; // for (Long tt = 0; tt < Mt; tt++) { // for (Long pp = 0; pp < Mp; pp++) { // for (Long t = 0; t < ORDER; t++) { // for (Long p = 0; p < ORDER; p++) { // Matrix Mnodes = Basis::Nodes(); // Real theta = (tt + Mnodes[0][t]) / Mt; // Real phi = (pp + Mnodes[0][p]) / Mp; // Long i = (Long)(theta * Nt); // Long j = (Long)(phi * Np); // Real x = theta * Nt - i; // Real y = phi * Np - j; // Vector Interp0(INTERP_ORDER); // Vector Interp1(INTERP_ORDER); // { // Set Interp0, Interp1 // auto node = [] (Long i) { // return (Real)i - (INTERP_ORDER-1)/2; // }; // for (Long i = 0; i < INTERP_ORDER; i++) { // Real wt_x = 1, wt_y = 1; // for (Long j = 0; j < INTERP_ORDER; j++) { // if (j != i) { // wt_x *= (x - node(j)) / (node(i) - node(j)); // wt_y *= (y - node(j)) / (node(i) - node(j)); // } // Interp0[i] = wt_x; // Interp1[i] = wt_y; // } // } // } // for (Long k = 0; k < dof; k++) { // Real X0 = 0; // for (Long ii = 0; ii < INTERP_ORDER; ii++) { // for (Long jj = 0; jj < INTERP_ORDER; jj++) { // Long idx_i = (i + ii-(INTERP_ORDER-1)/2 + Nt) % Nt; // Long idx_j = (j + jj-(INTERP_ORDER-1)/2 + Np) % Np; // X0 += Interp0[ii] * Interp1[jj] * Xf[(k*Nt+idx_i)*Np+idx_j]; // } // } // Long elem_idx = tt * Mp + pp; // Long node_idx = p * ORDER + t; // X[elem_idx*dof+k][node_idx] = X0; // } // } // } // } // } // return X; // }; // Long dof = f.Dim() / S.NElem(); // SCTL_ASSERT(f.Dim() == S.NElem() * dof); // for (Long i = 0; i < S.Nsurf(); i++) { // const Long Mt = S.NTor(i); // const Long Mp = S.NPol(i); // const Long Nelem = Mt * Mp; // const Long offset = S.ElemDsp(i); // const Long Nt = Mt * ORDER / 5; // const Long Np = Mp * ORDER / 5; // Vector f_(Nelem*dof, f.begin() + offset*dof, false); // Vector f_fourier = cheb2grid(f_, Mt, Mp, Nt, Np); // f_ = grid2cheb(f_fourier, Nt, Np, Mt, Mp); // } //}; } auto filter = [](const Stellarator& S, const Comm& comm, Vector& f, Real sigma) { auto cheb2grid = [] (const Vector& X, Long Mt, Long Mp, Long Nt, Long Np) { const Long dof = X.Dim() / (Mt * Mp); SCTL_ASSERT(X.Dim() == Mt * Mp *dof); Vector Xf(dof*Nt*Np); Xf = 0; const Long Nnodes = ElemBasis::Size(); const Matrix& Mnodes = Basis::Nodes(); for (Long t = 0; t < Nt; t++) { for (Long p = 0; p < Np; p++) { Real theta = t / (Real)Nt; Real phi = p / (Real)Np; Long i = (Long)(theta * Mt); Long j = (Long)(phi * Mp); Real x = theta * Mt - i; Real y = phi * Mp - j; Long elem_idx = i * Mp + j; Vector Interp0(ORDER); Vector Interp1(ORDER); { // Set Interp0, Interp1 auto node = [&Mnodes] (Long i) { return Mnodes[0][i]; }; for (Long i = 0; i < ORDER; i++) { Real wt_x = 1, wt_y = 1; for (Long j = 0; j < ORDER; j++) { if (j != i) { wt_x *= (x - node(j)) / (node(i) - node(j)); wt_y *= (y - node(j)) / (node(i) - node(j)); } Interp0[i] = wt_x; Interp1[i] = wt_y; } } } for (Long ii = 0; ii < ORDER; ii++) { for (Long jj = 0; jj < ORDER; jj++) { Long node_idx = jj * ORDER + ii; for (Long k = 0; k < dof; k++) { Xf[(k*Nt+t)*Np+p] += X[elem_idx*dof+k][node_idx] * Interp0[ii] * Interp1[jj]; } } } } } return Xf; }; auto grid2cheb = [] (const Vector& Xf, Long Nt, Long Np, Long Mt, Long Mp) { Long dof = Xf.Dim() / (Nt*Np); SCTL_ASSERT(Xf.Dim() == dof*Nt*Np); Vector X(Mt*Mp*dof); constexpr Integer INTERP_ORDER = 12; for (Long tt = 0; tt < Mt; tt++) { for (Long pp = 0; pp < Mp; pp++) { for (Long t = 0; t < ORDER; t++) { for (Long p = 0; p < ORDER; p++) { Matrix Mnodes = Basis::Nodes(); Real theta = (tt + Mnodes[0][t]) / Mt; Real phi = (pp + Mnodes[0][p]) / Mp; Long i = (Long)(theta * Nt); Long j = (Long)(phi * Np); Real x = theta * Nt - i; Real y = phi * Np - j; Vector Interp0(INTERP_ORDER); Vector Interp1(INTERP_ORDER); { // Set Interp0, Interp1 auto node = [] (Long i) { return (Real)i - (INTERP_ORDER-1)/2; }; for (Long i = 0; i < INTERP_ORDER; i++) { Real wt_x = 1, wt_y = 1; for (Long j = 0; j < INTERP_ORDER; j++) { if (j != i) { wt_x *= (x - node(j)) / (node(i) - node(j)); wt_y *= (y - node(j)) / (node(i) - node(j)); } Interp0[i] = wt_x; Interp1[i] = wt_y; } } } for (Long k = 0; k < dof; k++) { Real X0 = 0; for (Long ii = 0; ii < INTERP_ORDER; ii++) { for (Long jj = 0; jj < INTERP_ORDER; jj++) { Long idx_i = (i + ii-(INTERP_ORDER-1)/2 + Nt) % Nt; Long idx_j = (j + jj-(INTERP_ORDER-1)/2 + Np) % Np; X0 += Interp0[ii] * Interp1[jj] * Xf[(k*Nt+idx_i)*Np+idx_j]; } } Long elem_idx = tt * Mp + pp; Long node_idx = p * ORDER + t; X[elem_idx*dof+k][node_idx] = X0; } } } } } return X; }; auto fourier_filter = [](sctl::Vector& X, long Nt_, long Np_, Real sigma, const Comm& comm) { long dof = X.Dim() / (Nt_ * Np_); SCTL_ASSERT(X.Dim() == dof * Nt_ * Np_); sctl::FFT fft_r2c, fft_c2r; sctl::StaticArray fft_dim = {Nt_, Np_}; fft_r2c.Setup(sctl::FFT_Type::R2C, 1, sctl::Vector(2, fft_dim, false), omp_get_max_threads()); fft_c2r.Setup(sctl::FFT_Type::C2R, 1, sctl::Vector(2, fft_dim, false), omp_get_max_threads()); long Nt = Nt_; long Np = fft_r2c.Dim(1) / (Nt * 2); SCTL_ASSERT(fft_r2c.Dim(1) == Nt * Np * 2); //auto filter_fn = [](Real x2, Real sigma) {return exp(-x2/(2*sigma*sigma));}; auto filter_fn = [](Real x2, Real sigma) {return (x2 normal, gradX; biest::SurfaceOp op(comm, Nt_, Np_); sctl::Vector coeff(fft_r2c.Dim(1)); for (long k = 0; k < dof; k++) { sctl::Vector X_(Nt_*Np_, X.begin() + k*Nt_*Np_, false); fft_r2c.Execute(X_, coeff); for (long t = 0; t < Nt; t++) { for (long p = 0; p < Np; p++) { Real tt = (t - (t > Nt / 2 ? Nt : 0)) / (Real)(Nt / 2); Real pp = p / (Real)Np; Real f = filter_fn(tt*tt+pp*pp, sigma); coeff[(t * Np + p) * 2 + 0] *= f; coeff[(t * Np + p) * 2 + 1] *= f; } } fft_c2r.Execute(coeff, X_); } }; Long dof = f.Dim() / S.NElem(); SCTL_ASSERT(f.Dim() == S.NElem() * dof); for (Long i = 0; i < S.Nsurf(); i++) { const Long Mt = S.NTor(i); const Long Mp = S.NPol(i); const Long Nelem = Mt * Mp; const Long offset = S.ElemDsp(i); const Long Nt = Mt * ORDER * 4; const Long Np = Mp * ORDER * 4; Vector f_(Nelem*dof, f.begin() + offset*dof, false); Vector f_fourier = cheb2grid(f_, Mt, Mp, Nt, Np); fourier_filter(f_fourier, Nt, Np, 0.25 * sigma, comm); f_ = grid2cheb(f_fourier, Nt, Np, Mt, Mp); } }; Long iter = 0; Real dt = 0.1; while (1) { // time-step Vector dgdnu = compute_gradient(S, pressure, flux_tor, flux_pol)*(-1); //Vector dgdnu = compute_pressure_jump(S, pressure, flux_tor, flux_pol)*(-1); Vector dXdt(dgdnu.Dim()*COORD_DIM); { // Set dXdt dXdt = 0; const Long Nnodes = ElemBasis::Size(); Vector normal, area_elem; compute_norm_area_elem(S, normal, area_elem); for (Long i = 0; i < S.ElemDsp(S.Nsurf()-1); i++) { for (Long j = 0; j < Nnodes; j++) { dXdt[i*COORD_DIM+0][j] = normal[i*COORD_DIM+0][j] * dgdnu[i][j]; dXdt[i*COORD_DIM+1][j] = normal[i*COORD_DIM+1][j] * dgdnu[i][j]; dXdt[i*COORD_DIM+2][j] = normal[i*COORD_DIM+2][j] * dgdnu[i][j]; } } filter(S, comm, dXdt, 0.1); } { // Update dt const Long Nelem = S.NElem(); Stellarator S0 = S, S1 = S, S2 = S; for (Long i = 0; i < S.NElem(); i++) { S0.Elem(i, 0) += dXdt[i*COORD_DIM+0] * 0.0 * dt; S0.Elem(i, 1) += dXdt[i*COORD_DIM+1] * 0.0 * dt; S0.Elem(i, 2) += dXdt[i*COORD_DIM+2] * 0.0 * dt; S1.Elem(i, 0) += dXdt[i*COORD_DIM+0] * 0.5 * dt; S1.Elem(i, 1) += dXdt[i*COORD_DIM+1] * 0.5 * dt; S1.Elem(i, 2) += dXdt[i*COORD_DIM+2] * 0.5 * dt; S2.Elem(i, 0) += dXdt[i*COORD_DIM+0] * 1.0 * dt; S2.Elem(i, 1) += dXdt[i*COORD_DIM+1] * 1.0 * dt; S2.Elem(i, 2) += dXdt[i*COORD_DIM+2] * 1.0 * dt; } Real g0, g1, g2; compute_pressure_jump(S0, pressure, flux_tor, flux_pol, &g0); compute_pressure_jump(S1, pressure, flux_tor, flux_pol, &g1); compute_pressure_jump(S2, pressure, flux_tor, flux_pol, &g2); { // Calculate optimal step size dt Real a = 2*g0 - 4*g1 + 2*g2; Real b =-3*g0 + 4*g1 - g2; Real c = g0; Real s = -b/(2*a); dt *= s; Real g_ = a*s*s + b*s + c; std::cout<<"g = "< pressure_jump = compute_pressure_jump(S, pressure, flux_tor, flux_pol); VTUData vtu; vtu.AddElems(S.GetElemList(), pressure_jump, ORDER); vtu.WriteVTK("pressure_jump"+std::to_string(iter), comm); } { // Update S <-- filter(S + dXdt * dt) const Long Nelem = S.NElem(); Vector X(Nelem*COORD_DIM); for (Long i = 0; i < S.NElem(); i++) { X[i*COORD_DIM+0] = S.Elem(i, 0) + dXdt[i*COORD_DIM+0] * dt; X[i*COORD_DIM+1] = S.Elem(i, 1) + dXdt[i*COORD_DIM+1] * dt; X[i*COORD_DIM+2] = S.Elem(i, 2) + dXdt[i*COORD_DIM+2] * dt; } filter(S, comm, X, 0.3); for (Long i = 0; i < S.NElem(); i++) { S.Elem(i, 0) = X[i*COORD_DIM+0]; S.Elem(i, 1) = X[i*COORD_DIM+1]; S.Elem(i, 2) = X[i*COORD_DIM+2]; } } iter++; } return; } { // Verify using finite difference approximation Vector dgdnu = compute_gradient(S, pressure, flux_tor, flux_pol); { // Write VTU VTUData vtu; vtu.AddElems(S.GetElemList(), dgdnu, ORDER); vtu.WriteVTK("dgdnu", comm); } Real eps = 1e-4; const Long Nnodes = ElemBasis::Size(); Vector normal, area_elem; compute_norm_area_elem(S, normal, area_elem); Vector nu = area_elem; for (Long i = S.ElemDsp(S.Nsurf()-1); i < S.NElem(); i++) nu[i] = 0; Stellarator S0 = S, S1 = S; for (Long i = 0; i < S.NElem(); i++) { for (Long j = 0; j < Nnodes; j++) { S0.Elem(i, 0)[j] -= 0.5 * eps * normal[i*COORD_DIM+0][j] * nu[i][j]; S0.Elem(i, 1)[j] -= 0.5 * eps * normal[i*COORD_DIM+1][j] * nu[i][j]; S0.Elem(i, 2)[j] -= 0.5 * eps * normal[i*COORD_DIM+2][j] * nu[i][j]; S1.Elem(i, 0)[j] += 0.5 * eps * normal[i*COORD_DIM+0][j] * nu[i][j]; S1.Elem(i, 1)[j] += 0.5 * eps * normal[i*COORD_DIM+1][j] * nu[i][j]; S1.Elem(i, 2)[j] += 0.5 * eps * normal[i*COORD_DIM+2][j] * nu[i][j]; } } Real g0, g1; compute_pressure_jump(S0, pressure, flux_tor, flux_pol, &g0); compute_pressure_jump(S1, pressure, flux_tor, flux_pol, &g1); std::cout<<"g0 = "<& nu, Real eps, Real flux_tor, Real flux_pol, const Vector& pressure) { const Long Nelem = S.NElem(); const Long Nnodes = ElemBasis::Size(); Vector normal, area_elem; compute_norm_area_elem(S, normal, area_elem); Vector X_orig(Nelem*COORD_DIM); for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { X_orig[i*COORD_DIM+0][j] = S.Elem(i,0)[j]; X_orig[i*COORD_DIM+1][j] = S.Elem(i,1)[j]; X_orig[i*COORD_DIM+2][j] = S.Elem(i,2)[j]; S.Elem(i,0)[j] += eps*nu[i][j] * normal[i*COORD_DIM+0][j]; S.Elem(i,1)[j] += eps*nu[i][j] * normal[i*COORD_DIM+1][j]; S.Elem(i,2)[j] += eps*nu[i][j] * normal[i*COORD_DIM+2][j]; } } ///////////////////////////////////////////////////////////////////////////////////////////////////////////// SetupQuadrature(S.quadrature_BS , S, S.BiotSavart , order_singular, order_direct, -1.0, comm, -0.01 * pow<-2,Real>(ORDER)); SetupQuadrature(S.quadrature_FxU , S, S.Laplace_FxU , order_singular, order_direct, -1.0, comm); SetupQuadrature(S.quadrature_FxdU, S, S.Laplace_FxdU, order_singular, order_direct, -1.0, comm); Vector Jt, Jp; compute_harmonic_vector_potentials(Jt, Jp, S); EvalQuadrature(Bt0, S.quadrature_BS, S, Jp, S.BiotSavart); EvalQuadrature(Bp0, S.quadrature_BS, S, Jt, S.BiotSavart); Real alpha, beta; Vector sigma; compute_invA(sigma, alpha, beta, flux_tor, flux_pol); Vector B = compute_B(sigma, alpha, beta); compute_norm_area_elem(S, normal, area_elem); Real g = compute_inner_prod(area_elem, compute_gvec(S,B,pressure), area_elem*0+1); ///////////////////////////////////////////////////////////////////////////////////////////////////////////// for (Long i = 0; i < Nelem; i++) { for (Long j = 0; j < Nnodes; j++) { S.Elem(i,0)[j] = X_orig[i*COORD_DIM+0][j]; S.Elem(i,1)[j] = X_orig[i*COORD_DIM+1][j]; S.Elem(i,2)[j] = X_orig[i*COORD_DIM+2][j]; } } return g; }; Vector normal, area_elem; compute_norm_area_elem(S, normal, area_elem); const Long Nelem = S.NElem(); { Vector nu(Nelem); nu = area_elem; Real eps = 1e-4; Real g0 = compute_g(nu,-eps, flux_tor, flux_pol, pressure); Real g1 = compute_g(nu,eps, flux_tor, flux_pol, pressure); std::cout<<"g = "< M(S.NtNp_[0]*ORDER, S.NtNp_[1]*ORDER); // for (Long tt = 0; tt < S.NtNp_[0]; tt++) { // for (Long pp = 0; pp < S.NtNp_[1]; pp++) { // for (Long t = 0; t < ORDER; t++) { // for (Long p = 0; p < ORDER; p++) { // Long elem_idx = tt * S.NtNp_[1] + pp; // Long node_idx = p * ORDER + t; // M[tt*ORDER+t][pp*ORDER+p] = dg_dnu[elem_idx][node_idx]; // } // } // } // } // M.Write("dg_dnu.mat"); //} //if (0) { // filter dg_dnu and write VTU // const Long Nelem = S.NElem(); // const Long Nnodes = ElemBasis::Size(); // const Integer INTERP_ORDER = 12; // Long Nt = S.NtNp_[0]*ORDER/5, Np = S.NtNp_[1]*ORDER/5; // Matrix M(Nt, Np); M = 0; // const auto& quad_wts = ElemBasis::QuadWts(); // const Matrix& Mnodes = Basis::Nodes(); // for (Long tt = 0; tt < S.NtNp_[0]; tt++) { // for (Long pp = 0; pp < S.NtNp_[1]; pp++) { // for (Long t = 0; t < ORDER; t++) { // for (Long p = 0; p < ORDER; p++) { // Real theta = (tt + Mnodes[0][t]) / S.NtNp_[0]; // Real phi = (pp + Mnodes[0][p]) / S.NtNp_[1]; // Long i = (Long)(theta * Nt); // Long j = (Long)(phi * Np); // Real x = theta * Nt - i; // Real y = phi * Np - j; // Long elem_idx = tt * S.NtNp_[1] + pp; // Long node_idx = p * ORDER + t; // Vector Interp0(INTERP_ORDER); // Vector Interp1(INTERP_ORDER); // { // Set Interp0, Interp1 // auto node = [] (Long i) { // return (Real)i - (INTERP_ORDER-1)/2; // }; // for (Long i = 0; i < INTERP_ORDER; i++) { // Real wt_x = 1, wt_y = 1; // for (Long j = 0; j < INTERP_ORDER; j++) { // if (j != i) { // wt_x *= (x - node(j)) / (node(i) - node(j)); // wt_y *= (y - node(j)) / (node(i) - node(j)); // } // Interp0[i] = wt_x; // Interp1[i] = wt_y; // } // } // } // for (Long ii = 0; ii < INTERP_ORDER; ii++) { // for (Long jj = 0; jj < INTERP_ORDER; jj++) { // Long idx_i = (i + ii-(INTERP_ORDER-1)/2 + Nt) % Nt; // Long idx_j = (j + jj-(INTERP_ORDER-1)/2 + Np) % Np; // M[idx_i][idx_j] += dg_dnu[elem_idx][node_idx] * quad_wts[node_idx] * Interp0[ii] * Interp1[jj] / (S.NtNp_[0] * S.NtNp_[1]) * (Nt * Np); // } // } // } // } // } // } // Vector f(Nelem); // for (Long tt = 0; tt < S.NtNp_[0]; tt++) { // for (Long pp = 0; pp < S.NtNp_[1]; pp++) { // for (Long t = 0; t < ORDER; t++) { // for (Long p = 0; p < ORDER; p++) { // Matrix Mnodes = Basis::Nodes(); // Real theta = (tt + Mnodes[0][t]) / S.NtNp_[0]; // Real phi = (pp + Mnodes[0][p]) / S.NtNp_[1]; // Long i = (Long)(theta * Nt); // Long j = (Long)(phi * Np); // Real x = theta * Nt - i; // Real y = phi * Np - j; // Vector Interp0(INTERP_ORDER); // Vector Interp1(INTERP_ORDER); // { // Set Interp0, Interp1 // auto node = [] (Long i) { // return (Real)i - (INTERP_ORDER-1)/2; // }; // for (Long i = 0; i < INTERP_ORDER; i++) { // Real wt_x = 1, wt_y = 1; // for (Long j = 0; j < INTERP_ORDER; j++) { // if (j != i) { // wt_x *= (x - node(j)) / (node(i) - node(j)); // wt_y *= (y - node(j)) / (node(i) - node(j)); // } // Interp0[i] = wt_x; // Interp1[i] = wt_y; // } // } // } // Real f0 = 0; // for (Long ii = 0; ii < INTERP_ORDER; ii++) { // for (Long jj = 0; jj < INTERP_ORDER; jj++) { // Long idx_i = (i + ii-(INTERP_ORDER-1)/2 + Nt) % Nt; // Long idx_j = (j + jj-(INTERP_ORDER-1)/2 + Np) % Np; // f0 += Interp0[ii] * Interp1[jj] * M[idx_i][idx_j]; // } // } // Long elem_idx = tt * S.NtNp_[1] + pp; // Long node_idx = p * ORDER + t; // f[elem_idx][node_idx] = f0; // } // } // } // } // { // Write VTU // VTUData vtu; // vtu.AddElems(S.GetElemList(), f, ORDER); // vtu.WriteVTK("dg_dnu_filtered", comm); // } // dg_dnu = f; //} } static void FlipNormal(Vector& v) { for (Long i = 0; i < v.Dim(); i++) { const auto elem = v[i]; for (Long j0 = 0; j0 < ORDER; j0++) { for (Long j1 = 0; j1 < ORDER; j1++) { v[i][j0*ORDER+j1] = elem[j0*ORDER+(ORDER-j1-1)]; } } } } template static void SetupQuadrature(Quadrature& quadrature, const Stellarator& S, const Kernel& kernel, Integer order_singular, Integer order_direct, Real period_length, const Comm& comm, Real Rqbx = 0) { if (S.Nsurf() == 2) { Long Nelem0 = S.NTor(0)*S.NPol(0); ElemList elem_lst = S.GetElemList(); { // Update elem_lst Vector X = elem_lst.ElemVector(); Vector X0(Nelem0*COORD_DIM, X.begin(), false); FlipNormal(X0); elem_lst.ReInit(X); } quadrature.template Setup(elem_lst, kernel, order_singular, order_direct, period_length, comm, Rqbx); } else { quadrature.template Setup(S.GetElemList(), kernel, order_singular, order_direct, period_length, comm, Rqbx); } } template static void EvalQuadrature(Vector& potential, const Quadrature& quadrature, const Stellarator& S, const Vector& density, const Kernel& kernel) { if (S.Nsurf() == 2) { Long Nelem0 = S.NTor(0)*S.NPol(0); Vector potential_, density_ = density; ElemList elem_lst = S.GetElemList(); { // Update elem_lst Vector X = elem_lst.ElemVector(); Vector X0(Nelem0*COORD_DIM, X.begin(), false); FlipNormal(X0); elem_lst.ReInit(X); } { // Update density_ Long dof = density_.Dim() / S.NElem(); Vector density0(Nelem0*dof, density_.begin(), false); FlipNormal(density0); } quadrature.Eval(potential_, elem_lst, density_, kernel); { // Update potential_ Long dof = potential_.Dim() / S.NElem(); Vector potential0(Nelem0*dof, potential_.begin(), false); FlipNormal(potential0); } potential = potential_; } else { quadrature.Eval(potential, S.GetElemList(), density, kernel); } } void InitSurf(Long l, Long Nsurf) { const auto& nodes = ElemBasis::Nodes(); const Long Nt = NTor(l); const Long Np = NPol(l); for (Long i = 0; i < Nt; i++) { for (Long j = 0; j < Np; j++) { for (Long k = 0; k < ElemBasis::Size(); k++) { Real theta = (i + nodes[0][k]) * 2*const_pi()/Nt; Real phi = (j + nodes[1][k]) * 2*const_pi()/Np; Real X,Y,Z; SurfGeom(X,Y,Z,theta,phi, (2.0+l)/(1.0+Nsurf)); Elem(ElemIdx(l,i,j),0)[k] = X; Elem(ElemIdx(l,i,j),1)[k] = Y; Elem(ElemIdx(l,i,j),2)[k] = Z; } } } } static void SurfGeom(Real& X, Real& Y, Real& Z, Real theta, Real phi, Real s) { sctl::Integer Nperiod = 5; #if 0 Real Aspect_ratio = 10.27932548522949; Real coeffmat[21][21] = { 0.00000478813217, 0.00000000000000, 0.00000351611652, 0.00000135354389, 0.00000061357832, 0.00000220091101, 0.00000423862912, -0.00003000058678, 0.00000064187111, -0.00024228452821, 0.00003116775770, 0.00000176210710, 0.00000289141326, -0.00000150300525, 0.00000772853855, 0.00000098855242, 0.00000316606793, 0.00000002168364, 0.00000212047939, 0.00000299016097, 0.00000443224508, 0.00000028202930, 0.00000000000000, -0.00000249222421, -0.00000203136278, 0.00000131104809, 0.00000011987446, -0.00000370760154, 0.00004553918916, -0.00007711342914, -0.00004685295062, 0.00011049838213, -0.00000197486270, 0.00000395827146, 0.00000615046474, 0.00000755337123, 0.00000700606006, 0.00000922725030, -0.00000043310337, 0.00000107416383, 0.00000449787694, 0.00000305137178, 0.00001226376662, 0.00000000000000, 0.00000270820692, 0.00000208059305, 0.00000521478523, 0.00001779037302, 0.00000846544117, 0.00001120913385, -0.00065816845745, -0.00085107452469, -0.00013171190221, -0.00005540943675, -0.00001835885450, 0.00000101879823, 0.00000209222071, 0.00000091532502, -0.00000521515358, -0.00000209227142, -0.00000678545939, -0.00000034963549, -0.00000015111488, 0.00001560274177, 0.00000000000000, 0.00000350691471, -0.00001160475040, -0.00001763036562, 0.00003487367940, -0.00002787247831, -0.00000910982726, 0.00008818832430, -0.00524408789352, 0.00009378376126, 0.00004184526188, 0.00002849263365, -0.00002757280527, 0.00003388467667, 0.00000706207265, 0.00000625263419, -0.00003315929280, -0.00001181772132, 0.00000311426015, 0.00001875682574, -0.00000398287420, 0.00000000000000, -0.00001524541040, 0.00001724056165, 0.00002245173346, 0.00002806861812, -0.00000388776925, 0.00008143573359, -0.00005900909309, 0.00110496615525, 0.00134626252111, 0.00005128383054, -0.00001372421866, 0.00003612563887, 0.00002236580076, -0.00002728391883, 0.00001981237256, 0.00000655450458, 0.00000985319002, 0.00001347597299, 0.00000645987802, 0.00003304968050, 0.00000000000000, -0.00000530822217, 0.00001324870937, -0.00003610889689, -0.00005478735329, -0.00005818806312, -0.00037112057908, -0.00017812002625, -0.00093204283621, 0.00115969858598, -0.00033559172880, -0.00010441876657, -0.00001617923044, -0.00000555065844, 0.00007343527250, -0.00004408047607, 0.00000403802142, 0.00001843931204, 0.00001694047933, 0.00001213414362, -0.00000751115658, 0.00000000000000, 0.00005457974839, -0.00000334614515, 0.00005845565465, 0.00015000770509, 0.00021849104087, 0.00002724147635, 0.00167233624961, 0.00011666602222, 0.00276563479565, -0.00085952825611, -0.00030217235326, -0.00008841593808, 0.00000997664119, -0.00015285826521, 0.00002517224675, 0.00003009161810, 0.00001883217556, 0.00002146127554, 0.00001822445302, -0.00004128706860, 0.00000000000000, -0.00003496417776, 0.00001088761655, -0.00000298955979, -0.00005359326315, -0.00019021633489, -0.00017992728681, -0.00347794801928, 0.00064632791327, 0.00449698418379, -0.00017710507382, 0.00006126180233, 0.00018059254216, 0.00002354096432, 0.00008189838991, -0.00010060678323, -0.00017183290038, 0.00019413756672, 0.00021334811754, 0.00011263617489, 0.00000853522670, -0.00000000000000, -0.00006544789358, 0.00005424076880, -0.00000679056529, -0.00001249735487, -0.00053082982777, 0.00035396864405, -0.00115020677913, 0.05894451215863, 0.06573092192411, 0.01498018857092, 0.00278125284240, 0.00145188067108, 0.00033717858605, 0.00000800427370, -0.00009335305367, 0.00024286781263, -0.00023916347709, 0.00031213948387, 0.00018134393031, -0.00002521496390, -0.00000000000000, -0.00054337945767, 0.00012690725271, 0.00053313979879, 0.00064233405283, -0.00047686311882, 0.00176536326762, 0.00074157933705, -0.02684566564858, 1.00000000000000, 0.07176169008017, 0.00837037432939, -0.00000381640211, 0.00088998704450, -0.00049218931235, -0.00024546548957, -0.00036608282244, 0.00049480766756, 0.00031158892671, 0.00006898906577, 0.00021280418150, 0.00028127161204, -0.00070030166535, 0.00022237010126, -0.00028713891516, -0.00013800295710, 0.00005912094275, 0.00172126013786, -0.00618684850633, 0.03608432412148, Aspect_ratio , 0.49896776676178, 0.00091372377938, -0.00085712829605, -0.00124801427592, -0.00007427225501, -0.00005245858847, 0.00002841771493, 0.00020249813679, -0.00014303345233, 0.00001406490901, 0.00023699452868, 0.00008661757602, 0.00025744654704, -0.00022715188970, -0.00076146807987, 0.00055185536621, -0.00012325309217, -0.00072356045712, -0.00160693109501, 0.00246682553552, -0.14175094664097, -0.36207047104836, -0.04089594259858, 0.00060774467420, 0.00088646943914, 0.00004865296432, -0.00041878610500, -0.00023025234987, -0.00009676301852, -0.00000000000000, 0.00008409228758, 0.00011432896281, -0.00000707848403, 0.00004698805787, -0.00043642931269, 0.00081384339137, -0.00065635429928, -0.00011831733718, 0.00017413357273, 0.00224463525228, 0.00478497287259, 0.03294761106372, 0.01078986655921, 0.10731782764196, 0.00075034319889, -0.00009241879889, 0.00055023463210, 0.00006596000458, 0.00005045382932, 0.00014874986664, 0.00000000000000, -0.00015369028552, 0.00001037383754, 0.00009250180301, 0.00026204055757, 0.00007424291834, -0.00047751804232, 0.00029184055165, 0.00050921301590, -0.00004825839278, -0.00029933769838, 0.00279659987427, 0.00210463814437, -0.00618590926751, -0.02400829829276, -0.02316811867058, -0.00086368201301, -0.00032258985448, -0.00018304496189, 0.00008438774967, -0.00008305341908, 0.00000000000000, 0.00013047417451, -0.00001376930322, -0.00001723831701, -0.00011543079017, -0.00022646733851, 0.00013467084500, -0.00004661652201, -0.00008419520600, 0.00035772417323, -0.00011815709877, 0.00028718306567, 0.00092207465786, -0.00317224999890, 0.00061770365573, 0.01017294172198, 0.00294739892706, 0.00014669894881, 0.00015702951350, 0.00003432080121, -0.00008555022214, -0.00000000000000, 0.00000454909878, -0.00000196001542, -0.00003198397462, -0.00004425687075, -0.00004129848094, -0.00003789070615, -0.00027583551127, 0.00025874207495, -0.00002334945384, -0.00007259396807, -0.00008295358566, 0.00011360697681, -0.00101968157105, 0.00046784928418, -0.00208410434425, -0.00313158822246, -0.00046005158219, -0.00010552268213, -0.00005850767775, 0.00003971093611, 0.00000000000000, -0.00005275657168, -0.00001065901233, -0.00001934838656, -0.00001220186732, -0.00002060524639, -0.00000225423423, -0.00001894621164, -0.00001533334580, -0.00001791087379, 0.00008156246622, -0.00008441298269, 0.00021060956351, -0.00030303673702, 0.00075949780876, -0.00010539998038, 0.00109045265708, 0.00068949378328, 0.00009268362192, 0.00003471063246, 0.00001204656473, -0.00000000000000, 0.00001500743110, 0.00000105878155, -0.00000910870767, -0.00000172467264, -0.00000722095228, 0.00000699280463, -0.00002061720625, -0.00000889817693, -0.00001993474507, 0.00000370749740, -0.00000090311920, 0.00002677819793, 0.00043428712524, 0.00210293265991, 0.00018200518389, -0.00009621794743, -0.00035250501242, -0.00012996385340, -0.00002185157609, -0.00001116586463, -0.00000000000000, -0.00000451994811, 0.00000424055270, -0.00000463139304, 0.00000301006116, -0.00000123974939, 0.00000632465435, -0.00002090823000, 0.00001773388794, 0.00000121050368, 0.00001886057362, -0.00001043497195, -0.00002269273500, -0.00021979617304, -0.00001043962493, -0.00116343051195, -0.00004193381756, 0.00007944958634, 0.00007301353617, 0.00002082651736, -0.00000119863023, -0.00000000000000, -0.00001440504820, -0.00000391270805, -0.00000490489265, -0.00000504441778, -0.00000904507579, -0.00000111389932, 0.00000597532107, 0.00000047090245, -0.00001553130096, -0.00001524566323, -0.00000522222899, -0.00007707672921, -0.00004165665086, 0.00015764687851, 0.00035649110214, 0.00038701237645, 0.00002386798405, -0.00001946414341, -0.00000913835174, -0.00000489907188, 0.00000000000000, 0.00000172327657, -0.00000015388650, -0.00000603232729, -0.00000397650865, 0.00000280493782, 0.00000463132073, -0.00000788678426, -0.00000471605335, -0.00000283715985, -0.00000422824724, 0.00000366817630, -0.00001159603562, -0.00001625759251, 0.00049116823357, 0.00005048640014, -0.00020234247495, -0.00006341376866, -0.00000807822744, 0.00000070463199, 0.00000014041755, 0.00000000000000, -0.00000718306910}; #else Real Aspect_ratio = 5; Real coeffmat[21][21] = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, s, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, Aspect_ratio, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2*s, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 , 0, 0, 0, 0, 0, 0, 0, 0, 0}; #endif Z = 0; Real R = 0; for (long i = -10; i <= 10; i++) { for (long j = -10; j <= 10; j++) { R += coeffmat[i+10][j+10] * sctl::cos(-i*phi + Nperiod*j*theta); Z += coeffmat[i+10][j+10] * sctl::sin(-i*phi + Nperiod*j*theta); } } X = R * sctl::cos(theta); Y = R * sctl::sin(theta); } GenericKernel BiotSavart ; GenericKernel BiotSavartGrad; GenericKernel Laplace_FxU ; GenericKernel Laplace_FxdU; GenericKernel Laplace_dUxF; GenericKernel Laplace_dUxD; GenericKernel Laplace_Fxd2U; mutable Quadrature quadrature_BS ; mutable Quadrature quadrature_dBS ; mutable Quadrature quadrature_FxU ; mutable Quadrature quadrature_FxdU; mutable Quadrature quadrature_dUxF; mutable Quadrature quadrature_dUxD; mutable Quadrature quadrature_Fxd2U; mutable Vector Bt0, Bp0, dBt0, dBp0; mutable Vector sigma, B, gvec, dgdB; mutable Real alpha, beta; ElemLst elements; Vector NtNp_; Vector elem_dsp; }; template class MHDEquilib { static constexpr Integer fourier_upsample = 2; static constexpr Integer COORD_DIM = 3; static constexpr Integer ELEM_DIM = COORD_DIM-1; using ElemBasis = Basis; static Vector cheb2grid(const Vector& X, Long Mt, Long Mp, Long Nt, Long Np) { const Long dof = X.Dim() / (Mt * Mp); SCTL_ASSERT(X.Dim() == Mt * Mp *dof); Vector Xf(dof*Nt*Np); Xf = 0; const Long Nnodes = ElemBasis::Size(); const Matrix& Mnodes = Basis::Nodes(); for (Long t = 0; t < Nt; t++) { for (Long p = 0; p < Np; p++) { Real theta = t / (Real)Nt; Real phi = p / (Real)Np; Long i = (Long)(theta * Mt); Long j = (Long)(phi * Mp); Real x = theta * Mt - i; Real y = phi * Mp - j; Long elem_idx = i * Mp + j; Vector Interp0(ORDER); Vector Interp1(ORDER); { // Set Interp0, Interp1 auto node = [&Mnodes] (Long i) { return Mnodes[0][i]; }; for (Long i = 0; i < ORDER; i++) { Real wt_x = 1, wt_y = 1; for (Long j = 0; j < ORDER; j++) { if (j != i) { wt_x *= (x - node(j)) / (node(i) - node(j)); wt_y *= (y - node(j)) / (node(i) - node(j)); } Interp0[i] = wt_x; Interp1[i] = wt_y; } } } for (Long ii = 0; ii < ORDER; ii++) { for (Long jj = 0; jj < ORDER; jj++) { Long node_idx = jj * ORDER + ii; for (Long k = 0; k < dof; k++) { Xf[(k*Nt+t)*Np+p] += X[elem_idx*dof+k][node_idx] * Interp0[ii] * Interp1[jj]; } } } } } return Xf; } static Vector grid2cheb(const Vector& Xf, Long Nt, Long Np, Long Mt, Long Mp) { Long dof = Xf.Dim() / (Nt*Np); SCTL_ASSERT(Xf.Dim() == dof*Nt*Np); Vector X(Mt*Mp*dof); constexpr Integer INTERP_ORDER = 12; for (Long tt = 0; tt < Mt; tt++) { for (Long pp = 0; pp < Mp; pp++) { for (Long t = 0; t < ORDER; t++) { for (Long p = 0; p < ORDER; p++) { Matrix Mnodes = Basis::Nodes(); Real theta = (tt + Mnodes[0][t]) / Mt; Real phi = (pp + Mnodes[0][p]) / Mp; Long i = (Long)(theta * Nt); Long j = (Long)(phi * Np); Real x = theta * Nt - i; Real y = phi * Np - j; Vector Interp0(INTERP_ORDER); Vector Interp1(INTERP_ORDER); { // Set Interp0, Interp1 auto node = [] (Long i) { return (Real)i - (INTERP_ORDER-1)/2; }; for (Long i = 0; i < INTERP_ORDER; i++) { Real wt_x = 1, wt_y = 1; for (Long j = 0; j < INTERP_ORDER; j++) { if (j != i) { wt_x *= (x - node(j)) / (node(i) - node(j)); wt_y *= (y - node(j)) / (node(i) - node(j)); } Interp0[i] = wt_x; Interp1[i] = wt_y; } } } for (Long k = 0; k < dof; k++) { Real X0 = 0; for (Long ii = 0; ii < INTERP_ORDER; ii++) { for (Long jj = 0; jj < INTERP_ORDER; jj++) { Long idx_i = (i + ii-(INTERP_ORDER-1)/2 + Nt) % Nt; Long idx_j = (j + jj-(INTERP_ORDER-1)/2 + Np) % Np; X0 += Interp0[ii] * Interp1[jj] * Xf[(k*Nt+idx_i)*Np+idx_j]; } } Long elem_idx = tt * Mp + pp; Long node_idx = p * ORDER + t; X[elem_idx*dof+k][node_idx] = X0; } } } } } return X; } static void fourier_filter(sctl::Vector& X, long Nt_, long Np_, Real sigma, const Comm& comm) { long dof = X.Dim() / (Nt_ * Np_); SCTL_ASSERT(X.Dim() == dof * Nt_ * Np_); sctl::FFT fft_r2c, fft_c2r; sctl::StaticArray fft_dim = {Nt_, Np_}; fft_r2c.Setup(sctl::FFT_Type::R2C, 1, sctl::Vector(2, fft_dim, false), omp_get_max_threads()); fft_c2r.Setup(sctl::FFT_Type::C2R, 1, sctl::Vector(2, fft_dim, false), omp_get_max_threads()); long Nt = Nt_; long Np = fft_r2c.Dim(1) / (Nt * 2); SCTL_ASSERT(fft_r2c.Dim(1) == Nt * Np * 2); auto filter_fn = [](Real x2, Real sigma) {return exp(-x2/(2*sigma*sigma));}; //auto filter_fn = [](Real x2, Real sigma) {return (x2 normal, gradX; biest::SurfaceOp op(comm, Nt_, Np_); sctl::Vector coeff(fft_r2c.Dim(1)); for (long k = 0; k < dof; k++) { sctl::Vector X_(Nt_*Np_, X.begin() + k*Nt_*Np_, false); fft_r2c.Execute(X_, coeff); for (long t = 0; t < Nt; t++) { for (long p = 0; p < Np; p++) { Real tt = (t - (t > Nt / 2 ? Nt : 0)) / (Real)(Nt / 2); Real pp = p / (Real)Np; Real f = filter_fn(tt*tt+pp*pp, sigma); coeff[(t * Np + p) * 2 + 0] *= f; coeff[(t * Np + p) * 2 + 1] *= f; } } fft_c2r.Execute(coeff, X_); } }; static void filter(const Stellarator& S, const Comm& comm, Vector& f, Real sigma) { Long dof = f.Dim() / S.NElem(); SCTL_ASSERT(f.Dim() == S.NElem() * dof); for (Long i = 0; i < S.Nsurf()-1; i++) { const Long Mt = S.NTor(i); const Long Mp = S.NPol(i); const Long Nelem = Mt * Mp; const Long offset = S.ElemDsp(i); const Long Nt = Mt * ORDER * 4; const Long Np = Mp * ORDER * 4; Vector f_(Nelem*dof, f.begin() + offset*dof, false); Vector f_fourier = cheb2grid(f_, Mt, Mp, Nt, Np); fourier_filter(f_fourier, Nt, Np, 0.25 * sigma, comm); f_ = grid2cheb(f_fourier, Nt, Np, Mt, Mp); } }; static Vector cheb2grid(const Stellarator& S, const Vector& f) { const Long Nnodes = ElemBasis::Size(); const Long Nelem = S.NElem(); const Long dof = f.Dim() / Nelem; SCTL_ASSERT(Nelem * dof == f.Dim()); Vector f_fourier(dof * Nelem * (ORDER*ORDER*fourier_upsample*fourier_upsample)); for (Long i = 0; i < S.Nsurf(); i++) { const Long Mt = S.NTor(i); const Long Mp = S.NPol(i); const Long offset = S.ElemDsp(i); const Long Nt = Mt * ORDER * fourier_upsample; const Long Np = Mp * ORDER * fourier_upsample; const Vector f_(Mt*Mp*dof, (Iterator)f.begin() + offset*dof, false); Vector f_fourier_(dof*Nt*Np, f_fourier.begin() + dof*offset * (ORDER*ORDER*fourier_upsample*fourier_upsample), false); f_fourier_ = cheb2grid(f_, Mt, Mp, Nt, Np); SCTL_ASSERT(f_fourier_.Dim() == dof*Nt*Np); } return f_fourier; } static Vector grid2cheb(const Stellarator& S, const Vector& f_fourier) { const Long Nnodes = ElemBasis::Size(); const Long Nelem = S.NElem(); const Long dof = f_fourier.Dim() / (Nelem * (ORDER*ORDER*fourier_upsample*fourier_upsample)); SCTL_ASSERT(dof * Nelem * (ORDER*ORDER*fourier_upsample*fourier_upsample) == f_fourier.Dim()); Vector f(Nelem * dof); for (Long i = 0; i < S.Nsurf(); i++) { const Long Mt = S.NTor(i); const Long Mp = S.NPol(i); const Long offset = S.ElemDsp(i); const Long Nt = Mt * ORDER * fourier_upsample; const Long Np = Mp * ORDER * fourier_upsample; Vector f_(Mt*Mp*dof, f.begin() + offset*dof, false); const Vector f_fourier_(dof*Nt*Np, (Iterator)f_fourier.begin() + dof*offset * (ORDER*ORDER*fourier_upsample*fourier_upsample), false); f_ = grid2cheb(f_fourier_, Nt, Np, Mt, Mp); SCTL_ASSERT(f_.Dim() == Mt*Mp*dof); } return f; } template static Long GradientDescent(GradOp& grad_op, Eigen::VectorXd& x, Real& fx, Long max_iter, Real tol) { Real dt = 0.1; for (Long iter = 0; iter < max_iter; iter++) { Eigen::VectorXd grad(x.size()); fx = grad_op(x, grad); { // Update dt Eigen::VectorXd grad_(x.size()); Eigen::VectorXd x1 = x - grad * dt * 0.5; Eigen::VectorXd x2 = x - grad * dt * 1.0; Real fx1 = grad_op(x1, grad_); Real fx2 = grad_op(x2, grad_); { // Calculate optimal step size dt Real a = 2*fx - 4*fx1 + 2*fx2; Real b =-3*fx + 4*fx1 - fx2; Real c = fx; Real s = -b/(2*a); dt *= s; Real fx_ = a*s*s + b*s + c; std::cout<<"g = "<& S, const Vector& pressure, const Vector& flux_tor, const Vector& flux_pol) { S_ = S; pressure_ = pressure; flux_tor_ = flux_tor; flux_pol_ = flux_pol; iter = 0; } Real operator()(const Eigen::VectorXd& x, Eigen::VectorXd& grad) { const Comm comm = Comm::World(); const Long Nelem = S_.NElem(); const Long Nnodes = ElemBasis::Size(); const Long N = Nelem * COORD_DIM * Nnodes; Vector X_fourier(x.size()); for (Long i = 0; i < x.size(); i++) { // Set X_fourier X_fourier[i] = x(i); } Vector X = grid2cheb(S_, X_fourier); Real g; for (Long i = 0; i < Nelem; i++) { // Set S_ for (Long j = 0; j < Nnodes; j++) { S_.Elem(i,0)[j] = X[i*COORD_DIM+0][j]; S_.Elem(i,1)[j] = X[i*COORD_DIM+1][j]; S_.Elem(i,2)[j] = X[i*COORD_DIM+2][j]; } } Vector dgdnu = Stellarator::compute_gradient(S_, pressure_, flux_tor_, flux_pol_, &g); //Vector dgdnu = Stellarator::compute_pressure_jump(S_, pressure_, flux_tor_, flux_pol_, &g); Vector dXdt_fourier; { // Set grad //filter(S_, comm, dgdnu, 0.1); //Vector dgdnu_fourier = cheb2grid(S_, dgdnu); { // deprecate Vector dXdt(Nelem*COORD_DIM); { // Set dXdt dXdt = 0; const Long Nnodes = ElemBasis::Size(); Vector normal, area_elem; Stellarator::compute_norm_area_elem(S_, normal, area_elem); for (Long i = 0; i < S_.ElemDsp(S_.Nsurf()-1); i++) { for (Long j = 0; j < Nnodes; j++) { dXdt[i*COORD_DIM+0][j] = normal[i*COORD_DIM+0][j] * dgdnu[i][j]; dXdt[i*COORD_DIM+1][j] = normal[i*COORD_DIM+1][j] * dgdnu[i][j]; dXdt[i*COORD_DIM+2][j] = normal[i*COORD_DIM+2][j] * dgdnu[i][j]; } } } { // Filter dXdt filter(S_, comm, dXdt, 0.1); } dXdt_fourier = cheb2grid(S_, dXdt); SCTL_ASSERT(grad.size() == dXdt_fourier.Dim()); for (Long i = 0; i < grad.size(); i++) { // Set grad grad(i) = dXdt_fourier[i]; } } } if (1) { // Write VTU VTUData vtu; vtu.AddElems(S_.GetElemList(), dgdnu, ORDER); vtu.WriteVTK("dgdnu"+std::to_string(iter), comm); } if (1) { // Write VTU VTUData vtu; Vector dXdt = grid2cheb(S_, dXdt_fourier); vtu.AddElems(S_.GetElemList(), dXdt, ORDER); vtu.WriteVTK("dXdt"+std::to_string(iter), comm); } std::cout<<"iter = "< rhs; BIOpL(rhs, U); Vector sigma; { // Set sigma Long Nnode = DensityBasis::Size(); Long Nelem = S.GetElem().NElem(); typename sctl::ParallelSolver::ParallelOp A = [&S,&BIOp_half_S_D](sctl::Vector* Ax, const sctl::Vector& x) { Long Nnode = DensityBasis::Size(); Long Nelem = S.GetElem().NElem(); Ax->ReInit(Nelem*COORD_DIM*Nnode); Vector x_(Nelem*COORD_DIM), Ax_(Nelem*COORD_DIM); for (Long i = 0; i < Nelem*COORD_DIM; i++) { // Set x_ for (Long k = 0; k < Nnode; k++) { x_[i][k] = x[i*Nnode+k]; } } BIOp_half_S_D(Ax_, x_); for (Long i = 0; i < Nelem*COORD_DIM; i++) { // Set Ax for (Long k = 0; k < Nnode; k++) { (*Ax)[i*Nnode+k] = Ax_[i][k]; } } }; Vector sigma_(Nelem*COORD_DIM*Nnode), rhs_(Nelem*COORD_DIM*Nnode); for (Long i = 0; i < Nelem*COORD_DIM; i++) {// Set rhs_ for (Long k = 0; k < Nnode; k++) { rhs_[i*Nnode+k] = rhs[i][k]; } } sigma_ = 0; ParallelSolver linear_solver(comm, true); linear_solver(&sigma_, A, rhs_, 1e-6, 50); sigma.ReInit(Nelem * COORD_DIM); for (Long i = 0; i < Nelem*COORD_DIM; i++) {// Set sigma for (Long k = 0; k < Nnode; k++) { sigma[i][k] = sigma_[i*Nnode+k]; } } } Vector U1; BIOp_half_S_D(U1, sigma); { // Write VTU VTUData vtu_sigma; vtu_sigma.AddElems(S.elements, sigma, ORDER); vtu_sigma.WriteVTK("sphere-sigma1", comm); VTUData vtu_U; vtu_U.AddElems(S.elements, U1, ORDER); vtu_U.WriteVTK("sphere-U1", comm); } } Profile::print(&comm); } private: template void SurfInteg(Vector& I, const Vector& f) { static_assert(std::is_same::value, "FnBasis is different from CoordBasis"); const Long Nelem = elements.NElem(); const Long dof = f.Dim() / Nelem; SCTL_ASSERT(f.Dim() == Nelem * dof); auto nodes = FnBasis::Nodes(); auto quad_wts = FnBasis::QuadWts(); const Long Nnodes = FnBasis::Size(); auto EvalOp = CoordBasis::SetupEval(nodes); Vector dX; const auto& X = elements.ElemVector(); SCTL_ASSERT(X.Dim() == Nelem * COORD_DIM); CoordBasis::Grad(dX, X); Matrix I_(Nelem, dof); for (Long i = 0; i < Nelem; i++) { for (Long k = 0; k < dof; k++) { I_[i][k] = 0; } for (Long j = 0; j < Nnodes; j++) { Real dA = 0; StaticArray Xn; Xn[0] = dX[i*COORD_DIM*2+2][j] * dX[i*COORD_DIM*2+5][j] - dX[i*COORD_DIM*2+3][j] * dX[i*COORD_DIM*2+4][j]; Xn[1] = dX[i*COORD_DIM*2+4][j] * dX[i*COORD_DIM*2+1][j] - dX[i*COORD_DIM*2+5][j] * dX[i*COORD_DIM*2+0][j]; Xn[2] = dX[i*COORD_DIM*2+0][j] * dX[i*COORD_DIM*2+3][j] - dX[i*COORD_DIM*2+1][j] * dX[i*COORD_DIM*2+2][j]; dA += sqrt(Xn[0]*Xn[0] + Xn[1]*Xn[1] + Xn[2]*Xn[2]) * quad_wts[j]; for (Long k = 0; k < dof; k++) { I_[i][k] += dA * f[i*dof+k][j]; } } } Long Ns = elem_cnt.Dim(); if (I.Dim() != Ns * dof) I.ReInit(Ns * dof); I = 0; Long elem_itr = 0; for (Long i = 0; i < Ns; i++) { for (Long j = 0; j < elem_cnt[i]; j++) { for (Long k = 0; k < dof; k++) { I[i*dof+k] += I_[elem_itr][k]; } elem_itr++; } } } void InitSpheres(const Vector X, const Vector& R){ SCTL_ASSERT(X.Dim() == R.Dim() * COORD_DIM); Long N = R.Dim(); elements.ReInit(2*COORD_DIM*N); auto nodes = ElemLst::CoordBasis::Nodes(); for (Long l = 0; l < N; l++) { for (Integer i = 0; i < COORD_DIM; i++) { for (Integer j = 0; j < 2; j++) { for (int k = 0; k < ElemLst::CoordBasis::Size(); k++) { Real coord[COORD_DIM]; coord[(i+0)%COORD_DIM] = (j ? -1.0 : 1.0); coord[(i+1)%COORD_DIM] = 2.0 * nodes[j?1:0][k] - 1.0; coord[(i+2)%COORD_DIM] = 2.0 * nodes[j?0:1][k] - 1.0; Real R0 = sqrt(coord[0]*coord[0] + coord[1]*coord[1] + coord[2]*coord[2]); elements((l*COORD_DIM+i)*2+j,0)[k] = X[l*COORD_DIM+0] + R[l] * coord[0] / R0; elements((l*COORD_DIM+i)*2+j,1)[k] = X[l*COORD_DIM+1] + R[l] * coord[1] / R0; elements((l*COORD_DIM+i)*2+j,2)[k] = X[l*COORD_DIM+2] + R[l] * coord[2] / R0; } } } } elem_cnt.ReInit(N); elem_cnt = 6; } GenericKernel Stokes_DxU; GenericKernel Stokes_FxU; GenericKernel Stokes_FxT; Quadrature quadrature_DxU; Quadrature quadrature_FxU; Quadrature quadrature_FxT; ElemLst elements; Vector elem_cnt; }; } // end namespace #endif //_SCTL_BOUNDARY_QUADRATURE_HPP_