MFEM  v3.3.2
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
volta_solver.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
2 // the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
3 // reserved. See file COPYRIGHT for details.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability see http://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the GNU Lesser General Public License (as published by the Free
10 // Software Foundation) version 2.1 dated February 1999.
11 
12 #include "volta_solver.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 using namespace std;
17 namespace mfem
18 {
19 using namespace miniapps;
20 
21 namespace electromagnetics
22 {
23 
24 VoltaSolver::VoltaSolver(ParMesh & pmesh, int order,
25  Array<int> & dbcs, Vector & dbcv,
26  Array<int> & nbcs, Vector & nbcv,
27  Coefficient & epsCoef,
28  double (*phi_bc )(const Vector&),
29  double (*rho_src)(const Vector&),
30  void (*p_src )(const Vector&, Vector&))
31  : myid_(0),
32  num_procs_(1),
33  order_(order),
34  pmesh_(&pmesh),
35  dbcs_(&dbcs),
36  dbcv_(&dbcv),
37  nbcs_(&nbcs),
38  nbcv_(&nbcv),
39  visit_dc_(NULL),
40  H1FESpace_(NULL),
41  HCurlFESpace_(NULL),
42  HDivFESpace_(NULL),
43  divEpsGrad_(NULL),
44  h1Mass_(NULL),
45  h1SurfMass_(NULL),
46  hDivMass_(NULL),
47  hCurlHDivEps_(NULL),
48  hCurlHDiv_(NULL),
49  weakDiv_(NULL),
50  grad_(NULL),
51  phi_(NULL),
52  rho_(NULL),
53  rhod_(NULL),
54  sigma_(NULL),
55  e_(NULL),
56  d_(NULL),
57  p_(NULL),
58  epsCoef_(&epsCoef),
59  phiBCCoef_(NULL),
60  rhoCoef_(NULL),
61  pCoef_(NULL),
62  phi_bc_(phi_bc),
63  rho_src_(rho_src),
64  p_src_(p_src)
65 {
66  // Initialize MPI variables
67  MPI_Comm_size(pmesh_->GetComm(), &num_procs_);
68  MPI_Comm_rank(pmesh_->GetComm(), &myid_);
69 
70  // Define compatible parallel finite element spaces on the parallel
71  // mesh. Here we use arbitrary order H1, Nedelec, and Raviart-Thomas finite
72  // elements.
73  H1FESpace_ = new H1_ParFESpace(pmesh_,order,pmesh_->Dimension());
74  HCurlFESpace_ = new ND_ParFESpace(pmesh_,order,pmesh_->Dimension());
75  HDivFESpace_ = new RT_ParFESpace(pmesh_,order,pmesh_->Dimension());
76 
77  // Select surface attributes for Dirichlet BCs
78  ess_bdr_.SetSize(pmesh.bdr_attributes.Max());
79  ess_bdr_ = 0; // Deselect all outer surfaces
80  for (int i=0; i<dbcs_->Size(); i++)
81  {
82  ess_bdr_[(*dbcs_)[i]-1] = 1;
83  }
84 
85  // Setup various coefficients
86 
87  // Potential on outer surface
88  if ( phi_bc_ != NULL )
89  {
90  phiBCCoef_ = new FunctionCoefficient(*phi_bc_);
91  }
92 
93  // Volume Charge Density
94  if ( rho_src_ != NULL )
95  {
96  rhoCoef_ = new FunctionCoefficient(rho_src_);
97  }
98 
99  // Polarization
100  if ( p_src_ != NULL )
101  {
102  pCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
103  p_src_);
104  }
105 
106  // Bilinear Forms
107  divEpsGrad_ = new ParBilinearForm(H1FESpace_);
108  divEpsGrad_->AddDomainIntegrator(new DiffusionIntegrator(*epsCoef_));
109 
110  hDivMass_ = new ParBilinearForm(HDivFESpace_);
112 
113  hCurlHDivEps_ = new ParMixedBilinearForm(HCurlFESpace_,HDivFESpace_);
114  hCurlHDivEps_->AddDomainIntegrator(new VectorFEMassIntegrator(*epsCoef_));
115 
116  // Discrete Grad operator
117  grad_ = new ParDiscreteGradOperator(H1FESpace_, HCurlFESpace_);
118 
119  // Build grid functions
120  phi_ = new ParGridFunction(H1FESpace_);
121  rhod_ = new ParGridFunction(H1FESpace_);
122  d_ = new ParGridFunction(HDivFESpace_);
123  e_ = new ParGridFunction(HCurlFESpace_);
124 
125  if ( rho_src_ )
126  {
127  rho_ = new ParGridFunction(H1FESpace_);
128 
129  h1Mass_ = new ParBilinearForm(H1FESpace_);
130  h1Mass_->AddDomainIntegrator(new MassIntegrator);
131  }
132 
133  if ( p_src_ )
134  {
135  p_ = new ParGridFunction(HCurlFESpace_);
136 
137  hCurlHDiv_ = new ParMixedBilinearForm(HCurlFESpace_, HDivFESpace_);
139 
140  weakDiv_ = new ParMixedBilinearForm(HCurlFESpace_, H1FESpace_);
142 
143  }
144 
145  if ( nbcs_->Size() > 0 )
146  {
147  sigma_ = new ParGridFunction(H1FESpace_);
148 
149  h1SurfMass_ = new ParBilinearForm(H1FESpace_);
150  h1SurfMass_->AddBoundaryIntegrator(new MassIntegrator);
151  }
152 }
153 
155 {
156  delete phiBCCoef_;
157  delete rhoCoef_;
158  delete pCoef_;
159 
160  delete phi_;
161  delete rho_;
162  delete rhod_;
163  delete sigma_;
164  delete d_;
165  delete e_;
166  delete p_;
167 
168  delete grad_;
169 
170  delete divEpsGrad_;
171  delete h1Mass_;
172  delete h1SurfMass_;
173  delete hDivMass_;
174  delete hCurlHDivEps_;
175  delete hCurlHDiv_;
176  delete weakDiv_;
177 
178  delete H1FESpace_;
179  delete HCurlFESpace_;
180  delete HDivFESpace_;
181 
182  map<string,socketstream*>::iterator mit;
183  for (mit=socks_.begin(); mit!=socks_.end(); mit++)
184  {
185  delete mit->second;
186  }
187 }
188 
189 HYPRE_Int
191 {
192  return H1FESpace_->GlobalTrueVSize();
193 }
194 
195 void
197 {
198  HYPRE_Int size_h1 = H1FESpace_->GlobalTrueVSize();
199  HYPRE_Int size_nd = HCurlFESpace_->GlobalTrueVSize();
200  HYPRE_Int size_rt = HDivFESpace_->GlobalTrueVSize();
201  if (myid_ == 0)
202  {
203  cout << "Number of H1 unknowns: " << size_h1 << endl;
204  cout << "Number of H(Curl) unknowns: " << size_nd << endl;
205  cout << "Number of H(Div) unknowns: " << size_rt << endl;
206  }
207 }
208 
210 {
211  if (myid_ == 0) { cout << "Assembling ... " << flush; }
212 
213  divEpsGrad_->Assemble();
214  divEpsGrad_->Finalize();
215 
216  hDivMass_->Assemble();
217  hDivMass_->Finalize();
218 
219  hCurlHDivEps_->Assemble();
220  hCurlHDivEps_->Finalize();
221 
222  grad_->Assemble();
223  grad_->Finalize();
224 
225  if ( h1Mass_ )
226  {
227  h1Mass_->Assemble();
228  h1Mass_->Finalize();
229  }
230  if ( h1SurfMass_ )
231  {
232  h1SurfMass_->Assemble();
233  h1SurfMass_->Finalize();
234  }
235  if ( hCurlHDiv_ )
236  {
237  hCurlHDiv_->Assemble();
238  hCurlHDiv_->Finalize();
239  }
240  if ( weakDiv_ )
241  {
242  weakDiv_->Assemble();
243  weakDiv_->Finalize();
244  }
245 
246  if (myid_ == 0) { cout << "done." << endl << flush; }
247 }
248 
249 void
251 {
252  if (myid_ == 0) { cout << "Updating ..." << endl; }
253 
254  // Inform the spaces that the mesh has changed
255  // Note: we don't need to interpolate any GridFunctions on the new mesh
256  // so we pass 'false' to skip creation of any transformation matrices.
257  H1FESpace_->Update(false);
258  HCurlFESpace_->Update(false);
259  HDivFESpace_->Update(false);
260 
261  // Inform the grid functions that the space has changed.
262  phi_->Update();
263  rhod_->Update();
264  d_->Update();
265  e_->Update();
266  if ( rho_ ) { rho_->Update(); }
267  if ( sigma_ ) { sigma_->Update(); }
268  if ( p_ ) { p_->Update(); }
269 
270  // Inform the bilinear forms that the space has changed.
271  divEpsGrad_->Update();
272  hDivMass_->Update();
273  hCurlHDivEps_->Update();
274 
275  if ( h1Mass_ ) { h1Mass_->Update(); }
276  if ( h1SurfMass_ ) { h1SurfMass_->Update(); }
277  if ( hCurlHDiv_ ) { hCurlHDiv_->Update(); }
278  if ( weakDiv_ ) { weakDiv_->Update(); }
279 
280  // Inform the other objects that the space has changed.
281  grad_->Update();
282 }
283 
284 void
286 {
287  if (myid_ == 0) { cout << "Running solver ... " << endl; }
288 
289  // Initialize the electric potential with its boundary conditions
290  *phi_ = 0.0;
291 
292  // Initialize the charge density dual vector (rhs) to zero
293  *rhod_ = 0.0;
294 
295  if ( dbcs_->Size() > 0 )
296  {
297  if ( phiBCCoef_ )
298  {
299  // Apply gradient boundary condition
300  phi_->ProjectBdrCoefficient(*phiBCCoef_, ess_bdr_);
301  }
302  else
303  {
304  // Apply piecewise constant boundary condition
305  Array<int> dbc_bdr_attr(pmesh_->bdr_attributes.Max());
306  for (int i=0; i<dbcs_->Size(); i++)
307  {
308  ConstantCoefficient voltage((*dbcv_)[i]);
309  dbc_bdr_attr = 0;
310  dbc_bdr_attr[(*dbcs_)[i]-1] = 1;
311  phi_->ProjectBdrCoefficient(voltage, dbc_bdr_attr);
312  }
313  }
314  }
315 
316  // Initialize the volumetric charge density
317  if ( rho_ )
318  {
319  rho_->ProjectCoefficient(*rhoCoef_);
320  h1Mass_->AddMult(*rho_, *rhod_);
321  }
322 
323  // Initialize the Polarization
324  if ( p_ )
325  {
326  p_->ProjectCoefficient(*pCoef_);
327  weakDiv_->AddMult(*p_, *rhod_);
328  }
329 
330  // Initialize the surface charge density
331  if ( sigma_ )
332  {
333  *sigma_ = 0.0;
334 
335  Array<int> nbc_bdr_attr(pmesh_->bdr_attributes.Max());
336  for (int i=0; i<nbcs_->Size(); i++)
337  {
338  ConstantCoefficient sigma_coef((*nbcv_)[i]);
339  nbc_bdr_attr = 0;
340  nbc_bdr_attr[(*nbcs_)[i]-1] = 1;
341  sigma_->ProjectBdrCoefficient(sigma_coef, nbc_bdr_attr);
342  }
343  h1SurfMass_->AddMult(*sigma_, *rhod_);
344  }
345 
346  // Determine the essential BC degrees of freedom
347  if ( dbcs_->Size() > 0 )
348  {
349  // From user supplied boundary attributes
350  H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
351  }
352  else
353  {
354  // Use the first DoF on processor zero by default
355  if ( myid_ == 0 )
356  {
357  ess_bdr_tdofs_.SetSize(1);
358  ess_bdr_tdofs_[0] = 0;
359  }
360  }
361 
362  // Apply essential BC and form linear system
363  HypreParMatrix DivEpsGrad;
364  HypreParVector Phi(H1FESpace_);
365  HypreParVector RHS(H1FESpace_);
366 
367  divEpsGrad_->FormLinearSystem(ess_bdr_tdofs_, *phi_, *rhod_, DivEpsGrad,
368  Phi, RHS);
369 
370  // Define and apply a parallel PCG solver for AX=B with the AMG
371  // preconditioner from hypre.
372  HypreBoomerAMG amg(DivEpsGrad);
373  HyprePCG pcg(DivEpsGrad);
374  pcg.SetTol(1e-12);
375  pcg.SetMaxIter(500);
376  pcg.SetPrintLevel(2);
377  pcg.SetPreconditioner(amg);
378  pcg.Mult(RHS, Phi);
379 
380  // Extract the parallel grid function corresponding to the finite
381  // element approximation Phi. This is the local solution on each
382  // processor.
383  divEpsGrad_->RecoverFEMSolution(Phi, *rhod_, *phi_);
384 
385  // Compute the negative Gradient of the solution vector. This is
386  // the magnetic field corresponding to the scalar potential
387  // represented by phi.
388  grad_->Mult(*phi_, *e_); *e_ *= -1.0;
389 
390  // Compute electric displacement (D) from E and P (if present)
391  if (myid_ == 0) { cout << "Computing D ..." << flush; }
392 
393  ParGridFunction ed(HDivFESpace_);
394  hCurlHDivEps_->Mult(*e_, ed);
395  if ( p_ )
396  {
397  hCurlHDiv_->AddMult(*p_, ed, -1.0);
398  }
399 
400  HypreParMatrix MassHDiv;
401  Vector ED, D;
402 
403  Array<int> dbc_dofs_d;
404  hDivMass_->FormLinearSystem(dbc_dofs_d, *d_, ed, MassHDiv, D, ED);
405 
406  HyprePCG pcgM(MassHDiv);
407  pcgM.SetTol(1e-12);
408  pcgM.SetMaxIter(500);
409  pcgM.SetPrintLevel(0);
410  HypreDiagScale diagM;
411  pcgM.SetPreconditioner(diagM);
412  pcgM.Mult(ED, D);
413 
414  hDivMass_->RecoverFEMSolution(D, ed, *d_);
415 
416  if (myid_ == 0) { cout << "done." << flush; }
417 
418  if (myid_ == 0) { cout << "Solver done. " << endl; }
419 }
420 
421 void
423 {
424  if (myid_ == 0) { cout << "Estimating Error ... " << flush; }
425 
426  // Space for the discontinuous (original) flux
427  DiffusionIntegrator flux_integrator(*epsCoef_);
428  L2_FECollection flux_fec(order_, pmesh_->Dimension());
429  // ND_FECollection flux_fec(order_, pmesh_->Dimension());
430  ParFiniteElementSpace flux_fes(pmesh_, &flux_fec, pmesh_->SpaceDimension());
431 
432  // Space for the smoothed (conforming) flux
433  double norm_p = 1;
434  RT_FECollection smooth_flux_fec(order_-1, pmesh_->Dimension());
435  ParFiniteElementSpace smooth_flux_fes(pmesh_, &smooth_flux_fec);
436 
437  L2ZZErrorEstimator(flux_integrator, *phi_,
438  smooth_flux_fes, flux_fes, errors, norm_p);
439 
440  if (myid_ == 0) { cout << "done." << endl; }
441 }
442 
443 void
445 {
446  visit_dc_ = &visit_dc;
447 
448  visit_dc.RegisterField("Phi", phi_);
449  visit_dc.RegisterField("D", d_);
450  visit_dc.RegisterField("E", e_);
451  if ( rho_ ) { visit_dc.RegisterField("Rho", rho_); }
452  if ( p_ ) { visit_dc.RegisterField("P", p_); }
453  if ( sigma_ ) { visit_dc.RegisterField("Sigma", sigma_); }
454 }
455 
456 void
458 {
459  if ( visit_dc_ )
460  {
461  if (myid_ == 0) { cout << "Writing VisIt files ..." << flush; }
462 
463  HYPRE_Int prob_size = this->GetProblemSize();
464  visit_dc_->SetCycle(it);
465  visit_dc_->SetTime(prob_size);
466  visit_dc_->Save();
467 
468  if (myid_ == 0) { cout << " done." << endl; }
469  }
470 }
471 
472 void
474 {
475  if ( myid_ == 0 ) { cout << "Opening GLVis sockets." << endl; }
476 
477  socks_["Phi"] = new socketstream;
478  socks_["Phi"]->precision(8);
479 
480  socks_["D"] = new socketstream;
481  socks_["D"]->precision(8);
482 
483  socks_["E"] = new socketstream;
484  socks_["E"]->precision(8);
485 
486  if ( rho_)
487  {
488  socks_["Rho"] = new socketstream;
489  socks_["Rho"]->precision(8);
490  }
491  if ( p_)
492  {
493  socks_["P"] = new socketstream;
494  socks_["P"]->precision(8);
495  }
496  if ( sigma_)
497  {
498  socks_["Sigma"] = new socketstream;
499  socks_["Sigma"]->precision(8);
500  }
501 }
502 
503 void
505 {
506  if (myid_ == 0) { cout << "Sending data to GLVis ..." << flush; }
507 
508  char vishost[] = "localhost";
509  int visport = 19916;
510 
511  int Wx = 0, Wy = 0; // window position
512  int Ww = 350, Wh = 350; // window size
513  int offx = Ww+10, offy = Wh+45; // window offsets
514 
515  VisualizeField(*socks_["Phi"], vishost, visport,
516  *phi_, "Electric Potential (Phi)", Wx, Wy, Ww, Wh);
517  Wx += offx;
518 
519  VisualizeField(*socks_["D"], vishost, visport,
520  *d_, "Electric Displacement (D)", Wx, Wy, Ww, Wh);
521  Wx += offx;
522 
523  VisualizeField(*socks_["E"], vishost, visport,
524  *e_, "Electric Field (E)", Wx, Wy, Ww, Wh);
525 
526  Wx = 0; Wy += offy; // next line
527 
528  if ( rho_ )
529  {
530  VisualizeField(*socks_["Rho"], vishost, visport,
531  *rho_, "Charge Density (Rho)", Wx, Wy, Ww, Wh);
532  Wx += offx;
533  }
534  if ( p_ )
535  {
536  VisualizeField(*socks_["P"], vishost, visport,
537  *p_, "Electric Polarization (P)", Wx, Wy, Ww, Wh);
538  Wx += offx;
539  }
540  if ( sigma_ )
541  {
542  VisualizeField(*socks_["Sigma"], vishost, visport,
543  *sigma_, "Surface Charge Density (Sigma)", Wx, Wy, Ww, Wh);
544  Wx += offx;
545  }
546  if (myid_ == 0) { cout << " done." << endl; }
547 }
548 
549 } // namespace electromagnetics
550 
551 } // namespace mfem
552 
553 #endif // MFEM_USE_MPI
void SetTol(double tol)
Definition: hypre.cpp:2088
int Size() const
Logical size of the array.
Definition: array.hpp:110
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:529
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
void SetCycle(int c)
Set time cycle (for time-dependent simulations)
Subclass constant coefficient.
Definition: coefficient.hpp:57
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:2225
virtual void Save()
Save the collection and a VisIt root file.
Abstract parallel finite element space.
Definition: pfespace.hpp:31
void RegisterVisItFields(VisItDataCollection &visit_dc)
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:275
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2103
virtual void Update(FiniteElementSpace *nfes=NULL)
void GetErrorEstimates(Vector &errors)
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
The BoomerAMG solver in hypre.
Definition: hypre.hpp:783
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:65
Jacobi preconditioner in hypre.
Definition: hypre.hpp:744
Data collection with VisIt I/O routines.
void VisualizeField(socketstream &sock, const char *vishost, int visport, GridFunction &gf, const char *title, int x, int y, int w, int h, bool vec)
Definition: fem_extras.cpp:59
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:108
void Assemble(int skip_zeros=1)
HYPRE_Int GlobalTrueVSize() const
Definition: pfespace.hpp:178
void Assemble(int skip_zeros=1)
Assemble the local matrix.
int Dimension() const
Definition: mesh.hpp:641
void SetTime(double t)
Set physical time (for time-dependent simulations)
void SetMaxIter(int max_iter)
Definition: hypre.cpp:2093
Arbitrary order H(div)-conforming Raviart-Thomas finite elements.
Definition: fe_coll.hpp:239
int SpaceDimension() const
Definition: mesh.hpp:642
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:73
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition: mesh.hpp:172
MPI_Comm GetComm() const
Definition: pmesh.hpp:117
PCG solver in hypre.
Definition: hypre.hpp:643
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:503
double L2ZZErrorEstimator(BilinearFormIntegrator &flux_integrator, const ParGridFunction &x, ParFiniteElementSpace &smooth_flux_fes, ParFiniteElementSpace &flux_fes, Vector &errors, int norm_p, double solver_tol, int solver_max_it)
Definition: pgridfunc.cpp:591
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
virtual void RegisterField(const std::string &field_name, GridFunction *gf)
Add a grid function to the collection and update the root file.
Class for parallel bilinear form using different test and trial FE spaces.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
void AddBoundaryIntegrator(BilinearFormIntegrator *bfi)
Adds new Boundary Integrator.
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:2108
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator.
Class for parallel bilinear form.
class for C-function coefficient
Vector data type.
Definition: vector.hpp:41
Class for parallel grid function.
Definition: pgridfunc.hpp:32
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:175
virtual void Assemble(int skip_zeros=1)
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:2129
Class for parallel meshes.
Definition: pmesh.hpp:29
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Integrator for (Q u, v) for VectorFiniteElements.
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Definition: gridfunc.hpp:213
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:195