// // Created by blaz on 5/20/20. // // // Benchmark case for FDFD in 2D, scattering from a dielectric cylinder for a TE planewave. // The analytical solution is reasonably easy to find in terms of Hankel functions. // The solution uses the radial PML variant #include #include #include #include #include #include #include using namespace mm; // NOLINT constexpr std::complex operator"" _i(long double d) { return std::complex{0.0, static_cast(d)}; } int main() { XML conf("dielectric_cylinder.xml"); std::string outfile = conf.get("IO.outfile"); // solver int ffact = conf.get("solver.ffact"); int maxiter = conf.get("solver.maxiter"); double dtol = conf.get("solver.droptol"); double soltol = conf.get("solver.errtol"); // Physics double mu0 = 1.25663706212*1e-6; // [H/m = N/A^2] double eps0 = 8.8541878128*1e-12; // [F/m] double c0 = 1.0/sqrt(mu0*eps0); double epsr = conf.get("case.epsr"); double wavelength = conf.get("fourier.wavelength"); double k0 = 2*PI/wavelength; double omega = k0*c0; // Geometry double r = conf.get("geometry.r"); double a = conf.get("geometry.a"); double tol = 1e-5; double t0 = conf.get("case.theta0"); std::cout << "Solving for:" << "\nk0 = " << k0 << "\nomega = " << omega << "\nwavelength = " << wavelength << std::endl; // Approximation int supp = conf.get("approx.supp"); int num_monoms = conf.get("approx.num_monoms"); // Fill double dx_cyl = conf.get("case.dx_cylinder"); double dx_omega = conf.get("case.dx_omega"); double alpha = conf.get("case.alpha"); double power = conf.get("case.power"); auto dx = [=](const Vec2d &p){return (dx_cyl-dx_omega)/ (1.0+std::pow(std::sqrt(p[0]*p[0]+p[1]*p[1])/alpha, power))+dx_omega;}; // PML double d = conf.get("PML.d"); double m = conf.get("PML.m"); double sigma_max = conf.get("PML.sigma"); BallShape outbox(0, r); // outer BallShape cylinder(0, a); // inner auto outer_shape = outbox - cylinder; // fill outer boundary DomainDiscretization cyl_dom = cylinder.discretizeBoundaryWithDensity(dx); DomainDiscretization outer_dom = (outbox-cylinder).discretizeBoundaryWithDensity(dx); // Node bookkeeping // delete nodes of outer domain on shared boundary Range to_del = outer_dom.positions().filter([=] (const Vec2d& p){ return ((p).norm()-a < tol);}); outer_dom.removeNodes(to_del); // make index arrays Range interface_c_idx(0); // cylinder index Range interface_o_idx(0); // outer index for (int i : cyl_dom.boundary()) { int id = outer_dom.addBoundaryNode(cyl_dom.pos(i), -6, -1 * cyl_dom.normal(i)); interface_c_idx.push_back(i); interface_o_idx.push_back(id); } // outer boundary Range outer_bnd = outer_dom.positions().filter([=] (const Vec2d& p) {return r-p.norm() < tol;}); // fill domains GeneralFill fill; fill(cyl_dom, dx); fill(outer_dom, dx); int Ncyl = cyl_dom.size(); int Nouter = outer_dom.size(); std::cout << "Domain size = " << Ncyl + Nouter << "\n"; BasicRelax relax; relax(cyl_dom, dx); relax(outer_dom, dx); // filter for PML auto PML = (outer_dom.positions()).filter([=] (const Vec2d& p){ return (std::sqrt(r - d) < std::sqrt(p.norm()) && std::sqrt(p.norm()) < std::sqrt(r-tol));}); auto outer_inter = (outer_dom.positions()).filter([=] (const Vec2d& p) { return (std::sqrt(a + tol) < std::sqrt(p.norm()) && std::sqrt(p.norm()) <= std::sqrt(r-d));}); // find support cyl_dom.findSupport(FindClosest(supp)); outer_dom.findSupport(FindClosest(supp)); // make rbffd engine Polyharmonic phs; RBFFD, Vec2d, ScaleToClosest, Eigen::PartialPivLU> rbffd(phs, num_monoms); // make M and rhs Eigen::SparseMatrix , Eigen::RowMajor> M(Ncyl + Nouter, Ncyl + Nouter); Eigen::VectorXcd rhs(Ncyl + Nouter); rhs.setZero(); // make operators auto storage_inner = cyl_dom.computeShapes(rbffd); auto storage_outer = outer_dom.computeShapes(rbffd); auto op_inner = storage_inner.implicitOperators(M, rhs); auto op_outer = storage_outer.implicitOperators(M, rhs); // offset operator write op_inner.setRowOffset(Nouter); op_inner.setColOffset(Nouter); // reserve storage auto store = storage_inner.supportSizes() + storage_outer.supportSizes(); Range reserve = store; for (int& i : reserve) i *= 2; M.reserve(reserve); ExplicitOperators exop(storage_outer); Eigen::VectorXcd sw(Nouter); for (int i : outer_inter) { sw[i] = 1.0; } for (int i : PML) { Vec2d pos = outer_dom.pos(i); double rad = std::sqrt(std::pow(pos[0], 2) + std::pow(pos[1], 2)); sw[i] = 1.0-1.0_i*sigma_max/omega/eps0*std::pow((rad+d-r)/d, m); } // Construct system double l_in = epsr*std::pow(omega/c0, 2.0); double l_out = std::pow(omega/c0, 2.0); for (int i : cyl_dom.interior()) { op_inner.lap(i) + l_in*op_inner.value(i) = 0.0; } for (int i : outer_inter) { op_outer.lap(i) + l_out*op_outer.value(i) = 0.0; } for (int c = 0; c < interface_c_idx.size(); ++c) { int i = interface_c_idx[c]; // cyl index int j = interface_o_idx[c]; // outer index Vec2d pos = cyl_dom.pos(i); double x = pos[0]; double y = pos[1]; // calculate normal derivative of the source funciton at the i-th node // get normal, must point out of the cylinder Vec2d normal = cyl_dom.normal(i); std::complex incident = std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0))); std::complex din_dx = 1.0_i*k0*std::cos(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0))); std::complex din_dy = 1.0_i*k0*std::sin(t0)*std::exp(1.0_i*k0*(x*std::cos(t0)+y*std::sin(t0))); std::complex dui_dn = normal[0]*din_dx + normal[1]*din_dy; // continuity of the fields, and the incident field op_inner.value(i) + (-1)*op_outer.value(j, Nouter + i) = incident; // continuity of derivatives op_outer.neumann(j, outer_dom.normal(j)) + (1/epsr)*op_inner.neumann(i, cyl_dom.normal(i), j - Nouter) = dui_dn; } // SC-PML - where either x or y is constant for (int i : PML) { 1.0/(sw[i]*sw[i])*op_outer.lap(i) + l_out*op_outer.value(i) + ((-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 0, i)*op_outer.der1(i, 0) +(-1.0)/(sw[i]*sw[i]*sw[i])*exop.d1(sw, 1, i)*op_outer.der1(i, 1)) = 0.0; } for (int i : outer_bnd) { op_outer.value(i) = 0.0; } HDF out(outfile+".h5", HDF::DESTROY); out.writeXML("config", conf); out.writeDomain("cyl_dom", cyl_dom); out.writeDomain("outer_dom", outer_dom); Eigen::BiCGSTAB>> solver; // customize solver solver.preconditioner().setFillfactor(ffact); solver.preconditioner().setDroptol(dtol); solver.setMaxIterations(maxiter); solver.setTolerance(soltol); solver.compute(M); Eigen::VectorXcd sol = solver.solve(rhs); std::cout << "Iterations: " << solver.iterations() << "\n"; std::cout << "Error: " << solver.error() << "\n"; ScalarField rsol = sol.real(); ScalarField csol = sol.imag(); // out types for debugging out.writeIntArray("pml_nodes", PML); out.writeIntArray("outer_inter_nodes", outer_inter); out.writeDoubleArray("rsol", rsol); out.writeDoubleArray("csol", csol); out.writeDoubleAttribute("wavelength", wavelength); out.close(); return 0; }