1 // Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 #include "volta_solver.hpp"
13
14 #ifdef MFEM_USE_MPI
15
16 using namespace std;
17 namespace mfem
18 {
19
20 using namespace common;
21
22 namespace electromagnetics
23 {
24
VoltaSolver(ParMesh & pmesh,int order,Array<int> & dbcs,Vector & dbcv,Array<int> & nbcs,Vector & nbcv,Coefficient & epsCoef,double (* phi_bc)(const Vector &),double (* rho_src)(const Vector &),void (* p_src)(const Vector &,Vector &),Vector & point_charges)25 VoltaSolver::VoltaSolver(ParMesh & pmesh, int order,
26 Array<int> & dbcs, Vector & dbcv,
27 Array<int> & nbcs, Vector & nbcv,
28 Coefficient & epsCoef,
29 double (*phi_bc )(const Vector&),
30 double (*rho_src)(const Vector&),
31 void (*p_src )(const Vector&, Vector&),
32 Vector & point_charges)
33 : myid_(0),
34 num_procs_(1),
35 order_(order),
36 pmesh_(&pmesh),
37 dbcs_(&dbcs),
38 dbcv_(&dbcv),
39 nbcs_(&nbcs),
40 nbcv_(&nbcv),
41 visit_dc_(NULL),
42 H1FESpace_(NULL),
43 HCurlFESpace_(NULL),
44 HDivFESpace_(NULL),
45 L2FESpace_(NULL),
46 divEpsGrad_(NULL),
47 h1Mass_(NULL),
48 h1SurfMass_(NULL),
49 hDivMass_(NULL),
50 hCurlHDivEps_(NULL),
51 hCurlHDiv_(NULL),
52 weakDiv_(NULL),
53 rhod_(NULL),
54 l2_vol_int_(NULL),
55 rt_surf_int_(NULL),
56 grad_(NULL),
57 phi_(NULL),
58 rho_src_(NULL),
59 rho_(NULL),
60 sigma_src_(NULL),
61 e_(NULL),
62 d_(NULL),
63 p_src_(NULL),
64 oneCoef_(1.0),
65 epsCoef_(&epsCoef),
66 phiBCCoef_(NULL),
67 rhoCoef_(NULL),
68 pCoef_(NULL),
69 phi_bc_func_(phi_bc),
70 rho_src_func_(rho_src),
71 p_src_func_(p_src),
72 point_charge_params_(point_charges),
73 point_charges_(0)
74 {
75 // Initialize MPI variables
76 MPI_Comm_size(pmesh_->GetComm(), &num_procs_);
77 MPI_Comm_rank(pmesh_->GetComm(), &myid_);
78
79 // Define compatible parallel finite element spaces on the parallel
80 // mesh. Here we use arbitrary order H1, Nedelec, and Raviart-Thomas finite
81 // elements.
82 H1FESpace_ = new H1_ParFESpace(pmesh_,order,pmesh_->Dimension());
83 HCurlFESpace_ = new ND_ParFESpace(pmesh_,order,pmesh_->Dimension());
84 HDivFESpace_ = new RT_ParFESpace(pmesh_,order,pmesh_->Dimension());
85 L2FESpace_ = new L2_ParFESpace(pmesh_,order-1,pmesh_->Dimension());
86
87 // Select surface attributes for Dirichlet BCs
88 AttrToMarker(pmesh.bdr_attributes.Max(), *dbcs_, ess_bdr_);
89
90 // Setup various coefficients
91
92 // Potential on outer surface
93 if ( phi_bc_func_ != NULL )
94 {
95 phiBCCoef_ = new FunctionCoefficient(*phi_bc_func_);
96 }
97
98 // Volume Charge Density
99 if ( rho_src_func_ != NULL )
100 {
101 rhoCoef_ = new FunctionCoefficient(rho_src_func_);
102 }
103
104 // Polarization
105 if ( p_src_func_ != NULL )
106 {
107 pCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
108 p_src_func_);
109 }
110
111 // Bilinear Forms
112 divEpsGrad_ = new ParBilinearForm(H1FESpace_);
113 divEpsGrad_->AddDomainIntegrator(new DiffusionIntegrator(*epsCoef_));
114
115 hDivMass_ = new ParBilinearForm(HDivFESpace_);
116 hDivMass_->AddDomainIntegrator(new VectorFEMassIntegrator);
117
118 hCurlHDivEps_ = new ParMixedBilinearForm(HCurlFESpace_,HDivFESpace_);
119 hCurlHDivEps_->AddDomainIntegrator(new VectorFEMassIntegrator(*epsCoef_));
120
121 rhod_ = new ParLinearForm(H1FESpace_);
122
123 l2_vol_int_ = new ParLinearForm(L2FESpace_);
124 l2_vol_int_->AddDomainIntegrator(new DomainLFIntegrator(oneCoef_));
125
126 rt_surf_int_ = new ParLinearForm(HDivFESpace_);
127 rt_surf_int_->AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator);
128
129 // Discrete derivative operator
130 grad_ = new ParDiscreteGradOperator(H1FESpace_, HCurlFESpace_);
131 div_ = new ParDiscreteDivOperator(HDivFESpace_, L2FESpace_);
132
133 // Build grid functions
134 phi_ = new ParGridFunction(H1FESpace_);
135 d_ = new ParGridFunction(HDivFESpace_);
136 e_ = new ParGridFunction(HCurlFESpace_);
137 rho_ = new ParGridFunction(L2FESpace_);
138
139 if ( point_charge_params_.Size() > 0 )
140 {
141 int dim = pmesh_->Dimension();
142 int npts = point_charge_params_.Size() / (dim + 1);
143 point_charges_.resize(npts);
144
145 Vector cent(dim);
146 for (int i=0; i<npts; i++)
147 {
148 for (int d=0; d<dim; d++)
149 {
150 cent[d] = point_charge_params_[(dim + 1) * i + d];
151 }
152 double s = point_charge_params_[(dim + 1) * i + dim];
153
154 point_charges_[i] = new DeltaCoefficient();
155 point_charges_[i]->SetScale(s);
156 point_charges_[i]->SetDeltaCenter(cent);
157
158 rhod_->AddDomainIntegrator(new DomainLFIntegrator(*point_charges_[i]));
159 }
160 }
161
162 if ( rho_src_func_ )
163 {
164 rho_src_ = new ParGridFunction(H1FESpace_);
165
166 h1Mass_ = new ParBilinearForm(H1FESpace_);
167 h1Mass_->AddDomainIntegrator(new MassIntegrator);
168 }
169
170 if ( p_src_func_ )
171 {
172 p_src_ = new ParGridFunction(HCurlFESpace_);
173
174 hCurlHDiv_ = new ParMixedBilinearForm(HCurlFESpace_, HDivFESpace_);
175 hCurlHDiv_->AddDomainIntegrator(new VectorFEMassIntegrator);
176
177 weakDiv_ = new ParMixedBilinearForm(HCurlFESpace_, H1FESpace_);
178 weakDiv_->AddDomainIntegrator(new VectorFEWeakDivergenceIntegrator);
179 }
180
181 if ( nbcs_->Size() > 0 )
182 {
183 sigma_src_ = new ParGridFunction(H1FESpace_);
184
185 h1SurfMass_ = new ParBilinearForm(H1FESpace_);
186 h1SurfMass_->AddBoundaryIntegrator(new MassIntegrator);
187 }
188 }
189
~VoltaSolver()190 VoltaSolver::~VoltaSolver()
191 {
192 delete phiBCCoef_;
193 delete rhoCoef_;
194 delete pCoef_;
195
196 delete phi_;
197 delete rho_src_;
198 delete rho_;
199 delete rhod_;
200 delete l2_vol_int_;
201 delete rt_surf_int_;
202 delete sigma_src_;
203 delete d_;
204 delete e_;
205 delete p_src_;
206
207 delete grad_;
208 delete div_;
209
210 delete divEpsGrad_;
211 delete h1Mass_;
212 delete h1SurfMass_;
213 delete hDivMass_;
214 delete hCurlHDivEps_;
215 delete hCurlHDiv_;
216 delete weakDiv_;
217
218 delete H1FESpace_;
219 delete HCurlFESpace_;
220 delete HDivFESpace_;
221 delete L2FESpace_;
222
223 for (unsigned int i=0; i<point_charges_.size(); i++)
224 {
225 delete point_charges_[i];
226 }
227
228 map<string,socketstream*>::iterator mit;
229 for (mit=socks_.begin(); mit!=socks_.end(); mit++)
230 {
231 delete mit->second;
232 }
233 }
234
235 HYPRE_BigInt
GetProblemSize()236 VoltaSolver::GetProblemSize()
237 {
238 return H1FESpace_->GlobalTrueVSize();
239 }
240
241 void
PrintSizes()242 VoltaSolver::PrintSizes()
243 {
244 HYPRE_BigInt size_h1 = H1FESpace_->GlobalTrueVSize();
245 HYPRE_BigInt size_nd = HCurlFESpace_->GlobalTrueVSize();
246 HYPRE_BigInt size_rt = HDivFESpace_->GlobalTrueVSize();
247 HYPRE_BigInt size_l2 = L2FESpace_->GlobalTrueVSize();
248 if (myid_ == 0)
249 {
250 cout << "Number of H1 unknowns: " << size_h1 << endl;
251 cout << "Number of H(Curl) unknowns: " << size_nd << endl;
252 cout << "Number of H(Div) unknowns: " << size_rt << endl;
253 cout << "Number of L2 unknowns: " << size_l2 << endl;
254 }
255 }
256
Assemble()257 void VoltaSolver::Assemble()
258 {
259 if (myid_ == 0) { cout << "Assembling ... " << flush; }
260
261 divEpsGrad_->Assemble();
262 divEpsGrad_->Finalize();
263
264 hDivMass_->Assemble();
265 hDivMass_->Finalize();
266
267 hCurlHDivEps_->Assemble();
268 hCurlHDivEps_->Finalize();
269
270 *rhod_ = 0.0;
271 rhod_->Assemble();
272
273 l2_vol_int_->Assemble();
274 rt_surf_int_->Assemble();
275
276 grad_->Assemble();
277 grad_->Finalize();
278
279 div_->Assemble();
280 div_->Finalize();
281
282 if ( h1Mass_ )
283 {
284 h1Mass_->Assemble();
285 h1Mass_->Finalize();
286 }
287 if ( h1SurfMass_ )
288 {
289 h1SurfMass_->Assemble();
290 h1SurfMass_->Finalize();
291 }
292 if ( hCurlHDiv_ )
293 {
294 hCurlHDiv_->Assemble();
295 hCurlHDiv_->Finalize();
296 }
297 if ( weakDiv_ )
298 {
299 weakDiv_->Assemble();
300 weakDiv_->Finalize();
301 }
302
303 if (myid_ == 0) { cout << "done." << endl << flush; }
304 }
305
306 void
Update()307 VoltaSolver::Update()
308 {
309 if (myid_ == 0) { cout << "Updating ..." << endl; }
310
311 // Inform the spaces that the mesh has changed
312 // Note: we don't need to interpolate any GridFunctions on the new mesh
313 // so we pass 'false' to skip creation of any transformation matrices.
314 H1FESpace_->Update(false);
315 HCurlFESpace_->Update(false);
316 HDivFESpace_->Update(false);
317 L2FESpace_->Update(false);
318
319 // Inform the grid functions that the space has changed.
320 phi_->Update();
321 rhod_->Update();
322 l2_vol_int_->Update();
323 rt_surf_int_->Update();
324 d_->Update();
325 e_->Update();
326 rho_->Update();
327 if ( rho_src_ ) { rho_src_->Update(); }
328 if ( sigma_src_ ) { sigma_src_->Update(); }
329 if ( p_src_ ) { p_src_->Update(); }
330
331 // Inform the bilinear forms that the space has changed.
332 divEpsGrad_->Update();
333 hDivMass_->Update();
334 hCurlHDivEps_->Update();
335
336 if ( h1Mass_ ) { h1Mass_->Update(); }
337 if ( h1SurfMass_ ) { h1SurfMass_->Update(); }
338 if ( hCurlHDiv_ ) { hCurlHDiv_->Update(); }
339 if ( weakDiv_ ) { weakDiv_->Update(); }
340
341 // Inform the other objects that the space has changed.
342 grad_->Update();
343 div_->Update();
344 }
345
346 void
Solve()347 VoltaSolver::Solve()
348 {
349 if (myid_ == 0) { cout << "Running solver ... " << endl; }
350
351 // Initialize the electric potential with its boundary conditions
352 *phi_ = 0.0;
353
354 if ( dbcs_->Size() > 0 )
355 {
356 if ( phiBCCoef_ )
357 {
358 // Apply gradient boundary condition
359 phi_->ProjectBdrCoefficient(*phiBCCoef_, ess_bdr_);
360 }
361 else
362 {
363 // Apply piecewise constant boundary condition
364 Array<int> dbc_bdr_attr(pmesh_->bdr_attributes.Max());
365 for (int i=0; i<dbcs_->Size(); i++)
366 {
367 ConstantCoefficient voltage((*dbcv_)[i]);
368 dbc_bdr_attr = 0;
369 if ((*dbcs_)[i] <= dbc_bdr_attr.Size())
370 {
371 dbc_bdr_attr[(*dbcs_)[i]-1] = 1;
372 }
373 phi_->ProjectBdrCoefficient(voltage, dbc_bdr_attr);
374 }
375 }
376 }
377
378 // Initialize the volumetric charge density
379 if ( rho_src_ )
380 {
381 rho_src_->ProjectCoefficient(*rhoCoef_);
382 h1Mass_->AddMult(*rho_src_, *rhod_);
383 }
384
385 // Initialize the Polarization
386 if ( p_src_ )
387 {
388 p_src_->ProjectCoefficient(*pCoef_);
389 weakDiv_->AddMult(*p_src_, *rhod_);
390 }
391
392 // Initialize the surface charge density
393 if ( sigma_src_ )
394 {
395 *sigma_src_ = 0.0;
396
397 Array<int> nbc_bdr_attr(pmesh_->bdr_attributes.Max());
398 for (int i=0; i<nbcs_->Size(); i++)
399 {
400 ConstantCoefficient sigma_coef((*nbcv_)[i]);
401 nbc_bdr_attr = 0;
402 if ((*nbcs_)[i] <= nbc_bdr_attr.Size())
403 {
404 nbc_bdr_attr[(*nbcs_)[i]-1] = 1;
405 }
406 sigma_src_->ProjectBdrCoefficient(sigma_coef, nbc_bdr_attr);
407 }
408 h1SurfMass_->AddMult(*sigma_src_, *rhod_);
409 }
410
411 // Determine the essential BC degrees of freedom
412 if ( dbcs_->Size() > 0 )
413 {
414 // From user supplied boundary attributes
415 H1FESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
416 }
417 else
418 {
419 // Use the first DoF on processor zero by default
420 if ( myid_ == 0 )
421 {
422 ess_bdr_tdofs_.SetSize(1);
423 ess_bdr_tdofs_[0] = 0;
424 }
425 }
426
427 // Apply essential BC and form linear system
428 HypreParMatrix DivEpsGrad;
429 HypreParVector Phi(H1FESpace_);
430 HypreParVector RHS(H1FESpace_);
431
432 divEpsGrad_->FormLinearSystem(ess_bdr_tdofs_, *phi_, *rhod_, DivEpsGrad,
433 Phi, RHS);
434
435 // Define and apply a parallel PCG solver for AX=B with the AMG
436 // preconditioner from hypre.
437 HypreBoomerAMG amg(DivEpsGrad);
438 HyprePCG pcg(DivEpsGrad);
439 pcg.SetTol(1e-12);
440 pcg.SetMaxIter(500);
441 pcg.SetPrintLevel(2);
442 pcg.SetPreconditioner(amg);
443 pcg.Mult(RHS, Phi);
444
445 // Extract the parallel grid function corresponding to the finite
446 // element approximation Phi. This is the local solution on each
447 // processor.
448 divEpsGrad_->RecoverFEMSolution(Phi, *rhod_, *phi_);
449
450 // Compute the negative Gradient of the solution vector. This is
451 // the magnetic field corresponding to the scalar potential
452 // represented by phi.
453 grad_->Mult(*phi_, *e_); *e_ *= -1.0;
454
455 // Compute electric displacement (D) from E and P (if present)
456 if (myid_ == 0) { cout << "Computing D ..." << flush; }
457
458 ParGridFunction ed(HDivFESpace_);
459 hCurlHDivEps_->Mult(*e_, ed);
460 if ( p_src_ )
461 {
462 hCurlHDiv_->AddMult(*p_src_, ed, -1.0);
463 }
464
465 HypreParMatrix MassHDiv;
466 Vector ED, D;
467
468 Array<int> dbc_dofs_d;
469 hDivMass_->FormLinearSystem(dbc_dofs_d, *d_, ed, MassHDiv, D, ED);
470
471 HyprePCG pcgM(MassHDiv);
472 pcgM.SetTol(1e-12);
473 pcgM.SetMaxIter(500);
474 pcgM.SetPrintLevel(0);
475 HypreDiagScale diagM;
476 pcgM.SetPreconditioner(diagM);
477 pcgM.Mult(ED, D);
478
479 hDivMass_->RecoverFEMSolution(D, ed, *d_);
480
481 // Compute charge density from rho = Div(D)
482 div_->Mult(*d_, *rho_);
483
484 if (myid_ == 0) { cout << "done." << flush; }
485
486 {
487 // Compute total charge as volume integral of rho
488 double charge_rho = (*l2_vol_int_)(*rho_);
489
490 // Compute total charge as surface integral of D
491 double charge_D = (*rt_surf_int_)(*d_);
492
493 if (myid_ == 0)
494 {
495 cout << endl << "Total charge: \n"
496 << " Volume integral of charge density: " << charge_rho
497 << "\n Surface integral of dielectric flux: " << charge_D
498 << endl << flush;
499 }
500 }
501
502 if (myid_ == 0) { cout << "Solver done. " << endl; }
503 }
504
505 void
GetErrorEstimates(Vector & errors)506 VoltaSolver::GetErrorEstimates(Vector & errors)
507 {
508 if (myid_ == 0) { cout << "Estimating Error ... " << flush; }
509
510 // Space for the discontinuous (original) flux
511 DiffusionIntegrator flux_integrator(*epsCoef_);
512 L2_FECollection flux_fec(order_, pmesh_->Dimension());
513 // ND_FECollection flux_fec(order_, pmesh_->Dimension());
514 ParFiniteElementSpace flux_fes(pmesh_, &flux_fec, pmesh_->SpaceDimension());
515
516 // Space for the smoothed (conforming) flux
517 double norm_p = 1;
518 RT_FECollection smooth_flux_fec(order_-1, pmesh_->Dimension());
519 ParFiniteElementSpace smooth_flux_fes(pmesh_, &smooth_flux_fec);
520
521 L2ZZErrorEstimator(flux_integrator, *phi_,
522 smooth_flux_fes, flux_fes, errors, norm_p);
523
524 if (myid_ == 0) { cout << "done." << endl; }
525 }
526
527 void
RegisterVisItFields(VisItDataCollection & visit_dc)528 VoltaSolver::RegisterVisItFields(VisItDataCollection & visit_dc)
529 {
530 visit_dc_ = &visit_dc;
531
532 visit_dc.RegisterField("Phi", phi_);
533 visit_dc.RegisterField("D", d_);
534 visit_dc.RegisterField("E", e_);
535 visit_dc.RegisterField("Rho", rho_);
536 if ( rho_src_ ) { visit_dc.RegisterField("Rho Source", rho_src_); }
537 if ( p_src_ ) { visit_dc.RegisterField("P Source", p_src_); }
538 if ( sigma_src_ ) { visit_dc.RegisterField("Sigma Source", sigma_src_); }
539 }
540
541 void
WriteVisItFields(int it)542 VoltaSolver::WriteVisItFields(int it)
543 {
544 if ( visit_dc_ )
545 {
546 if (myid_ == 0) { cout << "Writing VisIt files ..." << flush; }
547
548 HYPRE_BigInt prob_size = this->GetProblemSize();
549 visit_dc_->SetCycle(it);
550 visit_dc_->SetTime(prob_size);
551 visit_dc_->Save();
552
553 if (myid_ == 0) { cout << " done." << endl; }
554 }
555 }
556
557 void
InitializeGLVis()558 VoltaSolver::InitializeGLVis()
559 {
560 if ( myid_ == 0 ) { cout << "Opening GLVis sockets." << endl; }
561
562 socks_["Phi"] = new socketstream;
563 socks_["Phi"]->precision(8);
564
565 socks_["D"] = new socketstream;
566 socks_["D"]->precision(8);
567
568 socks_["E"] = new socketstream;
569 socks_["E"]->precision(8);
570
571 socks_["Rho"] = new socketstream;
572 socks_["Rho"]->precision(8);
573
574 if ( rho_src_ )
575 {
576 socks_["RhoSrc"] = new socketstream;
577 socks_["RhoSrc"]->precision(8);
578 }
579 if ( p_src_ )
580 {
581 socks_["PSrc"] = new socketstream;
582 socks_["PSrc"]->precision(8);
583 }
584 if ( sigma_src_ )
585 {
586 socks_["SigmaSrc"] = new socketstream;
587 socks_["SigmaSrc"]->precision(8);
588 }
589 }
590
591 void
DisplayToGLVis()592 VoltaSolver::DisplayToGLVis()
593 {
594 if (myid_ == 0) { cout << "Sending data to GLVis ..." << flush; }
595
596 char vishost[] = "localhost";
597 int visport = 19916;
598
599 int Wx = 0, Wy = 0; // window position
600 int Ww = 350, Wh = 350; // window size
601 int offx = Ww+10, offy = Wh+45; // window offsets
602
603 VisualizeField(*socks_["Phi"], vishost, visport,
604 *phi_, "Electric Potential (Phi)", Wx, Wy, Ww, Wh);
605 Wx += offx;
606
607 VisualizeField(*socks_["E"], vishost, visport,
608 *e_, "Electric Field (E)", Wx, Wy, Ww, Wh);
609 Wx += offx;
610
611 VisualizeField(*socks_["D"], vishost, visport,
612 *d_, "Electric Displacement (D)", Wx, Wy, Ww, Wh);
613 Wx += offx;
614
615 VisualizeField(*socks_["Rho"], vishost, visport,
616 *rho_, "Charge Density", Wx, Wy, Ww, Wh);
617 Wx = 0; Wy += offy; // next line
618
619 if ( rho_src_ )
620 {
621 VisualizeField(*socks_["RhoSrc"], vishost, visport,
622 *rho_src_, "Charge Density Source (Rho)", Wx, Wy, Ww, Wh);
623 Wx += offx;
624 }
625 if ( p_src_ )
626 {
627 VisualizeField(*socks_["PSrc"], vishost, visport,
628 *p_src_, "Electric Polarization Source (P)",
629 Wx, Wy, Ww, Wh);
630 Wx += offx;
631 }
632 if ( sigma_src_ )
633 {
634 VisualizeField(*socks_["SigmaSrc"], vishost, visport,
635 *sigma_src_, "Surface Charge Density Source (Sigma)",
636 Wx, Wy, Ww, Wh);
637 // Wx += offx; // not used
638 }
639 if (myid_ == 0) { cout << " done." << endl; }
640 }
641
642 } // namespace electromagnetics
643
644 } // namespace mfem
645
646 #endif // MFEM_USE_MPI
647