From 56a05f47287745895d65ef8e454a3ef2bc12c2fb Mon Sep 17 00:00:00 2001 From: psocratis Date: Sun, 17 May 2020 18:31:16 -0700 Subject: [PATCH] Added new versions of DiagST as indicated in the paper 'Additive Overlapping DDM' (i.e, the source is calculated at the subdomain that receives it and it's multiplied by the indicator function of the non-overlapping partition) --- examples/maxwell-solver/DST/DSTnew.cpp | 467 ++++++ examples/maxwell-solver/DST/DSTnew.hpp | 52 + examples/maxwell-solver/DST/DiagST.cpp | 4 +- examples/maxwell-solver/DST/DiagSTnew.cpp | 1733 +++++++++++++++++++++ examples/maxwell-solver/DST/DiagSTnew.hpp | 120 ++ examples/maxwell-solver/DST/Utilities.cpp | 4 +- examples/maxwell-solver/DST/helmholtz.cpp | 10 +- 7 files changed, 2381 insertions(+), 9 deletions(-) create mode 100644 examples/maxwell-solver/DST/DSTnew.cpp create mode 100644 examples/maxwell-solver/DST/DSTnew.hpp create mode 100644 examples/maxwell-solver/DST/DiagSTnew.cpp create mode 100644 examples/maxwell-solver/DST/DiagSTnew.hpp diff --git a/examples/maxwell-solver/DST/DSTnew.cpp b/examples/maxwell-solver/DST/DSTnew.cpp new file mode 100644 index 0000000000..15d9a3b60c --- /dev/null +++ b/examples/maxwell-solver/DST/DSTnew.cpp @@ -0,0 +1,467 @@ +//Diagonal Source Transfer Preconditioner + +#include "DST.hpp" + + +DST::DST(SesquilinearForm * bf_, Array2D & Pmllength_, + double omega_, Coefficient * ws_, int nrlayers_) + : Solver(2*bf_->FESpace()->GetTrueVSize(), 2*bf_->FESpace()->GetTrueVSize()), + bf(bf_), Pmllength(Pmllength_), omega(omega_), ws(ws_), nrlayers(nrlayers_) +{ + Mesh * mesh = bf->FESpace()->GetMesh(); + dim = mesh->Dimension(); + + // ----------------- Step 1 -------------------- + // Introduce 2 layered partitios of the domain + // + int partition_kind; + + // 1. Ovelapping partition with overlap = 2h + partition_kind = 2; // Non Overlapping partition + int nx=4; + int ny=1; + int nz=1; + + povlp = new MeshPartition(mesh, partition_kind,nx,ny,nz, nrlayers); + nxyz[0] = povlp->nxyz[0]; + nxyz[1] = povlp->nxyz[1]; + nxyz[2] = povlp->nxyz[2]; + nrpatch = povlp->nrpatch; + subdomains = povlp->subdomains; + + // + // ----------------- Step 1a ------------------- + // Save the partition for visualization + // SaveMeshPartition(povlp->patch_mesh, "output/mesh_ovlp.", "output/sol_ovlp."); + + ovlp_prob = new DofMap(bf,povlp); + PmlMat.SetSize(nrpatch); + PmlMatInv.SetSize(nrpatch); + for (int ip=0; ipSetOperator(*PmlMat[ip]); + } + nsweeps = pow(2,dim); + sweeps.SetSize(nsweeps,dim); + // 2D + sweeps(0,0) = 1; sweeps(0,1) = 1; + sweeps(1,0) = -1; sweeps(1,1) = 1; + sweeps(2,0) = 1; sweeps(2,1) =-1; + sweeps(3,0) = -1; sweeps(3,1) =-1; + + // Set up src arrays size + f_orig.SetSize(nrpatch); + f_transf.SetSize(nrpatch); + // Construct a simple map used for directions of transfer + for (int ip=0; ipfespaces[ip]->GetTrueVSize(); // (x 2 for complex ) + f_orig[ip] = new Vector(n); *f_orig[ip] = 0.0; + f_transf[ip].SetSize(nsweeps); + for (int i=0;i * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; + r.GetSubVector(*Dof2GlobalDof,*f_orig[ip]); + } + + char vishost[] = "localhost"; + int visport = 19916; + z = 0.0; + Vector znew(z); + Vector z1(z); + Vector z2(z); + + // -------------------------------------------- + // Sweep in the direction (1,1) + // -------------------------------------------- + int nx = nxyz[0]; + int ny = nxyz[1]; + + int nsteps = nx + ny - 1; + + for (int l=0; l<1; l++) + { + for (int s = 0; s=ny) continue; + // cout << "Patch no: (" << i <<"," << j << ")" << endl; + + // find patch id + Array ij(2); ij[0] = i; ij[1]=j; + int ip = GetPatchId(ij); + // cout << "ip = " << ip << endl; + + // Solve the PML problem in patch ip with all sources + // Original and all transfered (maybe some of them) + Array * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; + int ndofs = Dof2GlobalDof->Size(); + + Vector sol_local(ndofs); sol_local = 0.0; + Vector res_local(ndofs); res_local = 0.0; + if (l==0) res_local += *f_orig[ip]; + // res_local += *f_orig[ip]; + res_local += *f_transf[ip][l]; + // Extend by zero to the PML mesh + // if (res_local.Norml2() < 1e-11) continue; + PmlMatInv[ip]->Mult(res_local, sol_local); + + TransferSources(l,ip, sol_local); + + // cut off the ip solution to all possible directions + Arraydirections(2); directions = 0; + + if (i+10) directions[0] = -1; + if (j>0) directions[1] = -1; + GetCutOffSolution(sol_local,cfsol_local,ip,directions,nrlayers,true); + znew = 0.0; + // znew.SetSubVector(*Dof2GlobalDof, cfsol_local); + znew.SetSubVector(*Dof2GlobalDof, sol_local); + z+=znew; + } + socketstream zsock(vishost, visport); + PlotSolution(z,zsock,0); cin.get(); + } + } +} + + +void DST::GetCutOffSolution(const Vector & sol, Vector & cfsol, + int ip, Array directions, int nlayers, bool local) const +{ + + int d = directions.Size(); + int directx = directions[0]; // 1,0,-1 + int directy = directions[1]; // 1,0,-1 + int directz; + if (d ==3) directz = directions[2]; + + Mesh * mesh = ovlp_prob->fespaces[ip]->GetMesh(); + + Vector pmin, pmax; + mesh->GetBoundingBox(pmin, pmax); + double h = GetUniformMeshElementSize(povlp->patch_mesh[ip]); + Array2D pmlh(dim,2); pmlh = 0.0; + + if (directions[0]==1) + { + pmlh[0][1] = h*nlayers; + } + if (directions[0]==-1) + { + pmlh[0][0] = h*nlayers; + } + if (directions[1]==1) + { + pmlh[1][1] = h*nlayers; + } + if (directions[1]==-1) + { + pmlh[1][0] = h*nlayers; + } + + CutOffFnCoefficient cf(CutOffFncn, pmin, pmax, pmlh); + + double * data = sol.GetData(); + + FiniteElementSpace * fes; + if (!local) + { + fes = bf->FESpace(); + } + else + { + fes = ovlp_prob->fespaces[ip]; + } + + int n = fes->GetTrueVSize(); + + GridFunction solgf_re(fes, data); + GridFunction solgf_im(fes, &data[n]); + + + GridFunctionCoefficient coeff1_re(&solgf_re); + GridFunctionCoefficient coeff1_im(&solgf_im); + + ProductCoefficient prod_re(coeff1_re, cf); + ProductCoefficient prod_im(coeff1_im, cf); + + ComplexGridFunction gf(fes); + gf.ProjectCoefficient(prod_re,prod_im); + + cfsol.SetSize(sol.Size()); + cfsol = gf; +} + +void DST::GetChiRes(const Vector & res, Vector & cfres, + int ip, Array directions, int nlayers) const +{ + // int l,k; + int d = directions.Size(); + int directx = directions[0]; // 1,0,-1 + int directy = directions[1]; // 1,0,-1 + int directz; + if (d ==3) directz = directions[2]; + + Mesh * mesh = ovlp_prob->fespaces[ip]->GetMesh(); + double h = GetUniformMeshElementSize(mesh); + + Vector pmin, pmax; + mesh->GetBoundingBox(pmin, pmax); + + Array2D pmlh(dim,2); pmlh = 0.0; + + if (directions[0]==1) + { + pmlh[0][1] = h*nlayers; + } + if (directions[0]==-1) + { + pmlh[0][0] = h*nlayers; + } + if (directions[1]==1) + { + pmlh[1][1] = h*nlayers; + } + if (directions[1]==-1) + { + pmlh[1][0] = h*nlayers; + } + + CutOffFnCoefficient cf(ChiFncn, pmin, pmax, pmlh); + + double * data = res.GetData(); + + FiniteElementSpace * fespace; + fespace = ovlp_prob->fespaces[ip]; + + int n = fespace->GetTrueVSize(); + + GridFunction solgf_re(fespace, data); + GridFunction solgf_im(fespace, &data[n]); + + GridFunctionCoefficient coeff1_re(&solgf_re); + GridFunctionCoefficient coeff1_im(&solgf_im); + + ProductCoefficient prod_re(coeff1_re, cf); + ProductCoefficient prod_im(coeff1_im, cf); + + ComplexGridFunction gf(fespace); + gf.ProjectCoefficient(prod_re,prod_im); + + cfres.SetSize(res.Size()); + cfres = gf; +} + + +DST::~DST() +{ +} + + +void DST::Getijk(int ip, int & i, int & j, int & k) const +{ + k = ip/(nxyz[0]*nxyz[1]); + j = (ip-k*nxyz[0]*nxyz[1])/nxyz[0]; + i = (ip-k*nxyz[0]*nxyz[1])%nxyz[0]; +} + +int DST::GetPatchId(const Array & ijk) const +{ + int d=ijk.Size(); + int z = (dim==2)? 0 : ijk[2]; + return subdomains(ijk[0],ijk[1],z); +} + + +void DST::TransferSources(int sweep, int ip0, Vector & sol0) const +{ + // Find all neighbors of patch ip0 + int nx = nxyz[0]; + int ny = nxyz[1]; + int i0, j0, k0; + Getijk(ip0, i0,j0,k0); + // cout << "Transfer to : " << endl; + // loop through possible directions + for (int i=-1; i<2; i++) + { + int i1 = i0 + i; + if (i1 <0 || i1>=nx) continue; + for (int j=-1; j<2; j++) + { + if (i==0 && j==0) continue; + int j1 = j0 + j; + if (j1 <0 || j1>=ny) continue; + // cout << "(" << i1 << "," << j1 <<"), "; + // Find ip 1 + Array ij1(2); ij1[0] = i1; ij1[1]=j1; + int ip1 = GetPatchId(ij1); + // cout << "ip1 = " << ip1; + // cout << " in the direction of (" << i <<", " < directions(2); + directions[0] = i; + directions[1] = j; + Vector cfsol0; + GetCutOffSolution(sol0,cfsol0,ip0,directions,nrlayers,true); + + // Transfer solution to ip1; + Array * Dof2GlobalDof0 = &ovlp_prob->Dof2GlobalDof[ip0]; + Array * Dof2GlobalDof1 = &ovlp_prob->Dof2GlobalDof[ip1]; + + Vector znew(2*bf->FESpace()->GetTrueVSize()); + znew = 0.0; + znew.SetSubVector(*Dof2GlobalDof0,sol0); + + Vector sol1(Dof2GlobalDof1->Size()); sol1 = 0.0; + Vector res1(Dof2GlobalDof1->Size()); res1 = 0.0; + znew.GetSubVector(*Dof2GlobalDof1,sol1); + PmlMat[ip1]->Mult(sol1,res1); + res1 *=-1.0; + + // remove the source in the pml restrict to the non overlapping subdomain + + Array direct(2); direct = 0; + if (i1>0) direct[0] = -1; + if (j1>0) direct[1] = -1; + Vector cfraux(res1.Size()); cfraux = 0.0; + GetChiRes(res1, cfraux,ip1,direct, nrlayers); + direct = 0; + if (i1+1patch_mesh[ip]); + Array2D length(dim,2); + length = h*(nrlayers); + + CartesianPML pml(povlp->patch_mesh[ip], length); + pml.SetOmega(omega); + + Array ess_tdof_list; + if (povlp->patch_mesh[ip]->bdr_attributes.Size()) + { + Array ess_bdr(povlp->patch_mesh[ip]->bdr_attributes.Max()); + ess_bdr = 1; + ovlp_prob->fespaces[ip]->GetEssentialTrueDofs(ess_bdr, ess_tdof_list); + } + + ConstantCoefficient one(1.0); + ConstantCoefficient sigma(-pow(omega, 2)); + PmlMatrixCoefficient c1_re(dim,pml_detJ_JT_J_inv_Re,&pml); + PmlMatrixCoefficient c1_im(dim,pml_detJ_JT_J_inv_Im,&pml); + PmlCoefficient detJ_re(pml_detJ_Re,&pml); + PmlCoefficient detJ_im(pml_detJ_Im,&pml); + ProductCoefficient c2_re0(sigma, detJ_re); + ProductCoefficient c2_im0(sigma, detJ_im); + ProductCoefficient c2_re(c2_re0, *ws); + ProductCoefficient c2_im(c2_im0, *ws); + SesquilinearForm a(ovlp_prob->fespaces[ip],ComplexOperator::HERMITIAN); + + a.AddDomainIntegrator(new DiffusionIntegrator(c1_re), + new DiffusionIntegrator(c1_im)); + a.AddDomainIntegrator(new MassIntegrator(c2_re), + new MassIntegrator(c2_im)); + a.Assemble(); + + OperatorPtr Alocal; + a.FormSystemMatrix(ess_tdof_list,Alocal); + ComplexSparseMatrix * AZ_ext = Alocal.As(); + SparseMatrix * Mat = AZ_ext->GetSystemMatrix(); + Mat->Threshold(0.0); + return Mat; +} + +void DST::PlotSolution(Vector & sol, socketstream & sol_sock, int ip) const +{ + FiniteElementSpace * fespace = bf->FESpace(); + Mesh * mesh = fespace->GetMesh(); + GridFunction gf(fespace); + double * data = sol.GetData(); + gf.SetData(data); + + string keys; + if (ip == 0) keys = "keys mrRljc\n"; + // sol_sock << "solution\n" << *mesh << gf << keys << "valuerange -0.1 0.1 \n" << flush; + sol_sock << "solution\n" << *mesh << gf << keys << flush; +} \ No newline at end of file diff --git a/examples/maxwell-solver/DST/DSTnew.hpp b/examples/maxwell-solver/DST/DSTnew.hpp new file mode 100644 index 0000000000..8fbae8c9ba --- /dev/null +++ b/examples/maxwell-solver/DST/DSTnew.hpp @@ -0,0 +1,52 @@ +#pragma once +#include "Utilities.hpp" +#include "PML.hpp" +using namespace std; +using namespace mfem; + +class DST : public Solver// +{ +private: + int nrpatch; + int dim; + SesquilinearForm *bf=nullptr; + MeshPartition * povlp=nullptr; + double omega = 0.5; + Coefficient * ws; + int nrlayers; + int nxyz[3]; + const Operator * A=nullptr; + DofMap * ovlp_prob = nullptr; + Array PmlMat; + Array PmlMatInv; + Array2D Pmllength; + Array3D subdomains; + mutable Array f_orig; + int ntransf_directions; + int nsweeps; + Array2D sweeps; + Array dirx; + Array diry; + Array dirz; + mutable Array> f_transf; + Array> usol; + + SparseMatrix * GetPmlSystemMatrix(int ip); + void PlotSolution(Vector & sol, socketstream & sol_sock, int ip) const; + void GetCutOffSolution(const Vector & sol, Vector & cfsol, + int ip, Array directions, int nlayers, bool local=false) const; + void GetChiRes(const Vector & res, Vector & cfres, + int ip, Array directions, int nlayers) const; + void TransferSources(int sweep, int ip, Vector & sol_ext) const; + int GetPatchId(const Array & ijk) const; + void Getijk(int ip, int & i, int & j, int & k ) const; + int SourceTransfer(const Vector & Psi0, Array direction, int ip, Vector & Psi1) const; +public: + DST(SesquilinearForm * bf_, Array2D & Pmllength_, + double omega_, Coefficient * ws_, int nrlayers_); + virtual void SetOperator(const Operator &op) {A = &op;} + virtual void Mult(const Vector &r, Vector &z) const; + virtual ~DST(); +}; + + diff --git a/examples/maxwell-solver/DST/DiagST.cpp b/examples/maxwell-solver/DST/DiagST.cpp index ca34b0234e..d14d107b04 100644 --- a/examples/maxwell-solver/DST/DiagST.cpp +++ b/examples/maxwell-solver/DST/DiagST.cpp @@ -19,7 +19,7 @@ DiagST::DiagST(SesquilinearForm * bf_, Array2D & Pmllength_, // 1. Ovelapping partition with overlap = 2h partition_kind = 2; // Non Overlapping partition int nx=8; - int ny=8; + int ny=1; int nz=1; ovlpnrlayers = 2; povlp = new MeshPartition(mesh, partition_kind,nx,ny,nz,ovlpnrlayers); @@ -173,7 +173,7 @@ void DiagST::Mult(const Vector &r, Vector &z) const int nsteps = nx + ny - 1; // Sweep number - for (int l=0; l<1; l++) + for (int l=0; l<4; l++) { for (int s = 0; s & Pmllength_, + double omega_, Coefficient * ws_, int nrlayers_) + : Solver(2*bf_->FESpace()->GetTrueVSize(), 2*bf_->FESpace()->GetTrueVSize()), + bf(bf_), Pmllength(Pmllength_), omega(omega_), ws(ws_), nrlayers(nrlayers_) +{ + Mesh * mesh = bf->FESpace()->GetMesh(); + dim = mesh->Dimension(); + + // ----------------- Step 1 -------------------- + // Introduce 2 layered partitios of the domain + // + int partition_kind; + + // 1. Ovelapping partition with overlap = 2h + partition_kind = 2; // Non Overlapping partition + int nx=4; + int ny=1; + int nz=1; + ovlpnrlayers = nrlayers; + povlp = new MeshPartition(mesh, partition_kind,nx,ny,nz,ovlpnrlayers); + nxyz[0] = povlp->nxyz[0]; + nxyz[1] = povlp->nxyz[1]; + nxyz[2] = povlp->nxyz[2]; + nrpatch = povlp->nrpatch; + subdomains = povlp->subdomains; + + // + // ----------------- Step 1a ------------------- + // Save the partition for visualization + SaveMeshPartition(povlp->patch_mesh, "output/mesh_ovlp.", "output/sol_ovlp."); + + // // // ------------------Step 2 -------------------- + // // // Construct the dof maps from subdomains to global (for the extended and not) + + ovlp_prob = new DofMap(bf,povlp,nrlayers); + + // ------------------Step 3 -------------------- + // Assemble the PML Problem matrices and factor them + PmlMat.SetSize(nrpatch); + PmlMatInv.SetSize(nrpatch); + for (int ip=0; ipSetOperator(*PmlMat[ip]); + } + + + nsweeps = pow(2,dim); + sweeps.SetSize(nsweeps,dim); + // 2D + sweeps(0,0) = 1; sweeps(0,1) = 1; + sweeps(1,0) = -1; sweeps(1,1) = 1; + sweeps(2,0) = 1; sweeps(2,1) =-1; + sweeps(3,0) = -1; sweeps(3,1) =-1; + + // Set up src arrays size + f_orig.SetSize(nrpatch); + f_transf.SetSize(nrpatch); + + + // Construct a simple map used for directions of transfer + ConstructDirectionsMap(); + for (int ip=0; ipfespaces[ip]->GetTrueVSize(); // (x 2 for complex ) + int npml = 2*ovlp_prob->PmlFespaces[ip]->GetTrueVSize(); // (x 2 for complex ) + f_orig[ip] = new Vector(n); *f_orig[ip] = 0.0; + f_transf[ip].SetSize(nsweeps); + for (int i=0;iPmlMeshes[ip]); + Array2D length(dim,2); + length = h*(nrlayers); + + CartesianPML pml(ovlp_prob->PmlMeshes[ip], length); + pml.SetOmega(omega); + + Array ess_tdof_list; + if (ovlp_prob->PmlMeshes[ip]->bdr_attributes.Size()) + { + Array ess_bdr(ovlp_prob->PmlMeshes[ip]->bdr_attributes.Max()); + ess_bdr = 1; + ovlp_prob->PmlFespaces[ip]->GetEssentialTrueDofs(ess_bdr, ess_tdof_list); + } + + ConstantCoefficient one(1.0); + ConstantCoefficient sigma(-pow(omega, 2)); + + PmlMatrixCoefficient c1_re(dim,pml_detJ_JT_J_inv_Re,&pml); + PmlMatrixCoefficient c1_im(dim,pml_detJ_JT_J_inv_Im,&pml); + + PmlCoefficient detJ_re(pml_detJ_Re,&pml); + PmlCoefficient detJ_im(pml_detJ_Im,&pml); + + ProductCoefficient c2_re0(sigma, detJ_re); + ProductCoefficient c2_im0(sigma, detJ_im); + + ProductCoefficient c2_re(c2_re0, *ws); + ProductCoefficient c2_im(c2_im0, *ws); + + SesquilinearForm a(ovlp_prob->PmlFespaces[ip],ComplexOperator::HERMITIAN); + + a.AddDomainIntegrator(new DiffusionIntegrator(c1_re), + new DiffusionIntegrator(c1_im)); + a.AddDomainIntegrator(new MassIntegrator(c2_re), + new MassIntegrator(c2_im)); + a.Assemble(); + + OperatorPtr Alocal; + a.FormSystemMatrix(ess_tdof_list,Alocal); + ComplexSparseMatrix * AZ_ext = Alocal.As(); + SparseMatrix * Mat = AZ_ext->GetSystemMatrix(); + Mat->Threshold(0.0); + return Mat; +} + + +void DiagST::Mult(const Vector &r, Vector &z) const +{ + // Step 0 + // Restrict original sources to the patches + for (int ip=0; ip * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; + r.GetSubVector(*Dof2GlobalDof,*f_orig[ip]); + } + + char vishost[] = "localhost"; + int visport = 19916; + z = 0.0; + Vector znew(z); + Vector rnew(r); + Vector raux(r); raux = 0.0; + Vector z1(z); + Vector z2(z); + Vector z3(z); + Vector z4(z); + znew = 0.0; + + // in 2D there are a total of 4 sweeps + // with nx + ny - 1 serial steps each + // -------------------------------------------- + // Sweep in the direction (1,1) + // -------------------------------------------- + int nx = nxyz[0]; + int ny = nxyz[1]; + + int nsteps = nx + ny - 1; + // Sweep number + for (int l=0; l<1; l++) + { + for (int s = 0; s=ny) continue; + + // find patch id + Array ij(2); ij[0] = i; ij[1]=j; + int ip = GetPatchId(ij); + + // Solve the PML problem in patch ip with all sources + // Original and all transfered (maybe some of them) + Array * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; + Array * Dof2PmlDof = &ovlp_prob->Dof2PmlDof[ip]; + int ndofs = Dof2GlobalDof->Size(); + Vector sol_local(ndofs); sol_local = 0.0; + Vector res_local(ndofs); res_local = 0.0; + if (l==0) res_local += *f_orig[ip]; + + res_local += *f_transf[ip][l]; + // if (res_local.Norml2() < 1e-12) continue; + // Extend by zero to the PML mesh + int nrdof_ext = PmlMat[ip]->Height(); + Vector res_ext(nrdof_ext); res_ext = 0.0; + Vector sol_ext(nrdof_ext); sol_ext = 0.0; + + res_ext.SetSubVector(*Dof2PmlDof,res_local); + PmlMatInv[ip]->Mult(res_ext, sol_ext); + + // socketstream resipsock(vishost, visport); + // PlotSolution(res_ext,resipsock,ip,true, true); cin.get(); + + // socketstream solipsock(vishost, visport); + // PlotSolution(sol_ext,solipsock,ip,true, true); cin.get(); + + // Multiply with the cutoff functions, find the new sources and + // and propagate to all neighboring subdomains + // (possible 8 in 2D, 26 in 3D) + // socketstream lsock(vishost, visport); + // PlotSolution(sol_ext,lsock,ip,true, true); cin.get(); + TransferSources(l,ip, sol_ext); + Vector cfsol_ext(sol_ext.Size()); + + // cut off the ip solution to all possible directions + Arraydirections(2); directions = 0; + + if (i+10) directions[0] = -1; + if (j>0) directions[1] = -1; + GetCutOffSolution(sol_ext,cfsol_ext,ip,directions,ovlpnrlayers,true); + cfsol_ext.GetSubVector(*Dof2PmlDof, sol_local); + znew = 0.0; + znew.SetSubVector(*Dof2GlobalDof, sol_local); + z+=znew; + } + socketstream zsock(vishost, visport); + PlotSolution(z,zsock,0); cin.get(); + } + } +} + +void DiagST::PlotSolution(Vector & sol, socketstream & sol_sock, int ip, + bool localdomain, bool pmldomain) const +{ + FiniteElementSpace * fes; + if (!localdomain) + { + fes = bf->FESpace(); + } + else if (!pmldomain) + { + fes = ovlp_prob->fespaces[ip]; + } + else + { + fes = ovlp_prob->PmlFespaces[ip]; + } + Mesh * mesh = fes->GetMesh(); + GridFunction gf(fes); + double * data = sol.GetData(); + gf.SetData(data); + + string keys; + if (ip == 0) keys = "keys mrRljc\n"; + // sol_sock << "solution\n" << *mesh << gf << keys << "valuerange -0.1 0.1 \n" << flush; + sol_sock << "solution\n" << *mesh << gf << keys << flush; +} + +void DiagST::GetCutOffSolution(const Vector & sol, Vector & cfsol, + int ip, Array directions, int ovlpnlayers, bool local) const +{ + // int l,k; + int d = directions.Size(); + int directx = directions[0]; // 1,0,-1 + int directy = directions[1]; // 1,0,-1 + int directz; + if (d ==3) directz = directions[2]; + + Mesh * mesh = ovlp_prob->fespaces[ip]->GetMesh(); + + Vector pmin, pmax; + mesh->GetBoundingBox(pmin, pmax); + double h = GetUniformMeshElementSize(mesh); + + Array2D pmlh(dim,2); pmlh = 0.0; + + if (directions[0]==1) + { + pmlh[0][1] = h*ovlpnlayers; + } + if (directions[0]==-1) + { + pmlh[0][0] = h*ovlpnlayers; + } + if (directions[1]==1) + { + pmlh[1][1] = h*ovlpnlayers; + } + if (directions[1]==-1) + { + pmlh[1][0] = h*ovlpnlayers; + } + + CutOffFnCoefficient cf(CutOffFncn, pmin, pmax, pmlh); + + double * data = sol.GetData(); + + FiniteElementSpace * fespace; + if (!local) + { + fespace = bf->FESpace(); + } + else + { + fespace = ovlp_prob->PmlFespaces[ip]; + } + + int n = fespace->GetTrueVSize(); + + GridFunction solgf_re(fespace, data); + GridFunction solgf_im(fespace, &data[n]); + + GridFunctionCoefficient coeff1_re(&solgf_re); + GridFunctionCoefficient coeff1_im(&solgf_im); + + ProductCoefficient prod_re(coeff1_re, cf); + ProductCoefficient prod_im(coeff1_im, cf); + + ComplexGridFunction gf(fespace); + gf.ProjectCoefficient(prod_re,prod_im); + + cfsol.SetSize(sol.Size()); + cfsol = gf; +} + + +void DiagST::GetChiRes(const Vector & res, Vector & cfres, + int ip, Array directions, int nlayers) const +{ + // int l,k; + int d = directions.Size(); + int directx = directions[0]; // 1,0,-1 + int directy = directions[1]; // 1,0,-1 + int directz; + if (d ==3) directz = directions[2]; + + Mesh * mesh = ovlp_prob->fespaces[ip]->GetMesh(); + double h = GetUniformMeshElementSize(mesh); + + Vector pmin, pmax; + mesh->GetBoundingBox(pmin, pmax); + + Array2D pmlh(dim,2); pmlh = 0.0; + + if (directions[0]==1) + { + pmlh[0][1] = h*nlayers; + } + if (directions[0]==-1) + { + pmlh[0][0] = h*nlayers; + } + if (directions[1]==1) + { + pmlh[1][1] = h*nlayers; + } + if (directions[1]==-1) + { + pmlh[1][0] = h*nlayers; + } + + CutOffFnCoefficient cf(ChiFncn, pmin, pmax, pmlh); + + double * data = res.GetData(); + + FiniteElementSpace * fespace; + fespace = ovlp_prob->fespaces[ip]; + + int n = fespace->GetTrueVSize(); + + GridFunction solgf_re(fespace, data); + GridFunction solgf_im(fespace, &data[n]); + + GridFunctionCoefficient coeff1_re(&solgf_re); + GridFunctionCoefficient coeff1_im(&solgf_im); + + ProductCoefficient prod_re(coeff1_re, cf); + ProductCoefficient prod_im(coeff1_im, cf); + + ComplexGridFunction gf(fespace); + gf.ProjectCoefficient(prod_re,prod_im); + + cfres.SetSize(res.Size()); + cfres = gf; +} + +DiagST::~DiagST() +{ + for (int ip = 0; ip & ijk) const +{ + int d=ijk.Size(); + if (d==2) + { + return subdomains(ijk[0],ijk[1],0); + } + else + { + return subdomains(ijk[0],ijk[1],ijk[2]); + } +} + +int DiagST::SourceTransfer(const Vector & Psi0, Array direction, int ip0, Vector & Psi1) const +{ + // For now 2D problems only + // Directions + // direction (1,1) + int i0,j0,k0; + Getijk(ip0,i0,j0,k0); + + int i1 = i0+direction[0]; + int j1 = j0+direction[1]; + Array ij(2); ij[0]=i1; ij[1]=j1; + int ip1 = GetPatchId(ij); + + MFEM_VERIFY(i1 < nxyz[0] && i1>=0, "SourceTransfer: i1 out of bounds"); + MFEM_VERIFY(j1 < nxyz[1] && j1>=0, "SourceTransfer: j1 out of bounds"); + + Array * Dof2GlobalDof0 = &ovlp_prob->Dof2GlobalDof[ip0]; + Array * Dof2GlobalDof1 = &ovlp_prob->Dof2GlobalDof[ip1]; + Psi1.SetSize(Dof2GlobalDof1->Size()); Psi1=0.0; + Vector r(2*bf->FESpace()->GetTrueVSize()); + r = 0.0; + r.SetSubVector(*Dof2GlobalDof0,Psi0); + + Vector zloc(Psi1.Size()); + r.GetSubVector(*Dof2GlobalDof1,zloc); + // extend + Vector psi_ext(PmlMat[ip1]->Height()); + Vector zloc_ext(PmlMat[ip1]->Height()); + Array * Dof2PmlDof1 = &ovlp_prob->Dof2PmlDof[ip1]; + zloc_ext.SetSubVector(*Dof2PmlDof1,zloc); + PmlMat[ip1]->Mult(zloc_ext,psi_ext); + psi_ext *=-1.0; + Vector Psi(Dof2GlobalDof1->Size()); Psi=0.0; + + psi_ext.GetSubVector(*Dof2PmlDof1,Psi); + + Psi1 = Psi; + + int nx = nxyz[0]; + int ny = nxyz[1]; + Array direct(2); direct = 0; + if (i1>0) direct[0] = -1; + if (j1>0) direct[1] = -1; + GetChiRes(Psi, Psi1,ip1,direct, ovlpnrlayers); + Psi = Psi1; + direct = 0; + if (i1+1 ijk(dim); + if (dim==2) + { + for (int i=-1; i<=1; i++) // directions x + { + for (int j=-1; j<=1; j++) // directions y + { + ijk[0]=i; + ijk[1]=j; + int k=GetDirectionId(ijk); + dirx[k]=i; + diry[k]=j; + } + } + } + else if (dim==3) + { + dirz.SetSize(ntransf_directions); + for (int i=-1; i<=1; i++) // directions x + { + for (int j=-1; j<=1; j++) // directions y + { + for (int k=-1; k<=1; k++) // directions zß + { + ijk[0]=i; + ijk[1]=j; + ijk[2]=k; + int l=GetDirectionId(ijk); + dirx[l]=i; + diry[l]=j; + dirz[l]=k; + } + } + } + } + + // cout << "dirx = " << endl; + // dirx.Print(cout,ntransf_directions); + // cout << "diry = " << endl; + // diry.Print(cout,ntransf_directions); + + // if (dim==2) + // { + // for (int id=0; id<9; id++) + // { + // GetDirectionijk(id,ijk); + // // cout << "for id = " << id << ": (" < & ijk) const +{ + int d = ijk.Size(); + int n=3; + if (d==2) + { + return (ijk[0]+1)*n+(ijk[1]+1); + } + else + { + return (ijk[0]+1)*n*n+(ijk[1]+1)*n+ijk[2]+1; + } +} + +void DiagST::GetDirectionijk(int id, Array & ijk) const +{ + int d = ijk.Size(); + int n=3; + if (d==2) + { + ijk[0]=id/n - 1; + ijk[1]=id%n - 1; + } + else + { + ijk[0]=id/(n*n)-1; + ijk[1]=(id-(ijk[0]+1)*n*n)/n - 1; + ijk[2]=(id-(ijk[0]+1)*n*n)%n - 1; + } + // cout << "ijk = " ; ijk.Print(); +} + + + + +void DiagST::TransferSources(int sweep, int ip0, Vector & sol_ext) const +{ + // Find all neighbors of patch ip + int nx = nxyz[0]; + int ny = nxyz[1]; + int i0, j0, k0; + Getijk(ip0, i0,j0,k0); + // cout << "Transfer to : " << endl; + // loop through possible directions + for (int i=-1; i<2; i++) + { + int i1 = i0 + i; + if (i1 <0 || i1>=nx) continue; + for (int j=-1; j<2; j++) + { + if (i==0 && j==0) continue; + int j1 = j0 + j; + if (j1 <0 || j1>=ny) continue; + // cout << "(" << i1 << "," << j1 <<"), "; + // Find ip 1 + Array ij1(2); ij1[0] = i1; ij1[1]=j1; + int ip1 = GetPatchId(ij1); + // cout << "ip1 = " << ip1; + // cout << " in the direction of (" << i <<", " < directions(2); + directions[0] = i; + directions[1] = j; + Vector cfsol_ext; + Vector res_ext(sol_ext.Size()); + GetCutOffSolution(sol_ext,cfsol_ext,ip0,directions,ovlpnrlayers,true); + // sol_ext = cfsol_ext; + // PmlMat[ip0]->Mult(cfsol_ext, res_ext); + + //--------------------------------------- + char vishost[] = "localhost"; + int visport = 19916; + // socketstream pmlsock(vishost, visport); + // PlotSolution(res_ext,pmlsock,ip0,true,true); cin.get(); + //--------------------------------------- + + // res_ext*= -1.0; + Array *Dof2PmlDof = &ovlp_prob->Dof2PmlDof[ip0]; + // Vector res_local(Dof2PmlDof->Size()); res_local = 0.0; + // res_ext.GetSubVector(*Dof2PmlDof,res_local); + + //----------------------------- + // pass to ip1 and calculate residual there + Vector sol_local(Dof2PmlDof->Size()); sol_local = 0.0; + cfsol_ext.GetSubVector(*Dof2PmlDof,sol_local); + //----------------------------- + + + + + + // Find the minumum sweep number that to transfer the source that + // satisfies the two rules + for (int l=sweep; lSize()==raux.Size(), + "Transfer Sources: inconsistent size"); + *f_transf[ip1][l]+=raux; + + break; + } + + } + } +} + + + + + + + + + + +// //Diagonal Source Transfer Preconditioner + +// #include "DiagST.hpp" + +// DiagST::DiagST(SesquilinearForm * bf_, Array2D & Pmllength_, +// double omega_, Coefficient * ws_, int nrlayers_) +// : Solver(2*bf_->FESpace()->GetTrueVSize(), 2*bf_->FESpace()->GetTrueVSize()), +// bf(bf_), Pmllength(Pmllength_), omega(omega_), ws(ws_), nrlayers(nrlayers_) +// { +// Mesh * mesh = bf->FESpace()->GetMesh(); +// dim = mesh->Dimension(); + +// // ----------------- Step 1 -------------------- +// // Introduce 2 layered partitios of the domain +// // +// int partition_kind; + +// // 1. Ovelapping partition with overlap = 2h +// partition_kind = 2; // Non Overlapping partition +// int nx=2; +// int ny=2; +// int nz=1; +// povlp = new MeshPartition(mesh, partition_kind,nx,ny,nz,2); +// nxyz[0] = povlp->nxyz[0]; +// nxyz[1] = povlp->nxyz[1]; +// nxyz[2] = povlp->nxyz[2]; +// nrpatch = povlp->nrpatch; +// subdomains = povlp->subdomains; + +// // +// // ----------------- Step 1a ------------------- +// // Save the partition for visualization +// // SaveMeshPartition(povlp->patch_mesh, "output/mesh_ovlp.", "output/sol_ovlp."); + +// // // // ------------------Step 2 -------------------- +// // // // Construct the dof maps from subdomains to global (for the extended and not) + +// ovlp_prob = new DofMap(bf,povlp,nrlayers); + +// // ------------------Step 3 -------------------- +// // Assemble the PML Problem matrices and factor them +// PmlMat.SetSize(nrpatch); +// PmlMatInv.SetSize(nrpatch); +// for (int ip=0; ipSetOperator(*PmlMat[ip]); +// } + + +// nsweeps = pow(2,dim); +// sweeps.SetSize(nsweeps,dim); +// // 2D +// sweeps(0,0) = 1; sweeps(0,1) = 1; +// sweeps(1,0) = -1; sweeps(1,1) = 1; +// sweeps(2,0) = 1; sweeps(2,1) =-1; +// sweeps(3,0) = -1; sweeps(3,1) =-1; + +// // Set up src arrays size +// f_orig.SetSize(nrpatch); +// f_transf.SetSize(nrpatch); +// usol.SetSize(nrpatch); + + +// // Construct a simple map used for directions of transfer +// ConstructDirectionsMap(); +// for (int ip=0; ipfespaces[ip]->GetTrueVSize(); // (x 2 for complex ) +// int npml = 2*ovlp_prob->PmlFespaces[ip]->GetTrueVSize(); // (x 2 for complex ) +// f_orig[ip] = new Vector(n); *f_orig[ip] = 0.0; +// f_transf[ip].SetSize(nsweeps); +// usol[ip].SetSize(nsweeps); +// for (int i=0;iPmlMeshes[ip]); +// Array2D length(dim,2); +// length = h*(nrlayers); + +// CartesianPML pml(ovlp_prob->PmlMeshes[ip], length); +// pml.SetOmega(omega); + +// Array ess_tdof_list; +// if (ovlp_prob->PmlMeshes[ip]->bdr_attributes.Size()) +// { +// Array ess_bdr(ovlp_prob->PmlMeshes[ip]->bdr_attributes.Max()); +// ess_bdr = 1; +// ovlp_prob->PmlFespaces[ip]->GetEssentialTrueDofs(ess_bdr, ess_tdof_list); +// } + +// ConstantCoefficient one(1.0); +// ConstantCoefficient sigma(-pow(omega, 2)); + +// PmlMatrixCoefficient c1_re(dim,pml_detJ_JT_J_inv_Re,&pml); +// PmlMatrixCoefficient c1_im(dim,pml_detJ_JT_J_inv_Im,&pml); + +// PmlCoefficient detJ_re(pml_detJ_Re,&pml); +// PmlCoefficient detJ_im(pml_detJ_Im,&pml); + +// ProductCoefficient c2_re0(sigma, detJ_re); +// ProductCoefficient c2_im0(sigma, detJ_im); + +// ProductCoefficient c2_re(c2_re0, *ws); +// ProductCoefficient c2_im(c2_im0, *ws); + +// SesquilinearForm a(ovlp_prob->PmlFespaces[ip],ComplexOperator::HERMITIAN); + +// a.AddDomainIntegrator(new DiffusionIntegrator(c1_re), +// new DiffusionIntegrator(c1_im)); +// a.AddDomainIntegrator(new MassIntegrator(c2_re), +// new MassIntegrator(c2_im)); +// a.Assemble(); + +// OperatorPtr Alocal; +// a.FormSystemMatrix(ess_tdof_list,Alocal); +// ComplexSparseMatrix * AZ_ext = Alocal.As(); +// SparseMatrix * Mat = AZ_ext->GetSystemMatrix(); +// Mat->Threshold(0.0); +// return Mat; +// } + + +// void DiagST::Mult(const Vector &r, Vector &z) const +// { +// // Step 0 +// // Restrict original sources to the patches +// for (int ip=0; ip * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; +// r.GetSubVector(*Dof2GlobalDof,*f_orig[ip]); +// } + +// char vishost[] = "localhost"; +// int visport = 19916; +// z = 0.0; +// Vector znew(z); +// Vector rnew(r); +// Vector raux(r); raux = 0.0; +// Vector z1(z); +// Vector z2(z); +// Vector z3(z); +// Vector z4(z); +// znew = 0.0; + +// // in 2D there are a total of 4 sweeps +// // with nx + ny - 1 serial steps each +// // -------------------------------------------- +// // Sweep in the direction (1,1) +// // -------------------------------------------- +// int nx = nxyz[0]; +// int ny = nxyz[1]; + +// int nsteps = nx + ny - 1; +// // Sweep number +// for (int l=0; l<1; l++) +// { +// for (int s = 0; s=ny) continue; + +// // find patch id +// Array ij(2); ij[0] = i; ij[1]=j; +// int ip = GetPatchId(ij); + +// // Solve the PML problem in patch ip with all sources +// // Original and all transfered (maybe some of them) +// Array * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; +// Array * Dof2PmlDof = &ovlp_prob->Dof2PmlDof[ip]; +// int ndofs = Dof2GlobalDof->Size(); +// Vector sol_local(ndofs); +// Vector res_local(ndofs); +// res_local = *f_orig[ip]; + +// res_local += *f_transf[ip][l]; +// // Extend by zero to the PML mesh +// int nrdof_ext = PmlMat[ip]->Height(); + +// Vector res_ext(nrdof_ext); res_ext = 0.0; +// Vector sol_ext(nrdof_ext); sol_ext = 0.0; + +// res_ext.SetSubVector(*Dof2PmlDof,res_local); +// PmlMatInv[ip]->Mult(res_ext, sol_ext); +// // *usol[l][ip] = sol_ext; +// // Multiply with the cutoff functions, find the new sources and +// // and propagate to all neighboring subdomains +// // (possible 8 in 2D, 26 in 3D) +// TransferSources(l,ip, sol_ext); +// Vector cfsol_ext(sol_ext.Size()); + +// // cut off the ip solution to all possible directions +// Arraydirections(2); directions = 0; +// if (i+10) directions[0] = -1; +// // if (j>0) directions[1] = -1; +// // GetCutOffSolution(sol_ext,cfsol_ext,ip,directions,true); +// cfsol_ext.GetSubVector(*Dof2PmlDof, sol_local); +// znew = 0.0; +// znew.SetSubVector(*Dof2GlobalDof, sol_local); +// // z1.AddElementVector(*Dof2GlobalDof, sol_local); +// z1+=znew; +// socketstream sub1_sock1(vishost, visport); +// PlotSolution(z1,sub1_sock1,0); cin.get(); +// } +// } +// } + +// // PlotSolution(z1,sub_sock1, 0); +// // cin.get(); + +// z +=z1; +// // A->Mult(z,raux); rnew = r; rnew -=raux; + +// // for (int ip=0; ip * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; +// // rnew.GetSubVector(*Dof2GlobalDof,*f_orig[ip]); +// // } + + +// // for (int l=1; l<2; l++) +// // { +// // for (int s = 0; s=ny) continue; +// // // cout << "2:Patch no: (" << i <<"," << j << ")" << endl; +// // // cin.get(); +// // // find patch id +// // Array ij(2); ij[0] = i; ij[1]=j; +// // int ip = GetPatchId(ij); +// // // cout << "ip = " << ip << endl; + +// // // Solve the PML problem in patch ip with all sources +// // // Original and all transfered (maybe some of them) +// // Array * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; +// // Array * Dof2PmlDof = &ovlp_prob->Dof2PmlDof[ip]; +// // int ndofs = Dof2GlobalDof->Size(); +// // Vector sol_local(ndofs); +// // Vector res_local(ndofs); +// // res_local = *f_orig[ip]; +// // // res_local = 0.0; +// // res_local += *f_transf[ip][l]; +// // // Extend by zero to the PML mesh +// // int nrdof_ext = PmlMat[ip]->Height(); + +// // Vector res_ext(nrdof_ext); res_ext = 0.0; +// // Vector sol_ext(nrdof_ext); sol_ext = 0.0; + +// // res_ext.SetSubVector(*Dof2PmlDof,res_local); +// // PmlMatInv[ip]->Mult(res_ext, sol_ext); + +// // // Multiply with the cutoff functions, find the new sources and +// // // and propagate to all neighboring subdomains +// // // (possible 8 in 2D, 26 in 3D) +// // TransferSources(l,ip, sol_ext); +// // Vector cfsol_ext(sol_ext.Size()); + +// // // cut off the ip solution to all possible directions +// // Arraydirections(2); directions = 0; +// // if (i>0) directions[0] = -1; +// // if (j+10) directions[1] = -1; +// // // GetCutOffSolution(sol_ext,cfsol_ext,ip,directions,true); + +// // cfsol_ext.GetSubVector(*Dof2PmlDof, sol_local); +// // znew = 0.0; +// // znew.SetSubVector(*Dof2GlobalDof, sol_local); +// // // z2.AddElementVector(*Dof2GlobalDof, sol_local); +// // z2+=znew; +// // } +// // } +// // } + +// // // // // socketstream sub_sock2(vishost, visport); +// // // // // PlotSolution(z2,sub_sock2, 0); +// // // // // cin.get(); + + +// // z +=z2; +// // A->Mult(z,raux); rnew = r; rnew -=raux; +// // for (int ip=0; ip * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; +// // rnew.GetSubVector(*Dof2GlobalDof,*f_orig[ip]); +// // } + +// // for (int l=2; l<3; l++) +// // { +// // for (int s = 0; s=ny) continue; +// // // cout << "2:Patch no: (" << i <<"," << j << ")" << endl; +// // // cin.get(); +// // // find patch id +// // Array ij(2); ij[0] = i; ij[1]=j; +// // int ip = GetPatchId(ij); +// // // cout << "ip = " << ip << endl; + +// // // Solve the PML problem in patch ip with all sources +// // // Original and all transfered (maybe some of them) +// // Array * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; +// // Array * Dof2PmlDof = &ovlp_prob->Dof2PmlDof[ip]; +// // int ndofs = Dof2GlobalDof->Size(); +// // Vector sol_local(ndofs); +// // Vector res_local(ndofs); +// // res_local = *f_orig[ip]; +// // // res_local = 0.0; +// // res_local += *f_transf[ip][l]; +// // // Extend by zero to the PML mesh +// // int nrdof_ext = PmlMat[ip]->Height(); + +// // Vector res_ext(nrdof_ext); res_ext = 0.0; +// // Vector sol_ext(nrdof_ext); sol_ext = 0.0; + +// // res_ext.SetSubVector(*Dof2PmlDof,res_local); +// // PmlMatInv[ip]->Mult(res_ext, sol_ext); + +// // // Multiply with the cutoff functions, find the new sources and +// // // and propagate to all neighboring subdomains +// // // (possible 8 in 2D, 26 in 3D) +// // TransferSources(l,ip, sol_ext); +// // Vector cfsol_ext(sol_ext.Size()); + +// // // cut off the ip solution to all possible directions +// // Arraydirections(2); directions = 0; +// // if (i+10) directions[1] = -1; +// // GetCutOffSolution(sol_ext,cfsol_ext,ip,directions,true); +// // // sol_ext = cfsol_ext; +// // // directions = 0; +// // // if (i>0) directions[0] = -1; +// // // if (j+1Mult(z,raux); rnew=r; rnew -=raux; + +// // for (int ip=0; ip * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; +// // rnew.GetSubVector(*Dof2GlobalDof,*f_orig[ip]); +// // } + + +// // for (int l=3; l<4; l++) +// // { +// // for (int s = 0; s=ny) continue; +// // // cout << "2:Patch no: (" << i <<"," << j << ")" << endl; +// // // cin.get(); +// // // find patch id +// // Array ij(2); ij[0] = i; ij[1]=j; +// // int ip = GetPatchId(ij); +// // // cout << "ip = " << ip << endl; + +// // // Solve the PML problem in patch ip with all sources +// // // Original and all transfered (maybe some of them) +// // Array * Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip]; +// // Array * Dof2PmlDof = &ovlp_prob->Dof2PmlDof[ip]; +// // int ndofs = Dof2GlobalDof->Size(); +// // Vector sol_local(ndofs); +// // Vector res_local(ndofs); +// // res_local = *f_orig[ip]; +// // // res_local = 0.0; +// // res_local += *f_transf[ip][l]; +// // // Extend by zero to the PML mesh +// // int nrdof_ext = PmlMat[ip]->Height(); + +// // Vector res_ext(nrdof_ext); res_ext = 0.0; +// // Vector sol_ext(nrdof_ext); sol_ext = 0.0; + +// // res_ext.SetSubVector(*Dof2PmlDof,res_local); +// // PmlMatInv[ip]->Mult(res_ext, sol_ext); + +// // // Multiply with the cutoff functions, find the new sources and +// // // and propagate to all neighboring subdomains +// // // (possible 8 in 2D, 26 in 3D) +// // TransferSources(l,ip, sol_ext); +// // Vector cfsol_ext(sol_ext.Size()); + +// // // cut off the ip solution to all possible directions +// // Arraydirections(2); directions = 0; +// // if (i>0) directions[0] = -1; +// // if (j>0) directions[1] = -1; +// // GetCutOffSolution(sol_ext,cfsol_ext,ip,directions,true); +// // // sol_ext = cfsol_ext; +// // // directions = 0; +// // // if (i+1FESpace(); +// Mesh * mesh = fespace->GetMesh(); +// GridFunction gf(fespace); +// double * data = sol.GetData(); +// // gf.SetData(&data[fespace->GetTrueVSize()]); +// gf.SetData(data); + +// string keys; +// if (ip == 0) keys = "keys mrRljc\n"; +// sol_sock << "solution\n" << *mesh << gf << keys << "valuerange -0.1 0.1 \n" << flush; +// // sol_sock << "solution\n" << *mesh << gf << keys << flush; +// } + +// void DiagST::GetCutOffSolution(const Vector & sol, Vector & cfsol, +// int ip0, Array directions, bool local) const +// { +// // int l,k; +// int d = directions.Size(); +// int directx = directions[0]; // 1,0,-1 +// int directy = directions[1]; // 1,0,-1 +// int directz; +// if (d ==3) directz = directions[2]; + +// // cout << "ip0 = " << ip0 << endl; + +// int i0, j0, k0; +// Getijk(ip0,i0, j0, k0); +// // cout << "(i0,j0) = " << "(" <=0, "GetCutOffSolution: i1 out of bounds"); +// MFEM_VERIFY(j1 < nxyz[1] && j1>=0, "GetCutOffSolution: j1 out of bounds"); + +// Array ijk(d); +// ijk[0] = i1; +// ijk[1] = j1; +// int ip1 = GetPatchId(ijk); + +// // cout << "ip1 = " << ip1 << endl; +// // cout << "(i1,j1) = " << "(" << i1 <<","<fespaces[ip0]->GetMesh(); +// Mesh * mesh1 = ovlp_prob->fespaces[ip1]->GetMesh(); + +// Vector pmin0, pmax0; +// Vector pmin1, pmax1; +// mesh0->GetBoundingBox(pmin0, pmax0); +// mesh1->GetBoundingBox(pmin1, pmax1); + +// Array2D h(dim,2); h = 0.0; + +// if (directions[0]==1) +// { +// h[0][1] = pmax0[0] - pmin1[0]; +// } +// if (directions[0]==-1) +// { +// h[0][0] = pmax1[0] - pmin0[0]; +// } +// if (directions[1]==1) +// { +// h[1][1] = pmax0[1] - pmin1[1]; +// } +// if (directions[1]==-1) +// { +// h[1][0] = pmax1[1] - pmin0[1]; +// } + +// pmin0.Print(); +// pmax0.Print(); +// h.Print(); + +// CutOffFnCoefficient cf(CutOffFncn, pmin0, pmax0, h); + +// double * data = sol.GetData(); + +// FiniteElementSpace * fespace; +// if (!local) +// { +// fespace = bf->FESpace(); +// } +// else +// { +// fespace = ovlp_prob->PmlFespaces[ip0]; +// } + +// int n = fespace->GetTrueVSize(); +// // GridFunction cutF(fespace); +// // cutF.ProjectCoefficient(cf); +// // char vishost[] = "localhost"; +// // int visport = 19916; + +// // socketstream sub_sock1(vishost, visport); +// // sub_sock1 << "solution\n" << *fespace->GetMesh() << cutF << flush; +// // cin.get(); + + +// GridFunction solgf_re(fespace, data); +// GridFunction solgf_im(fespace, &data[n]); + +// // socketstream sub_sock(vishost, visport); +// // sub_sock << "solution\n" << *fespace->GetMesh() << solgf_re << flush; +// // cin.get(); + + +// GridFunctionCoefficient coeff1_re(&solgf_re); +// GridFunctionCoefficient coeff1_im(&solgf_im); + +// ProductCoefficient prod_re(coeff1_re, cf); +// ProductCoefficient prod_im(coeff1_im, cf); + +// ComplexGridFunction gf(fespace); +// gf.ProjectCoefficient(prod_re,prod_im); + +// cfsol.SetSize(sol.Size()); +// cfsol = gf; +// // socketstream sub_sock2(vishost, visport); +// // sub_sock2 << "solution\n" << *fespace->GetMesh() << gf.real() << flush; +// // cin.get(); +// } + + +// void DiagST::GetChiRes(const Vector & res, Vector & cfres, +// int ip, Array directions, int nlayers) const +// { +// // int l,k; +// int d = directions.Size(); +// int directx = directions[0]; // 1,0,-1 +// int directy = directions[1]; // 1,0,-1 +// int directz; +// if (d ==3) directz = directions[2]; + +// Mesh * mesh = ovlp_prob->fespaces[ip]->GetMesh(); +// double h = GetUniformMeshElementSize(mesh); + +// Vector pmin, pmax; +// mesh->GetBoundingBox(pmin, pmax); + +// Array2D pmlh(dim,2); pmlh = 0.0; + +// if (directions[0]==1) +// { +// pmlh[0][1] = h*nlayers; +// } +// if (directions[0]==-1) +// { +// pmlh[0][0] = h*nlayers; +// } +// if (directions[1]==1) +// { +// pmlh[1][1] = h*nlayers; +// } +// if (directions[1]==-1) +// { +// pmlh[1][0] = h*nlayers; +// } + +// CutOffFnCoefficient cf(CutOffFncn, pmin, pmax, pmlh); + +// double * data = res.GetData(); + +// FiniteElementSpace * fespace; +// fespace = ovlp_prob->PmlFespaces[ip]; + +// int n = fespace->GetTrueVSize(); + +// GridFunction solgf_re(fespace, data); +// GridFunction solgf_im(fespace, &data[n]); + +// GridFunctionCoefficient coeff1_re(&solgf_re); +// GridFunctionCoefficient coeff1_im(&solgf_im); + +// ProductCoefficient prod_re(coeff1_re, cf); +// ProductCoefficient prod_im(coeff1_im, cf); + +// ComplexGridFunction gf(fespace); +// gf.ProjectCoefficient(prod_re,prod_im); + +// cfres.SetSize(res.Size()); +// cfres = gf; +// } + +// DiagST::~DiagST() +// { +// for (int ip = 0; ip & ijk) const +// { +// int d=ijk.Size(); +// if (d==2) +// { +// return subdomains(ijk[0],ijk[1],0); +// } +// else +// { +// return subdomains(ijk[0],ijk[1],ijk[2]); +// } +// } + +// int DiagST::SourceTransfer(const Vector & Psi0, Array direction, int ip0, Vector & Psi1) const +// { +// // For now 2D problems only +// // Directions +// // direction (1,1) +// int i0,j0,k0; +// Getijk(ip0,i0,j0,k0); + +// int i1 = i0+direction[0]; +// int j1 = j0+direction[1]; +// Array ij(2); ij[0]=i1; ij[1]=j1; +// int ip1 = GetPatchId(ij); + +// MFEM_VERIFY(i1 < nxyz[0] && i1>=0, "SourceTransfer: i1 out of bounds"); +// MFEM_VERIFY(j1 < nxyz[1] && j1>=0, "SourceTransfer: j1 out of bounds"); + +// Array * Dof2GlobalDof0 = &ovlp_prob->Dof2GlobalDof[ip0]; +// Array * Dof2GlobalDof1 = &ovlp_prob->Dof2GlobalDof[ip1]; +// Psi1.SetSize(Dof2GlobalDof1->Size()); Psi1=0.0; +// Vector r(2*bf->FESpace()->GetTrueVSize()); +// r = 0.0; +// r.SetSubVector(*Dof2GlobalDof0,Psi0); + +// // if (direction[0] == 1 && direction[1] == 1) +// // { +// // Vector zloc(Psi1.Size()); +// // r.GetSubVector(*Dof2GlobalDof1,zloc); +// // // extend +// // Vector psi_ext(PmlMat[ip1]->Height()); +// // Vector zloc_ext(PmlMat[ip1]->Height()); +// // Array * Dof2PmlDof1 = &ovlp_prob->Dof2PmlDof[ip1]; +// // zloc_ext.SetSubVector(*Dof2PmlDof1,zloc); +// // PmlMat[ip1]->Mult(zloc_ext,psi_ext); +// // psi_ext *=-1.0; +// // psi_ext.GetSubVector(*Dof2PmlDof1,Psi1); +// // // } +// // else +// // { +// r.GetSubVector(*Dof2GlobalDof1,Psi1); +// // } +// return ip1; +// } + +// void DiagST::ConstructDirectionsMap() +// { +// // total of 8 possible directions of transfer (2D) +// // form left ( 1 , 0) +// // form left-above ( 1 , -1) +// // form left-below ( 1 , 1) +// // form right (-1 , 0) +// // form right-below (-1 , 1) +// // form right-above (-1 , -1) +// // form above ( 0 , -1) +// // form below ( 0 , 1) +// ntransf_directions = pow(3,dim); + +// dirx.SetSize(ntransf_directions); +// diry.SetSize(ntransf_directions); +// int n=3; +// Array ijk(dim); +// if (dim==2) +// { +// for (int i=-1; i<=1; i++) // directions x +// { +// for (int j=-1; j<=1; j++) // directions y +// { +// ijk[0]=i; +// ijk[1]=j; +// int k=GetDirectionId(ijk); +// dirx[k]=i; +// diry[k]=j; +// } +// } +// } +// else if (dim==3) +// { +// dirz.SetSize(ntransf_directions); +// for (int i=-1; i<=1; i++) // directions x +// { +// for (int j=-1; j<=1; j++) // directions y +// { +// for (int k=-1; k<=1; k++) // directions zß +// { +// ijk[0]=i; +// ijk[1]=j; +// ijk[2]=k; +// int l=GetDirectionId(ijk); +// dirx[l]=i; +// diry[l]=j; +// dirz[l]=k; +// } +// } +// } +// } + +// // cout << "dirx = " << endl; +// // dirx.Print(cout,ntransf_directions); +// // cout << "diry = " << endl; +// // diry.Print(cout,ntransf_directions); + +// // if (dim==2) +// // { +// // for (int id=0; id<9; id++) +// // { +// // GetDirectionijk(id,ijk); +// // // cout << "for id = " << id << ": (" < & ijk) const +// { +// int d = ijk.Size(); +// int n=3; +// if (d==2) +// { +// return (ijk[0]+1)*n+(ijk[1]+1); +// } +// else +// { +// return (ijk[0]+1)*n*n+(ijk[1]+1)*n+ijk[2]+1; +// } +// } + +// void DiagST::GetDirectionijk(int id, Array & ijk) const +// { +// int d = ijk.Size(); +// int n=3; +// if (d==2) +// { +// ijk[0]=id/n - 1; +// ijk[1]=id%n - 1; +// } +// else +// { +// ijk[0]=id/(n*n)-1; +// ijk[1]=(id-(ijk[0]+1)*n*n)/n - 1; +// ijk[2]=(id-(ijk[0]+1)*n*n)%n - 1; +// } +// // cout << "ijk = " ; ijk.Print(); +// } + + + + +// void DiagST::TransferSources(int sweep, int ip0, Vector & sol_ext) const +// { +// // Find all neighbors of patch ip +// int nx = nxyz[0]; +// int ny = nxyz[1]; +// int i0, j0, k0; +// Getijk(ip0, i0,j0,k0); +// // cout << "Transfer to : " << endl; +// // loop through possible directions +// for (int i=-1; i<2; i++) +// { +// int i1 = i0 + i; +// if (i1 <0 || i1>=nx) continue; +// for (int j=-1; j<2; j++) +// { +// if (i==0 && j==0) continue; +// int j1 = j0 + j; +// if (j1 <0 || j1>=ny) continue; +// // cout << "(" << i1 << "," << j1 <<"), "; +// // Find ip 1 +// Array ij1(2); ij1[0] = i1; ij1[1]=j1; +// int ip1 = GetPatchId(ij1); +// // cout << "ip1 = " << ip1; +// // cout << " in the direction of (" << i <<", " < directions(2); +// directions[0] = i; +// directions[1] = j; +// Vector cfsol_ext; +// Vector res_ext(sol_ext.Size()); +// GetCutOffSolution(sol_ext,cfsol_ext,ip0,directions,true); +// // sol_ext = cfsol_ext; +// // Calculate source to be transfered + + + + + +// PmlMat[ip0]->Mult(cfsol_ext, res_ext); + +// //--------------------------------------- +// // FiniteElementSpace * fes = ovlp_prob->PmlFespaces[ip0]; +// // Mesh * mesh = fes->GetMesh(); +// // GridFunction gf(fes); +// // double * data = res_ext.GetData(); +// // // gf.SetData(&data[fespace->GetTrueVSize()]); +// // gf.SetData(data); +// // char vishost[] = "localhost"; +// // int visport = 19916; +// // socketstream pmlsock(vishost, visport); + +// // string keys; +// // keys = "keys mrRljc\n"; +// // pmlsock << "solution\n" << *mesh << gf << keys << "valuerange -0.1 0.1 \n" << flush; +// // // sol_sock << "solution\n" << *mesh << gf << keys << flush; +// // cin.get(); +// //--------------------------------------- + + + +// res_ext*= -1.0; +// Array *Dof2PmlDof = &ovlp_prob->Dof2PmlDof[ip0]; +// Vector res_local(Dof2PmlDof->Size()); res_local = 0.0; +// res_ext.GetSubVector(*Dof2PmlDof,res_local); + + +// // Vector sol_local(Dof2PmlDof->Size()); sol_local = 0.0; +// // cfsol_ext.GetSubVector(*Dof2PmlDof,sol_local); +// // Vector znew(A->Height()); znew = 0.0; +// // Vector rnew(A->Height()); rnew = 0.0; +// // Array *Dof2GlobalDof = &ovlp_prob->Dof2GlobalDof[ip0]; +// // znew.SetSubVector(*Dof2GlobalDof,sol_local); +// // A->Mult(znew,rnew); rnew *=-1.0; +// // rnew.GetSubVector(*Dof2GlobalDof, res_local); + +// //----------------------------- +// // pass to ip1 and calculate residual there +// // Vector sol_local(Dof2PmlDof->Size()); sol_local = 0.0; +// // cfsol_ext.GetSubVector(*Dof2PmlDof,sol_local); +// //----------------------------- + + + + + +// // Find the minumum sweep number that to transfer the source that +// // satisfies the two rules +// for (int l=sweep; lSize()==raux.Size(), +// "Transfer Sources: inconsistent size"); +// *f_transf[ip1][l]+=raux; + + +// // FiniteElementSpace * fes1 = ovlp_prob->fespaces[ip1]; +// // Mesh * mesh1 = fes1->GetMesh(); +// // GridFunction gf1(fes1); +// // double * data1 = raux.GetData(); +// // gf1.SetData(data1); +// // socketstream sock(vishost, visport); +// // sock << "solution\n" << *mesh1 << gf1 << keys << "valuerange -0.1 0.1 \n" << flush; +// // // sol_sock << "solution\n" << *mesh << gf << keys << flush; +// // cin.get(); +// break; +// } + +// } +// } +// } + + diff --git a/examples/maxwell-solver/DST/DiagSTnew.hpp b/examples/maxwell-solver/DST/DiagSTnew.hpp new file mode 100644 index 0000000000..f8e5d6bad3 --- /dev/null +++ b/examples/maxwell-solver/DST/DiagSTnew.hpp @@ -0,0 +1,120 @@ +// #pragma once +// #include "Utilities.hpp" +// #include "PML.hpp" +// using namespace std; +// using namespace mfem; + +// class DiagST : public Solver// +// { +// private: +// int nrpatch; +// int dim; +// SesquilinearForm *bf=nullptr; +// MeshPartition * povlp=nullptr; +// double omega = 0.5; +// Coefficient * ws; +// int nrlayers; +// int ovlpnrlayers; +// int nxyz[3]; +// const Operator * A=nullptr; +// Vector B; +// DofMap * ovlp_prob = nullptr; +// Array PmlMat; +// Array PmlMatInv; +// Array2D Pmllength; +// Array3D subdomains; +// mutable Array f_orig; +// int ntransf_directions; +// int nsweeps; +// Array2D sweeps; +// Array dirx; +// Array diry; +// Array dirz; +// mutable Array> f_transf; +// Array> usol; + +// SparseMatrix * GetPmlSystemMatrix(int ip); +// void PlotSolution(Vector & sol, socketstream & sol_sock, int ip) const; +// // void GetCutOffSolution(const Vector & sol, Vector & cfsol, +// // int ip, Array directions, bool local=false) const; +// void GetCutOffSolution(const Vector & sol, Vector & cfsol, +// int ip, Array directions, int ovlpnlayers, bool local=false) const; +// void GetChiRes(const Vector & res, Vector & cfres, +// int ip, Array directions, int nlayers) const; +// void TransferSources(int sweep, int ip, Vector & sol_ext) const; +// int GetDirectionId(const Array & ijk) const; +// void GetDirectionijk(int id, Array & ijk) const; +// void ConstructDirectionsMap(); +// int GetPatchId(const Array & ijk) const; +// void Getijk(int ip, int & i, int & j, int & k ) const; +// int SourceTransfer(const Vector & Psi0, Array direction, int ip, Vector & Psi1) const; +// public: +// DiagST(SesquilinearForm * bf_, Array2D & Pmllength_, +// double omega_, Coefficient * ws_, int nrlayers_); +// void SetLoadVector(Vector load) { B = load;} +// virtual void SetOperator(const Operator &op) {A = &op;} +// virtual void Mult(const Vector &r, Vector &z) const; +// virtual ~DiagST(); +// }; + +#pragma once +#include "Utilities.hpp" +#include "PML.hpp" +using namespace std; +using namespace mfem; + +class DiagST : public Solver// +{ +private: + int nrpatch; + int dim; + SesquilinearForm *bf=nullptr; + MeshPartition * povlp=nullptr; + double omega = 0.5; + Coefficient * ws; + int nrlayers; + int ovlpnrlayers; + int nxyz[3]; + const Operator * A=nullptr; + Vector B; + DofMap * ovlp_prob = nullptr; + Array PmlMat; + Array PmlMatInv; + Array2D Pmllength; + Array3D subdomains; + mutable Array f_orig; + int ntransf_directions; + int nsweeps; + Array2D sweeps; + Array dirx; + Array diry; + Array dirz; + mutable Array> f_transf; + Array> usol; + + SparseMatrix * GetPmlSystemMatrix(int ip); + void PlotSolution(Vector & sol, socketstream & sol_sock, int ip, bool localdomain = false, bool pmldomain = false) const; + // void GetCutOffSolution(const Vector & sol, Vector & cfsol, + // int ip, Array directions, bool local=false) const; + void GetCutOffSolution(const Vector & sol, Vector & cfsol, + int ip, Array directions, int ovlpnlayers, bool local=false) const; + void GetChiRes(const Vector & res, Vector & cfres, + int ip, Array directions, int nlayers) const; + void TransferSources(int sweep, int ip, Vector & sol_ext) const; + int GetDirectionId(const Array & ijk) const; + void GetDirectionijk(int id, Array & ijk) const; + void ConstructDirectionsMap(); + int GetPatchId(const Array & ijk) const; + void Getijk(int ip, int & i, int & j, int & k ) const; + int SourceTransfer(const Vector & Psi0, Array direction, int ip, Vector & Psi1) const; +public: + DiagST(SesquilinearForm * bf_, Array2D & Pmllength_, + double omega_, Coefficient * ws_, int nrlayers_); + void SetLoadVector(Vector load) { B = load;} + virtual void SetOperator(const Operator &op) {A = &op;} + virtual void Mult(const Vector &r, Vector &z) const; + virtual ~DiagST(); +}; + + + diff --git a/examples/maxwell-solver/DST/Utilities.cpp b/examples/maxwell-solver/DST/Utilities.cpp index c0a5b35283..ff5a41d1ce 100644 --- a/examples/maxwell-solver/DST/Utilities.cpp +++ b/examples/maxwell-solver/DST/Utilities.cpp @@ -21,7 +21,7 @@ double CutOffFncn(const Vector &x, const Vector & pmin, const Vector & pmax, con for (int i = 0; i pmax(i) || x(i) < pmin(i)) + if( x(i) >= pmax(i) || x(i) <= pmin(i)) { val = 0.0; } @@ -72,7 +72,7 @@ double ChiFncn(const Vector &x, const Vector & pmin, const Vector & pmax, const for (int i = 0; i pmax(i) || x(i) < pmin(i)) + if( x(i) >= pmax(i) || x(i) <= pmin(i)) { val = 0.0; } diff --git a/examples/maxwell-solver/DST/helmholtz.cpp b/examples/maxwell-solver/DST/helmholtz.cpp index bfbbdeb380..fdfbf91ff2 100644 --- a/examples/maxwell-solver/DST/helmholtz.cpp +++ b/examples/maxwell-solver/DST/helmholtz.cpp @@ -231,7 +231,7 @@ int main(int argc, char *argv[]) - int n= 50; + int n= 1; X = 0.0; Vector z(X.Size()); z = 0.0; Vector r(B); @@ -326,10 +326,10 @@ double f_exact_Re(const Vector &x) double x2 = length/2.0; // x0 = 0.59; // x0 = 0.19; - x0 = 0.15; + x0 = 0.2; // x1 = 0.768; // x1 = 0.168; - x1 = 0.15; + x1 = 0.5; double alpha,beta; // double n = 5.0*omega/M_PI; double n = 4.0*omega/M_PI; @@ -343,8 +343,8 @@ double f_exact_Re(const Vector &x) alpha = -pow(n,2) * beta; f_re = coeff*exp(alpha); - x0 = 0.7; - x1 = 0.7; + x0 = 0.85; + x1 = 0.85; beta = pow(x0-x(0),2) + pow(x1-x(1),2); if (dim == 3) { beta += pow(x2-x(2),2); } alpha = -pow(n,2) * beta;