MFEM  v4.5.1
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
joule_solver.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2022, 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 "joule_solver.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 using namespace std;
17 
18 namespace mfem
19 {
20 
21 using namespace common;
22 
23 namespace electromagnetics
24 {
25 
26 MagneticDiffusionEOperator::MagneticDiffusionEOperator(
27  int stateVectorLen,
28  ParFiniteElementSpace &L2FES,
29  ParFiniteElementSpace &HCurlFES,
30  ParFiniteElementSpace &HDivFES,
31  ParFiniteElementSpace &HGradFES,
32  Array<int> &ess_bdr_arg,
33  Array<int> &thermal_ess_bdr_arg,
34  Array<int> &poisson_ess_bdr_arg,
35  double mu_coef,
36  std::map<int, double> sigmaAttMap,
37  std::map<int, double> TcapacityAttMap,
38  std::map<int, double> InvTcapAttMap,
39  std::map<int, double> InvTcondAttMap)
40 
41  : TimeDependentOperator(stateVectorLen, 0.0),
42  L2FESpace(L2FES), HCurlFESpace(HCurlFES), HDivFESpace(HDivFES),
43  HGradFESpace(HGradFES),
44  a0(NULL), a1(NULL), a2(NULL), m1(NULL), m2(NULL), m3(NULL),
45  s1(NULL), s2(NULL), grad(NULL), curl(NULL), weakDiv(NULL), weakDivC(NULL),
46  weakCurl(NULL),
47  A0(NULL), A1(NULL), A2(NULL), M1(NULL), M2(NULL), M3(NULL),
48  X0(NULL), X1(NULL), X2(NULL), B0(NULL), B1(NULL), B2(NULL), B3(NULL),
49  v0(NULL), v1(NULL), v2(NULL),
50  amg_a0(NULL), pcg_a0(NULL), ads_a2(NULL), pcg_a2(NULL), ams_a1(NULL),
51  pcg_a1(NULL), dsp_m3(NULL),pcg_m3(NULL),
52  dsp_m1(NULL), pcg_m1(NULL), dsp_m2(NULL), pcg_m2(NULL),
53  mu(mu_coef), dt_A1(-1.0), dt_A2(-1.0)
54 {
55  ess_bdr.SetSize(ess_bdr_arg.Size());
56  for (int i=0; i<ess_bdr_arg.Size(); i++)
57  {
58  ess_bdr[i] = ess_bdr_arg[i];
59  }
60  thermal_ess_bdr.SetSize(thermal_ess_bdr_arg.Size());
61  for (int i=0; i<thermal_ess_bdr_arg.Size(); i++)
62  {
63  thermal_ess_bdr[i] = thermal_ess_bdr_arg[i];
64  }
65  poisson_ess_bdr.SetSize(poisson_ess_bdr_arg.Size());
66  for (int i=0; i<poisson_ess_bdr_arg.Size(); i++)
67  {
68  poisson_ess_bdr[i] = poisson_ess_bdr_arg[i];
69  }
70 
71  sigma = new MeshDependentCoefficient(sigmaAttMap);
72  Tcapacity = new MeshDependentCoefficient(TcapacityAttMap);
73  InvTcap = new MeshDependentCoefficient(InvTcapAttMap);
74  InvTcond = new MeshDependentCoefficient(InvTcondAttMap);
75 
76  this->buildA0(*sigma);
77  this->buildM3(*Tcapacity);
78  this->buildM1(*sigma);
79  this->buildM2(*InvTcond);
80  this->buildS2(*InvTcap);
81  this->buildS1(1.0/mu);
82  this->buildCurl(1.0/mu);
83  this->buildDiv(*InvTcap);
84  this->buildGrad();
85 
89  A0 = new HypreParMatrix;
90  A1 = new HypreParMatrix;
91  A2 = new HypreParMatrix;
92  X0 = new Vector;
93  X1 = new Vector;
94  X2 = new Vector;
95  B0 = new Vector;
96  B1 = new Vector;
97  B2 = new Vector;
98  B3 = new Vector;
99 }
100 
102 {
103  Vector zero_vec(3); zero_vec = 0.0;
104  VectorConstantCoefficient Zero_vec(zero_vec);
105  ConstantCoefficient Zero(0.0);
106 
107  // The big BlockVector stores the fields as follows:
108  // Temperature
109  // Temperature Flux
110  // P field
111  // E field
112  // B field
113  // Joule Heating
114 
115  int Vsize_l2 = L2FESpace.GetVSize();
116  int Vsize_nd = HCurlFESpace.GetVSize();
117  int Vsize_rt = HDivFESpace.GetVSize();
118  int Vsize_h1 = HGradFESpace.GetVSize();
119 
120  Array<int> true_offset(7);
121  true_offset[0] = 0;
122  true_offset[1] = true_offset[0] + Vsize_l2;
123  true_offset[2] = true_offset[1] + Vsize_rt;
124  true_offset[3] = true_offset[2] + Vsize_h1;
125  true_offset[4] = true_offset[3] + Vsize_nd;
126  true_offset[5] = true_offset[4] + Vsize_rt;
127  true_offset[6] = true_offset[5] + Vsize_l2;
128 
129  Vector* xptr = (Vector*) &X;
130  ParGridFunction E, B, T, F, W, P;
131  T.MakeRef(&L2FESpace, *xptr,true_offset[0]);
132  F.MakeRef(&HDivFESpace, *xptr,true_offset[1]);
133  P.MakeRef(&HGradFESpace,*xptr,true_offset[2]);
134  E.MakeRef(&HCurlFESpace,*xptr,true_offset[3]);
135  B.MakeRef(&HDivFESpace, *xptr,true_offset[4]);
136  W.MakeRef(&L2FESpace, *xptr,true_offset[5]);
137 
138  E.ProjectCoefficient(Zero_vec);
139  B.ProjectCoefficient(Zero_vec);
140  F.ProjectCoefficient(Zero_vec);
141  T.ProjectCoefficient(Zero);
142  P.ProjectCoefficient(Zero);
143  W.ProjectCoefficient(Zero);
144 }
145 
146 /*
147 This is an experimental Mult() method for explicit integration.
148 Not recommended for actual use.
149 
150 S0 P = 0
151 M1 E = WeakCurl^T B + Grad P
152  dB = -Curl E
153 M2 F = WeakDiv^T T
154 M3 dT = WeakDiv F + W
155 
156 where W is the Joule heating.
157 
158 Boundary conditions are applied to E. No boundary conditions are applied to B.
159 Since we are using Hdiv, zero flux is an essential BC on F. P is given by Div
160 sigma Grad P = 0 with appropriate BC's.
161 */
162 void MagneticDiffusionEOperator::Mult(const Vector &X, Vector &dX_dt) const
163 {
164  dX_dt = 0.0;
165 
166  // The big BlockVector stores the fields as follows:
167  // Temperature
168  // Temperature Flux
169  // P field
170  // E field
171  // B field
172  // Joule Heating
173 
174  int Vsize_l2 = L2FESpace.GetVSize();
175  int Vsize_nd = HCurlFESpace.GetVSize();
176  int Vsize_rt = HDivFESpace.GetVSize();
177  int Vsize_h1 = HGradFESpace.GetVSize();
178 
179  Array<int> true_offset(7);
180  true_offset[0] = 0;
181  true_offset[1] = true_offset[0] + Vsize_l2;
182  true_offset[2] = true_offset[1] + Vsize_rt;
183  true_offset[3] = true_offset[2] + Vsize_h1;
184  true_offset[4] = true_offset[3] + Vsize_nd;
185  true_offset[5] = true_offset[4] + Vsize_rt;
186  true_offset[6] = true_offset[5] + Vsize_l2;
187 
188  Vector* xptr = (Vector*) &X;
189  ParGridFunction E, B, T, F, W, P;
190  T.MakeRef(&L2FESpace, *xptr,true_offset[0]);
191  F.MakeRef(&HDivFESpace, *xptr,true_offset[1]);
192  P.MakeRef(&HGradFESpace,*xptr,true_offset[2]);
193  E.MakeRef(&HCurlFESpace,*xptr,true_offset[3]);
194  B.MakeRef(&HDivFESpace, *xptr,true_offset[4]);
195  W.MakeRef(&L2FESpace, *xptr,true_offset[5]);
196 
197  ParGridFunction dE, dB, dT, dF, dW, dP;
198  dT.MakeRef(&L2FESpace, dX_dt,true_offset[0]);
199  dF.MakeRef(&HDivFESpace, dX_dt,true_offset[1]);
200  dP.MakeRef(&HGradFESpace,dX_dt,true_offset[2]);
201  dE.MakeRef(&HCurlFESpace,dX_dt,true_offset[3]);
202  dB.MakeRef(&HDivFESpace, dX_dt,true_offset[4]);
203  dW.MakeRef(&L2FESpace, dX_dt,true_offset[5]);
204 
205  // db = - Curl E
206  curl->Mult(E, dB);
207  dB *= -1.0;
208 
209  // form the Laplacian and solve it
210  ParGridFunction Phi_gf(&HGradFESpace);
211 
212  // p_bc is given function defining electrostatic potential on surface
213  FunctionCoefficient voltage(p_bc);
214  voltage.SetTime(this->GetTime());
215  Phi_gf = 0.0;
216 
217  // the line below is currently not fully supported on AMR meshes
218  // Phi_gf.ProjectBdrCoefficient(voltage,poisson_ess_bdr);
219 
220  // this is a hack to get around the above issue
221  Phi_gf.ProjectCoefficient(voltage);
222  // end of hack
223 
224  // apply essential BC's and apply static condensation, the new system to
225  // solve is A0 X0 = B0
226  Array<int> poisson_ess_tdof_list;
227  HGradFESpace.GetEssentialTrueDofs(poisson_ess_bdr, poisson_ess_tdof_list);
228 
229  *v0 = 0.0;
230  a0->FormLinearSystem(poisson_ess_tdof_list,Phi_gf,*v0,*A0,*X0,*B0);
231 
232  if (amg_a0 == NULL) { amg_a0 = new HypreBoomerAMG(*A0); }
233  if (pcg_a0 == NULL)
234  {
235  pcg_a0 = new HyprePCG(*A0);
240  }
241  // pcg "Mult" operation is a solve
242  // X0 = A0^-1 * B0
243  pcg_a0->Mult(*B0, *X0);
244 
245  // "undo" the static condensation using dP as a temporary variable, dP stores
246  // Pnew
247  a0->RecoverFEMSolution(*X0,*v0,P);
248  dP = 0.0;
249 
250  // v1 = <1/mu v, curl u> B
251  // B is a grid function but weakCurl is not parallel assembled so is OK
252  weakCurl->MultTranspose(B, *v1);
253 
254  // now add Grad dPhi/dt term
255  // use E as a temporary, E = Grad P
256  // v1 = curl 1/mu B + M1 * Grad P
257  // note: these two steps could be replaced by one step if we have the
258  // bilinear form <sigma gradP, E>
259  grad->Mult(P,E);
260  m1->AddMult(E,*v1,1.0);
261 
262  // OK now v1 is the right hand side, just need to add essential BC's
263 
265 
266  // edot_bc is time-derivative E-field on a boundary surface and then it is
267  // used as a Dirichlet BC.
268 
270  J_gf = 0.0;
272 
273  // apply essential BC's and apply static condensation
274  // the new system to solve is M1 X1 = B1
275  Array<int> ess_tdof_list;
276  HCurlFESpace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
277 
278  m1->FormLinearSystem(ess_tdof_list,J_gf,*v1,*A1,*X1,*B1);
279 
280  if (dsp_m1 == NULL) { dsp_m1 = new HypreDiagScale(*A1); }
281  if (pcg_m1 == NULL)
282  {
283  pcg_m1 = new HyprePCG(*A1);
288  }
289  // pcg "Mult" operation is a solve
290  // X1 = M1^-1 * B1 = M1^-1 (-S1 E)
291  pcg_m1->Mult(*B1, *X1);
292 
293  // "undo" the static condensation and fill in grid function dE
294  m1->RecoverFEMSolution(*X1,*v1,E);
295  dE = 0.0;
296 
297  // the total field is E_tot = E_ind - Grad Phi
298  // so we need to subtract out Grad Phi
299  // E = E - grad (P)
300  grad->AddMult(P,E,-1.0);
301 
302  // Compute Joule heating using the previous value of E
303  this->GetJouleHeating(E,W);
304  dW = 0.0;
305 
306  // Mult(x,y,alpha=1,beta=0)
307  // y = alpha*A*x + beta*y
308  // giving
309  // v2 = <v, div u> * T
310  weakDiv->Mult(T, *v2);
311 
312  // apply the thermal BC
313  // recall for Hdiv formulation the essential BC is on the flux
314  Vector zero_vec(3); zero_vec = 0.0;
315  VectorConstantCoefficient Zero_vec(zero_vec);
317  F_gf = 0.0;
319 
320  // apply essential BC's and apply static condensation
321  // the new system to solve is M2 X2 = B2
322  Array<int> thermal_ess_tdof_list;
323  HDivFESpace.GetEssentialTrueDofs(thermal_ess_bdr, thermal_ess_tdof_list);
324 
325  m2->FormLinearSystem(thermal_ess_tdof_list,F_gf,*v2,*A2,*X2,*B2);
326 
327  if (dsp_m2 == NULL) { dsp_m2 = new HypreDiagScale(*A2); }
328  if (pcg_m2 == NULL)
329  {
330  pcg_m2 = new HyprePCG(*A2);
335  }
336  // X2 = m2^-1 * B2
337  pcg_m2->Mult(*B2, *X2);
338 
339  // "undo" the static condensation and fill in grid function dF
340  m2->RecoverFEMSolution(*X2,*v2,F);
341 
342  // Compute dT using previous value of flux
343  // dT = [w - div F]
344  //
345  // <u,u> dT = <1/c W,u> - <1/c div v,u> F
346  //
347  // where W is Joule heating and F is the flux that we just computed
348  //
349  // note: if div is a BilinearForm, then W should be converted to a LoadVector
350 
351  GridFunctionCoefficient Wcoeff(&W);
352  ParLinearForm temp_lf(&L2FESpace);
353 
354  // compute load vector < W, u>
355  temp_lf.AddDomainIntegrator(new DomainLFIntegrator(Wcoeff));
356  temp_lf.Assemble();
357  // lf = lf - div F
358  weakDiv->AddMult(F, temp_lf, -1.0);
359 
360  // if div is a BilinearForm, need to perform mass matrix solve to convert
361  // energy cT to temperature T
362 
363  if (dsp_m3 == NULL) { dsp_m3 = new HypreDiagScale(*M3); }
364  if (pcg_m3 == NULL)
365  {
366  pcg_m3 = new HyprePCG(*M3);
371  }
372  // solve for dT from M3 dT = lf
373  // no boundary conditions on this solve
374  pcg_m3->Mult(temp_lf, dT);
375 }
376 
377 /*
378 This is the main computational code that computes dX/dt implicitly
379 where X is the state vector containing P, E, B, F, T, and W
380 
381  S0 P = 0
382 (M1+dt S1) E = WeakCurl^T B + Grad P
383  dB = -Curl E
384 (M2+dt S2) F = WeakDiv^T T
385  M3 dT = WeakDiv F + W
386 
387 where W is the Joule heating.
388 
389 Boundary conditions are applied to E. Boundary conditions are applied to F. No
390 boundary conditions are applied to B or T.
391 
392 The W term in the left hand side is the Joule heating which is a nonlinear
393 (quadratic) function of E.
394 
395 P is solution of Div sigma Grad dP = 0.
396 
397 The total E-field is given by E_tot = E_ind - Grad P, the big equation for E
398 above is really for E_ind (the induced, or solenoidal, component) and this is
399 corrected for.
400 */
402  const Vector &X, Vector &dX_dt)
403 {
404  if ( A2 == NULL || fabs(dt-dt_A2) > 1.0e-12*dt )
405  {
406  this->buildA2(*InvTcond, *InvTcap, dt);
407  }
408  if ( A1 == NULL || fabs(dt-dt_A1) > 1.0e-12*dt )
409  {
410  this->buildA1(1.0/mu, *sigma, dt);
411  }
412 
413  dX_dt = 0.0;
414 
415  // The big BlockVector stores the fields as follows:
416  // Temperature
417  // Temperature Flux
418  // P field
419  // E field
420  // B field
421  // Joule Heating
422 
423  int Vsize_l2 = L2FESpace.GetVSize();
424  int Vsize_nd = HCurlFESpace.GetVSize();
425  int Vsize_rt = HDivFESpace.GetVSize();
426  int Vsize_h1 = HGradFESpace.GetVSize();
427 
428  Array<int> true_offset(7);
429  true_offset[0] = 0;
430  true_offset[1] = true_offset[0] + Vsize_l2;
431  true_offset[2] = true_offset[1] + Vsize_rt;
432  true_offset[3] = true_offset[2] + Vsize_h1;
433  true_offset[4] = true_offset[3] + Vsize_nd;
434  true_offset[5] = true_offset[4] + Vsize_rt;
435  true_offset[6] = true_offset[5] + Vsize_l2;
436 
437  Vector* xptr = (Vector*) &X;
438  ParGridFunction E, B, T, F, W, P;
439  T.MakeRef(&L2FESpace, *xptr,true_offset[0]);
440  F.MakeRef(&HDivFESpace, *xptr,true_offset[1]);
441  P.MakeRef(&HGradFESpace,*xptr,true_offset[2]);
442  E.MakeRef(&HCurlFESpace,*xptr,true_offset[3]);
443  B.MakeRef(&HDivFESpace, *xptr,true_offset[4]);
444  W.MakeRef(&L2FESpace, *xptr,true_offset[5]);
445 
446  ParGridFunction dE, dB, dT, dF, dW, dP;
447  dT.MakeRef(&L2FESpace, dX_dt,true_offset[0]);
448  dF.MakeRef(&HDivFESpace, dX_dt,true_offset[1]);
449  dP.MakeRef(&HGradFESpace,dX_dt,true_offset[2]);
450  dE.MakeRef(&HCurlFESpace,dX_dt,true_offset[3]);
451  dB.MakeRef(&HDivFESpace, dX_dt,true_offset[4]);
452  dW.MakeRef(&L2FESpace, dX_dt,true_offset[5]);
453 
454  // form the Laplacian and solve it
455  ParGridFunction Phi_gf(&HGradFESpace);
456 
457  // p_bc is given function defining electrostatic potential on surface
458  FunctionCoefficient voltage(p_bc);
459  voltage.SetTime(this->GetTime());
460  Phi_gf = 0.0;
461 
462  // the function below is currently not fully supported on AMR meshes
463  // Phi_gf.ProjectBdrCoefficient(voltage,poisson_ess_bdr);
464 
465  // this is a hack to get around the above issue
466  Phi_gf.ProjectCoefficient(voltage);
467  // end of hack
468 
469  // apply essential BC's and apply static condensation, the new system to
470  // solve is A0 X0 = B0
471  Array<int> poisson_ess_tdof_list;
472  HGradFESpace.GetEssentialTrueDofs(poisson_ess_bdr, poisson_ess_tdof_list);
473 
474  *v0 = 0.0;
475  a0->FormLinearSystem(poisson_ess_tdof_list,Phi_gf,*v0,*A0,*X0,*B0);
476 
477  if (amg_a0 == NULL) { amg_a0 = new HypreBoomerAMG(*A0); }
478  if (pcg_a0 == NULL)
479  {
480  pcg_a0 = new HyprePCG(*A0);
485  }
486  // pcg "Mult" operation is a solve
487  // X0 = A0^-1 * B0
488  pcg_a0->Mult(*B0, *X0);
489 
490  // "undo" the static condensation saving result in grid function dP
491  a0->RecoverFEMSolution(*X0,*v0,P);
492  dP = 0.0;
493 
494  // v1 = <1/mu v, curl u> B
495  // B is a grid function but weakCurl is not parallel assembled so is OK
496  weakCurl->MultTranspose(B, *v1);
497 
498  // now add Grad dPhi/dt term
499  // use E as a temporary, E = Grad P
500  // v1 = curl 1/mu B + M1 * Grad P
501  grad->Mult(P,E);
502  m1->AddMult(E,*v1,1.0);
503 
505 
506  // edot_bc is time-derivative E-field on a boundary surface
507  // and then it is used as a Dirichlet BC
508  // the vector v1 will be modified by the values Jtmp and
509  // the part of the matrix m1 that hs been eliminated (but stored).
511  J_gf = 0.0;
513 
514  // form the linear system, including eliminating essential BC's and applying
515  // static condensation. The system to solve is A1 X1 = B1
516  Array<int> ess_tdof_list;
517  HCurlFESpace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
518 
519  a1->FormLinearSystem(ess_tdof_list,J_gf,*v1,*A1,*X1,*B1);
520 
521  // We only need to create the solver and preconditioner once
522  if ( ams_a1 == NULL )
523  {
524  ParFiniteElementSpace *prec_fespace =
526  ams_a1 = new HypreAMS(*A1, prec_fespace);
527  }
528  if ( pcg_a1 == NULL )
529  {
530  pcg_a1 = new HyprePCG(*A1);
535  }
536  // solve the system
537  // dE = (A1)^-1 [-S1 E]
538  pcg_a1->Mult(*B1, *X1);
539 
540  // this is required because of static condensation, E is a grid function
541  a1->RecoverFEMSolution(*X1,*v1,E);
542  dE = 0.0;
543 
544  // the total field is E_tot = E_ind - Grad Phi
545  // so we need to subtract out Grad Phi
546  // E = E - grad (P)
547  // note grad maps GF to GF
548  grad->AddMult(P,E,-1.0);
549 
550  // Compute dB/dt = -Curl(E_{n+1})
551  // note curl maps GF to GF
552  curl->Mult(E, dB);
553  dB *= -1.0;
554 
555  // Compute Energy Deposition
556  this->GetJouleHeating(E,W);
557 
558  // v2 = Div^T * W, where W is the Joule heating computed above, and
559  // Div is the matrix <div u, v>
560  weakDivC->MultTranspose(W, *v2);
561  *v2 *= dt;
562 
563  // v2 = <v, div u> T + (1.0)*v2
564  weakDiv->AddMultTranspose(T, *v2, 1.0);
565 
566  // apply the thermal BC
567  Vector zero_vec(3); zero_vec = 0.0;
568  VectorConstantCoefficient Zero_vec(zero_vec);
570  F_gf = 0.0;
572 
573  // form the linear system, including eliminating essential BC's and applying
574  // static condensation. The system to solve is A2 X2 = B2
575  Array<int> thermal_ess_tdof_list;
576  HDivFESpace.GetEssentialTrueDofs(thermal_ess_bdr, thermal_ess_tdof_list);
577  a2->FormLinearSystem(thermal_ess_tdof_list,F_gf,*v2,*A2,*X2,*B2);
578 
579  // We only need to create the solver and preconditioner once
580  if ( ads_a2 == NULL )
581  {
582  ParFiniteElementSpace *prec_fespace =
584  ads_a2 = new HypreADS(*A2, prec_fespace);
585  }
586  if ( pcg_a2 == NULL )
587  {
588  pcg_a2 = new HyprePCG(*A2);
593  }
594  // solve for dF from a2 dF = v2
595  // dF = (A2)^-1 [S2*F + rhs]
596  pcg_a2->Mult(*B2, *X2);
597 
598  // this is required because of static condensation
599  a2->RecoverFEMSolution(*X2,*v2,F);
600 
601  // c dT = [W - div F]
602  //
603  // <u,u> dT = <1/c W,u> - <1/c div v,u>
604  //
605  // where W is Joule heating and F is the flux that we just computed
606  //
607  // note: if div is a BilinearForm, then W should be converted to a LoadVector
608  // compute load vector <1/c W, u> where W is the Joule heating GF
609 
610  // create the Coefficient 1/c W
611  //ScaledGFCoefficient Wcoeff(&W, *InvTcap);
612  GridFunctionCoefficient Wcoeff(&W);
613 
614  // compute <W,u>
615  ParLinearForm temp_lf(&L2FESpace);
616  temp_lf.AddDomainIntegrator(new DomainLFIntegrator(Wcoeff));
617  temp_lf.Assemble();
618 
619  // lf = lf - div F
620  weakDiv->AddMult(F, temp_lf, -1.0);
621 
622  // need to perform mass matrix solve to get temperature T
623  // <c u, u> Tdot = -<div v, u> F + <1/c W, u>
624  // NOTE: supposedly we can just invert any L2 matrix, could do that here
625  // instead of a solve
626 
627  if (dsp_m3 == NULL) { dsp_m3 = new HypreDiagScale(*M3); }
628  if (pcg_m3 == NULL)
629  {
630  pcg_m3 = new HyprePCG(*M3);
635  }
636 
637  // solve for dT from M3 dT = lf
638  // no boundary conditions on this solve
639  pcg_m3->Mult(temp_lf, dT);
640 }
641 
643 {
644  if ( a0 != NULL ) { delete a0; }
645 
646  // First create and assemble the bilinear form. For now we assume the mesh
647  // isn't moving, the materials are time independent, and dt is constant. So
648  // we only need to do this once.
649 
650  // ConstantCoefficient Sigma(sigma);
653  if (STATIC_COND == 1) { a0->EnableStaticCondensation(); }
654  a0->Assemble();
655 
656  // Don't finalize or parallel assemble this is done in FormLinearSystem.
657 }
658 
661  double dt)
662 {
663  if ( a1 != NULL ) { delete a1; }
664 
665  // First create and assemble the bilinear form. For now we assume the mesh
666  // isn't moving, the materials are time independent, and dt is constant. So
667  // we only need to do this once.
668 
669  ConstantCoefficient dtMuInv(dt*muInv);
673  if (STATIC_COND == 1) { a1->EnableStaticCondensation(); }
674  a1->Assemble();
675 
676  // Don't finalize or parallel assemble this is done in FormLinearSystem.
677 
678  dt_A1 = dt;
679 }
680 
682  MeshDependentCoefficient &InvTcap_,
683  double dt_)
684 {
685  if ( a2 != NULL ) { delete a2; }
686 
687  InvTcap_.SetScaleFactor(dt_);
690  a2->AddDomainIntegrator(new DivDivIntegrator(InvTcap_));
691  if (STATIC_COND == 1) { a2->EnableStaticCondensation(); }
692  a2->Assemble();
693 
694  // Don't finalize or parallel assemble this is done in FormLinearSystem.
695 
696  dt_A2 = dt_;
697 }
698 
700 {
701  if ( m1 != NULL ) { delete m1; }
702 
705  m1->Assemble();
706 
707  // Don't finalize or parallel assemble this is done in FormLinearSystem.
708 }
709 
711 {
712  if ( m2 != NULL ) { delete m2; }
713 
714  // ConstantCoefficient MuInv(muInv);
717  m2->Assemble();
718 
719  // Don't finalize or parallel assemble this is done in FormLinearSystem.
720 }
721 
723 {
724  if ( m3 != NULL ) { delete m3; }
725 
726  // ConstantCoefficient Sigma(sigma);
727  m3 = new ParBilinearForm(&L2FESpace);
728  m3->AddDomainIntegrator(new MassIntegrator(Tcapacity_));
729  m3->Assemble();
730  m3->Finalize();
731  M3 = m3->ParallelAssemble();
732 }
733 
735 {
736  if ( s1 != NULL ) { delete s1; }
737 
738  ConstantCoefficient MuInv(muInv);
741  s1->Assemble();
742 }
743 
745 {
746  if ( s2 != NULL ) { delete s2; }
747 
748  // ConstantCoefficient param(a);
750  s2->AddDomainIntegrator(new DivDivIntegrator(InvTcap_));
751  s2->Assemble();
752 }
753 
755 {
756  if ( curl != NULL ) { delete curl; }
757  if ( weakCurl != NULL ) { delete weakCurl; }
758 
761  curl->Assemble();
762 
763  ConstantCoefficient MuInv(muInv);
766  weakCurl->Assemble();
767 
768  // no ParallelAssemble since this will be applied to GridFunctions
769 }
770 
772 {
773  if ( weakDiv != NULL ) { delete weakDiv; }
774  if ( weakDivC != NULL ) { delete weakDivC; }
775 
778  weakDivC->Assemble();
779 
782  weakDiv->Assemble();
783 
784  // no ParallelAssemble since this will be applied to GridFunctions
785 }
786 
788 {
789  if ( grad != NULL ) { delete grad; }
790 
793  grad->Assemble();
794 
795  // no ParallelAssemble since this will be applied to GridFunctions
796 }
797 
799 {
800  double el = m1->InnerProduct(E_gf,E_gf);
801 
802  double global_el;
803  MPI_Allreduce(&el, &global_el, 1, MPI_DOUBLE, MPI_SUM,
804  m2->ParFESpace()->GetComm());
805 
806  return el;
807 }
808 
809 // E is the input GF, w is the output GF which is assumed to be an L2 scalar
810 // representing the Joule heating
812  ParGridFunction &w_gf) const
813 {
814  // The w_coeff object stashes a reference to sigma and E, and it has
815  // an Eval method that will be used by ProjectCoefficient.
816  JouleHeatingCoefficient w_coeff(*sigma, E_gf);
817 
818  // This applies the definition of the finite element degrees-of-freedom
819  // to convert the function to a set of discrete values
820  w_gf.ProjectCoefficient(w_coeff);
821 }
822 
824 { t = t_; }
825 
827 {
828  if ( ams_a1 != NULL ) { delete ams_a1; }
829  if ( pcg_a1 != NULL ) { delete pcg_a1; }
830 
831  if ( dsp_m1 != NULL ) { delete dsp_m1; }
832  if ( pcg_m1 != NULL ) { delete pcg_m1; }
833 
834  if ( dsp_m2 != NULL ) { delete dsp_m2; }
835  if ( pcg_m2 != NULL ) { delete pcg_m2; }
836 
837  if ( curl != NULL ) { delete curl; }
838  if ( weakDiv != NULL ) { delete weakDiv; }
839  if ( weakDivC != NULL ) { delete weakDivC; }
840  if ( weakCurl != NULL ) { delete weakCurl; }
841  if ( grad != NULL ) { delete grad; }
842 
843  if ( a0 != NULL ) { delete a0; }
844  if ( a1 != NULL ) { delete a1; }
845  if ( a2 != NULL ) { delete a2; }
846  if ( m1 != NULL ) { delete m1; }
847  if ( m2 != NULL ) { delete m2; }
848  if ( s1 != NULL ) { delete s1; }
849  if ( s2 != NULL ) { delete s2; }
850 
851  if ( A0 != NULL ) { delete A0; }
852  if ( X0 != NULL ) { delete X0; }
853  if ( B0 != NULL ) { delete B0; }
854 
855  if ( A1 != NULL ) { delete A1; }
856  if ( X1 != NULL ) { delete X1; }
857  if ( B1 != NULL ) { delete B1; }
858 
859  if ( A2 != NULL ) { delete A2; }
860  if ( X2 != NULL ) { delete X2; }
861  if ( B2 != NULL ) { delete B2; }
862 
863  if ( v1 != NULL ) { delete v1; }
864  if ( v2 != NULL ) { delete v2; }
865 
866  if (sigma != NULL) { delete sigma; }
867  if (Tcapacity != NULL) { delete Tcapacity; }
868  if (InvTcap != NULL) { delete InvTcap; }
869  if (InvTcond != NULL) { delete InvTcond; }
870 
871  delete amg_a0;
872  delete pcg_a0;
873  delete pcg_a2;
874  delete ads_a2;
875  delete m3;
876  delete dsp_m3;
877  delete pcg_m3;
878  delete M1;
879  delete M2;
880  delete M3;
881  delete v0;
882  delete B3;
883 }
884 
885 void MagneticDiffusionEOperator::Debug(const char *base, double)
886 {
887  {
888  hypre_ParCSRMatrixPrint(*A1,"A1_");
889  HypreParVector tempB1(A1->GetComm(),A1->N(),B1->GetData(),A1->ColPart());
890  tempB1.Print("B1_");
891  HypreParVector tempX1(A1->GetComm(),A1->N(),X1->GetData(),A1->ColPart());
892  tempX1.Print("X1_");
893  }
894 
895  {
896  hypre_ParCSRMatrixPrint(*A2,"A2_");
897  HypreParVector tempB2(A2->GetComm(),A2->N(),B2->GetData(),A2->ColPart());
898  tempB2.Print("B2_");
899  HypreParVector tempX2(A2->GetComm(),A2->N(),X2->GetData(),A2->ColPart());
900  tempX2.Print("X2_");
901  }
902 }
903 
905  const IntegrationPoint &ip)
906 {
907  Vector E;
908  double thisSigma;
909  E_gf.GetVectorValue(T, ip, E);
910  thisSigma = sigma.Eval(T, ip);
911  return thisSigma*(E*E);
912 }
913 
915  const std::map<int, double> &inputMap, double scale)
916  : Coefficient()
917 {
918  // make a copy of the magic attribute-value map for later use
919  materialMap = new std::map<int, double>(inputMap);
920  scaleFactor = scale;
921 }
922 
924  const MeshDependentCoefficient &cloneMe)
925  : Coefficient()
926 {
927  // make a copy of the magic attribute-value map for later use
928  materialMap = new std::map<int, double>(*(cloneMe.materialMap));
929  scaleFactor = cloneMe.scaleFactor;
930 }
931 
933  const IntegrationPoint &ip)
934 {
935  // given the attribute, extract the coefficient value from the map
936  std::map<int, double>::iterator it;
937  int thisAtt = T.Attribute;
938  double value;
939  it = materialMap->find(thisAtt);
940  if (it != materialMap->end())
941  {
942  value = it->second;
943  }
944  else
945  {
946  value = 0.0; // avoid compile warning
947  std::cerr << "MeshDependentCoefficient attribute " << thisAtt
948  << " not found" << std::endl;
949  mfem_error();
950  }
951 
952  return value*scaleFactor;
953 }
954 
956  MeshDependentCoefficient &input_mdc)
957  : GridFunctionCoefficient(gf), mdc(input_mdc) {}
958 
960  const IntegrationPoint &ip)
961 {
962  return mdc.Eval(T,ip) * GridFunctionCoefficient::Eval(T,ip);
963 }
964 
965 } // namespace electromagnetics
966 
967 } // namespace mfem
968 
969 #endif // MFEM_USE_MPI
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
void buildDiv(MeshDependentCoefficient &InvTcap)
void SetTol(double tol)
Definition: hypre.cpp:3995
virtual double GetTime() const
Read the currently set time.
Definition: operator.hpp:330
int Size() const
Return the logical size of the array.
Definition: array.hpp:138
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:108
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1032
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:534
ParFiniteElementSpace * SCParFESpace() const
Return the parallel trace FE space associated with static condensation.
void buildM2(MeshDependentCoefficient &alpha)
The Auxiliary-space Maxwell Solver in hypre.
Definition: hypre.hpp:1732
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
The Auxiliary-space Divergence Solver in hypre.
Definition: hypre.hpp:1809
A coefficient that is constant across space and time.
Definition: coefficient.hpp:84
virtual void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
Form the linear system A X = B, corresponding to this bilinear form and the linear form b(...
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Add the matrix vector multiple to a vector: .
Integrator for (curl u, curl v) for Nedelec elements.
Base abstract class for first order time dependent operators.
Definition: operator.hpp:289
Vector coefficient that is constant in space and time.
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
virtual void ImplicitSolve(const double dt, const Vector &x, Vector &k)
Solve the equation: k = f(x + dt k, t), for the unknown k at the current time t.
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:525
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
virtual void SetTime(double t)
Set the time for time dependent coefficients.
Definition: coefficient.hpp:50
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:209
(Q div u, div v) for RT elements
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:4015
void AddDomainInterpolator(DiscreteInterpolator *di)
Adds a domain interpolator. Assumes ownership of di.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Debug(const char *basefilename, double time)
HYPRE_BigInt N() const
Returns the global number of columns.
Definition: hypre.hpp:585
double t
Current time.
Definition: operator.hpp:313
void GetJouleHeating(ParGridFunction &E_gf, ParGridFunction &w_gf) const
void Print(const char *fname) const
Prints the locally owned rows in parallel.
Definition: hypre.cpp:387
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds a domain integrator. Assumes ownership of bfi.
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Definition: pgridfunc.cpp:102
void buildA1(double muInv, MeshDependentCoefficient &sigma, double dt)
ScaledGFCoefficient(GridFunction *gf, MeshDependentCoefficient &input_mdc)
double p_bc(const Vector &x, double t)
Definition: joule.cpp:758
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1579
Class for parallel linear form.
Definition: plinearform.hpp:26
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Definition: gridfunc.cpp:2678
MPI_Comm GetComm() const
Definition: pfespace.hpp:273
virtual void AddMultTranspose(const Vector &x, Vector &y, const double a=1.0) const
void buildS2(MeshDependentCoefficient &alpha)
void buildM3(MeshDependentCoefficient &Tcap)
Jacobi preconditioner in hypre.
Definition: hypre.hpp:1410
ParFiniteElementSpace * ParFESpace() const
Return the parallel FE space associated with the ParBilinearForm.
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
void buildM1(MeshDependentCoefficient &sigma)
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Definition: plinearform.cpp:46
void Assemble(int skip_zeros=1)
double muInv(const Vector &x)
Definition: maxwell.cpp:96
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetMaxIter(int max_iter)
Definition: hypre.cpp:4005
A general vector function coefficient.
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: pgridfunc.cpp:321
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Definition: linearform.cpp:41
double ElectricLosses(ParGridFunction &E_gf) const
bool StaticCondensationIsEnabled() const
Check if static condensation was actually enabled by a previous call to EnableStaticCondensation().
MeshDependentCoefficient(const std::map< int, double > &inputMap, double scale=1.0)
PCG solver in hypre.
Definition: hypre.hpp:1204
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:679
double InnerProduct(const Vector &x, const Vector &y) const
Compute .
void SetTime(const double t_)
Set the current time.
void buildA2(MeshDependentCoefficient &InvTcond, MeshDependentCoefficient &InvTcap, double dt)
virtual void RecoverFEMSolution(const Vector &X, const Vector &b, Vector &x)
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
Class for integration point with weight.
Definition: intrules.hpp:25
Class for parallel bilinear form using different test and trial FE spaces.
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization.
HYPRE_BigInt * ColPart()
Returns the column partitioning.
Definition: hypre.hpp:573
double mu
Definition: ex25.cpp:139
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attri...
Definition: pgridfunc.cpp:671
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:4020
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Class for parallel bilinear form.
void edot_bc(const Vector &x, Vector &E)
Definition: joule.cpp:733
A general function coefficient.
Vector data type.
Definition: vector.hpp:60
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Class for parallel grid function.
Definition: pgridfunc.hpp:32
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient at ip.
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
virtual void Assemble(int skip_zeros=1)
Construct the internal matrix representation of the discrete linear operator.
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:4043
void buildA0(MeshDependentCoefficient &sigma)
virtual void Mult(const Vector &x, Vector &y) const
Matrix multiplication: .
virtual void Mult(const Vector &vx, Vector &dvx_dt) const
Perform the action of the operator: y = k = f(x, t), where k solves the algebraic equation F(x...
void EnableStaticCondensation()
Enable the use of static condensation. For details see the description for class StaticCondensation i...