MFEM  v3.3.2
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, 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 "joule_solver.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 using namespace std;
17 
18 namespace mfem
19 {
20 
21 using namespace miniapps;
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), X0(NULL), X1(NULL), X2(NULL), B0(NULL),
48  B1(NULL), B2(NULL), B3(NULL),
49  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  {
222  Array<int> my_ess_dof_list;
225  homer.ProjectCoefficient(voltage);
226  for (int i = 0; i < my_ess_dof_list.Size(); i++)
227  {
228  int dof = my_ess_dof_list[i];
229  if (dof < 0) { dof = -1 - dof; }
230  Phi_gf[i] = homer[i];
231  }
232  }
233  // end of hack
234 
235  // apply essential BC's and apply static condensation, the new system to
236  // solve is A0 X0 = B0
237  Array<int> poisson_ess_tdof_list;
238  HGradFESpace.GetEssentialTrueDofs(poisson_ess_bdr, poisson_ess_tdof_list);
239 
240  *v0 = 0.0;
241  a0->FormLinearSystem(poisson_ess_tdof_list,Phi_gf,*v0,*A0,*X0,*B0);
242 
243  if (amg_a0 == NULL) { amg_a0 = new HypreBoomerAMG(*A0); }
244  if (pcg_a0 == NULL)
245  {
246  pcg_a0 = new HyprePCG(*A0);
247  pcg_a0->SetTol(SOLVER_TOL);
248  pcg_a0->SetMaxIter(SOLVER_MAX_IT);
249  pcg_a0->SetPrintLevel(SOLVER_PRINT_LEVEL);
251  }
252  // pcg "Mult" operation is a solve
253  // X0 = A0^-1 * B0
254  pcg_a0->Mult(*B0, *X0);
255 
256  // "undo" the static condensation using dP as a temporary variable, dP stores
257  // Pnew
258  a0->RecoverFEMSolution(*X0,*v0,P);
259  dP = 0.0;
260 
261  // v1 = <1/mu v, curl u> B
262  // B is a grid function but weakCurl is not parallel assembled so is OK
263  weakCurl->MultTranspose(B, *v1);
264 
265  // now add Grad dPhi/dt term
266  // use E as a temporary, E = Grad P
267  // v1 = curl 1/mu B + M1 * Grad P
268  // note: these two steps could be replaced by one step if we have the
269  // bilinear form <sigma gradP, E>
270  grad->Mult(P,E);
271  m1->AddMult(E,*v1,1.0);
272 
273  // OK now v1 is the right hand side, just need to add essential BC's
274 
276 
277  // edot_bc is time-derivative E-field on a boundary surface and then it is
278  // used as a Dirichlet BC.
279 
281  J_gf = 0.0;
283 
284  // apply essential BC's and apply static condensation
285  // the new system to solve is M1 X1 = B1
286  Array<int> ess_tdof_list;
287  HCurlFESpace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
288 
289  m1->FormLinearSystem(ess_tdof_list,J_gf,*v1,*A1,*X1,*B1);
290 
291  if (dsp_m1 == NULL) { dsp_m1 = new HypreDiagScale(*A1); }
292  if (pcg_m1 == NULL)
293  {
294  pcg_m1 = new HyprePCG(*A1);
295  pcg_m1->SetTol(SOLVER_TOL);
296  pcg_m1->SetMaxIter(SOLVER_MAX_IT);
297  pcg_m1->SetPrintLevel(SOLVER_PRINT_LEVEL);
299  }
300  // pcg "Mult" operation is a solve
301  // X1 = M1^-1 * B1 = M1^-1 (-S1 E)
302  pcg_m1->Mult(*B1, *X1);
303 
304  // "undo" the static condensation and fill in grid function dE
305  m1->RecoverFEMSolution(*X1,*v1,E);
306  dE = 0.0;
307 
308  // the total field is E_tot = E_ind - Grad Phi
309  // so we need to subtract out Grad Phi
310  // E = E - grad (P)
311  grad->AddMult(P,E,-1.0);
312 
313  // Compute Joule heating using the previous value of E
314  this->GetJouleHeating(E,W);
315  dW = 0.0;
316 
317  // Mult(x,y,alpha=1,beta=0)
318  // y = alpha*A*x + beta*y
319  // giving
320  // v2 = <v, div u> * T
321  weakDiv->Mult(T, *v2);
322 
323  // apply the thermal BC
324  // recall for Hdiv formulation the essential BC is on the flux
325  Vector zero_vec(3); zero_vec = 0.0;
326  VectorConstantCoefficient Zero_vec(zero_vec);
328  F_gf = 0.0;
330 
331  // apply essential BC's and apply static condensation
332  // the new system to solve is M2 X2 = B2
333  Array<int> thermal_ess_tdof_list;
334  HDivFESpace.GetEssentialTrueDofs(thermal_ess_bdr, thermal_ess_tdof_list);
335 
336  m2->FormLinearSystem(thermal_ess_tdof_list,F_gf,*v2,*A2,*X2,*B2);
337 
338  if (dsp_m2 == NULL) { dsp_m2 = new HypreDiagScale(*A2); }
339  if (pcg_m2 == NULL)
340  {
341  pcg_m2 = new HyprePCG(*A2);
342  pcg_m2->SetTol(SOLVER_TOL);
343  pcg_m2->SetMaxIter(SOLVER_MAX_IT);
344  pcg_m2->SetPrintLevel(SOLVER_PRINT_LEVEL);
346  }
347  // X2 = m2^-1 * B2
348  pcg_m2->Mult(*B2, *X2);
349 
350  // "undo" the static condensation and fill in grid function dF
351  m2->RecoverFEMSolution(*X2,*v2,F);
352 
353  // Compute dT using previous value of flux
354  // dT = [w - div F]
355  //
356  // <u,u> dT = <1/c W,u> - <1/c div v,u> F
357  //
358  // where W is Joule heating and F is the flux that we just computed
359  //
360  // note: if div is a BilinearForm, then W should be converted to a LoadVector
361 
362  GridFunctionCoefficient Wcoeff(&W);
363  ParLinearForm temp_lf(&L2FESpace);
364 
365  // compute load vector < W, u>
366  temp_lf.AddDomainIntegrator(new DomainLFIntegrator(Wcoeff));
367  temp_lf.Assemble();
368  // lf = lf - div F
369  weakDiv->AddMult(F, temp_lf, -1.0);
370 
371  // if div is a BilinearForm, need to perform mass matrix solve to convert
372  // energy cT to temperature T
373 
374  if (dsp_m3 == NULL) { dsp_m3 = new HypreDiagScale(*M3); }
375  if (pcg_m3 == NULL)
376  {
377  pcg_m3 = new HyprePCG(*M3);
378  pcg_m3->SetTol(SOLVER_TOL);
379  pcg_m3->SetMaxIter(SOLVER_MAX_IT);
380  pcg_m3->SetPrintLevel(SOLVER_PRINT_LEVEL);
382  }
383  // solve for dT from M3 dT = lf
384  // no boundary conditions on this solve
385  pcg_m3->Mult(temp_lf, dT);
386 }
387 
388 /*
389 This is the main computational code that computes dX/dt implicitly
390 where X is the state vector containing P, E, B, F, T, and W
391 
392  S0 P = 0
393 (M1+dt S1) E = WeakCurl^T B + Grad P
394  dB = -Curl E
395 (M2+dt S2) F = WeakDiv^T T
396  M3 dT = WeakDiv F + W
397 
398 where W is the Joule heating.
399 
400 Boundary conditions are applied to E. Boundary conditions are applied to F. No
401 boundary conditions are applied to B or T.
402 
403 The W term in the left hand side is the Joule heating which is a nonlinear
404 (quadratic) function of E.
405 
406 P is solution of Div sigma Grad dP = 0.
407 
408 The total E-field is given by E_tot = E_ind - Grad P, the big equation for E
409 above is really for E_ind (the induced, or solenoidal, component) and this is
410 corrected for.
411 */
413  const Vector &X, Vector &dX_dt)
414 {
415  if ( A2 == NULL || fabs(dt-dt_A2) > 1.0e-12*dt )
416  {
417  this->buildA2(*InvTcond, *InvTcap, dt);
418  }
419  if ( A1 == NULL || fabs(dt-dt_A1) > 1.0e-12*dt )
420  {
421  this->buildA1(1.0/mu, *sigma, dt);
422  }
423 
424  dX_dt = 0.0;
425 
426  // The big BlockVector stores the fields as follows:
427  // Temperature
428  // Temperature Flux
429  // P field
430  // E field
431  // B field
432  // Joule Heating
433 
434  int Vsize_l2 = L2FESpace.GetVSize();
435  int Vsize_nd = HCurlFESpace.GetVSize();
436  int Vsize_rt = HDivFESpace.GetVSize();
437  int Vsize_h1 = HGradFESpace.GetVSize();
438 
439  Array<int> true_offset(7);
440  true_offset[0] = 0;
441  true_offset[1] = true_offset[0] + Vsize_l2;
442  true_offset[2] = true_offset[1] + Vsize_rt;
443  true_offset[3] = true_offset[2] + Vsize_h1;
444  true_offset[4] = true_offset[3] + Vsize_nd;
445  true_offset[5] = true_offset[4] + Vsize_rt;
446  true_offset[6] = true_offset[5] + Vsize_l2;
447 
448  Vector* xptr = (Vector*) &X;
449  ParGridFunction E, B, T, F, W, P;
450  T.MakeRef(&L2FESpace, *xptr,true_offset[0]);
451  F.MakeRef(&HDivFESpace, *xptr,true_offset[1]);
452  P.MakeRef(&HGradFESpace,*xptr,true_offset[2]);
453  E.MakeRef(&HCurlFESpace,*xptr,true_offset[3]);
454  B.MakeRef(&HDivFESpace, *xptr,true_offset[4]);
455  W.MakeRef(&L2FESpace, *xptr,true_offset[5]);
456 
457  ParGridFunction dE, dB, dT, dF, dW, dP;
458  dT.MakeRef(&L2FESpace, dX_dt,true_offset[0]);
459  dF.MakeRef(&HDivFESpace, dX_dt,true_offset[1]);
460  dP.MakeRef(&HGradFESpace,dX_dt,true_offset[2]);
461  dE.MakeRef(&HCurlFESpace,dX_dt,true_offset[3]);
462  dB.MakeRef(&HDivFESpace, dX_dt,true_offset[4]);
463  dW.MakeRef(&L2FESpace, dX_dt,true_offset[5]);
464 
465  // form the Laplacian and solve it
466  ParGridFunction Phi_gf(&HGradFESpace);
467 
468  // p_bc is given function defining electrostatic potential on surface
469  FunctionCoefficient voltage(p_bc);
470  voltage.SetTime(this->GetTime());
471  Phi_gf = 0.0;
472 
473  // the function below is currently not fully supported on AMR meshes
474  // Phi_gf.ProjectBdrCoefficient(voltage,poisson_ess_bdr);
475 
476  // this is a hack to get around the above issue
477  {
478  Array<int> my_ess_dof_list;
481  homer.ProjectCoefficient(voltage);
482  for (int i = 0; i < my_ess_dof_list.Size(); i++)
483  {
484  int dof = my_ess_dof_list[i];
485  if (dof < 0) { dof = -1 - dof; }
486  Phi_gf[i] = homer[i];
487  }
488  }
489  // end of hack
490 
491  // apply essential BC's and apply static condensation, the new system to
492  // solve is A0 X0 = B0
493  Array<int> poisson_ess_tdof_list;
494  HGradFESpace.GetEssentialTrueDofs(poisson_ess_bdr, poisson_ess_tdof_list);
495 
496  *v0 = 0.0;
497  a0->FormLinearSystem(poisson_ess_tdof_list,Phi_gf,*v0,*A0,*X0,*B0);
498 
499  if (amg_a0 == NULL) { amg_a0 = new HypreBoomerAMG(*A0); }
500  if (pcg_a0 == NULL)
501  {
502  pcg_a0 = new HyprePCG(*A0);
503  pcg_a0->SetTol(SOLVER_TOL);
504  pcg_a0->SetMaxIter(SOLVER_MAX_IT);
505  pcg_a0->SetPrintLevel(SOLVER_PRINT_LEVEL);
507  }
508  // pcg "Mult" operation is a solve
509  // X0 = A0^-1 * B0
510  pcg_a0->Mult(*B0, *X0);
511 
512  // "undo" the static condensation saving result in grid function dP
513  a0->RecoverFEMSolution(*X0,*v0,P);
514  dP = 0.0;
515 
516  // v1 = <1/mu v, curl u> B
517  // B is a grid function but weakCurl is not parallel assembled so is OK
518  weakCurl->MultTranspose(B, *v1);
519 
520  // now add Grad dPhi/dt term
521  // use E as a temporary, E = Grad P
522  // v1 = curl 1/mu B + M1 * Grad P
523  grad->Mult(P,E);
524  m1->AddMult(E,*v1,1.0);
525 
527 
528  // edot_bc is time-derivative E-field on a boundary surface
529  // and then it is used as a Dirichlet BC
530  // the vector v1 will be modified by the values Jtmp and
531  // the part of the matrix m1 that hs been eliminated (but stored).
533  J_gf = 0.0;
535 
536  // form the linear system, including eliminating essential BC's and applying
537  // static condensation. The system to solve is A1 X1 = B1
538  Array<int> ess_tdof_list;
539  HCurlFESpace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
540 
541  a1->FormLinearSystem(ess_tdof_list,J_gf,*v1,*A1,*X1,*B1);
542 
543  // We only need to create the solver and preconditioner once
544  if ( ams_a1 == NULL )
545  {
546  ParFiniteElementSpace *prec_fespace =
548  ams_a1 = new HypreAMS(*A1, prec_fespace);
549  }
550  if ( pcg_a1 == NULL )
551  {
552  pcg_a1 = new HyprePCG(*A1);
553  pcg_a1->SetTol(SOLVER_TOL);
554  pcg_a1->SetMaxIter(SOLVER_MAX_IT);
555  pcg_a1->SetPrintLevel(SOLVER_PRINT_LEVEL);
557  }
558  // solve the system
559  // dE = (A1)^-1 [-S1 E]
560  pcg_a1->Mult(*B1, *X1);
561 
562  // this is required because of static condensation, E is a grid function
563  a1->RecoverFEMSolution(*X1,*v1,E);
564  dE = 0.0;
565 
566  // the total field is E_tot = E_ind - Grad Phi
567  // so we need to subtract out Grad Phi
568  // E = E - grad (P)
569  // note grad maps GF to GF
570  grad->AddMult(P,E,-1.0);
571 
572  // Compute dB/dt = -Curl(E_{n+1})
573  // note curl maps GF to GF
574  curl->Mult(E, dB);
575  dB *= -1.0;
576 
577  // Compute Energy Deposition
578  this->GetJouleHeating(E,W);
579 
580  // v2 = Div^T * W, where W is the Joule heating computed above, and
581  // Div is the matrix <div u, v>
582  weakDivC->MultTranspose(W, *v2);
583  *v2 *= dt;
584 
585  // v2 = <v, div u> T + (1.0)*v2
586  weakDiv->AddMultTranspose(T, *v2, 1.0);
587 
588  // apply the thermal BC
589  Vector zero_vec(3); zero_vec = 0.0;
590  VectorConstantCoefficient Zero_vec(zero_vec);
592  F_gf = 0.0;
594 
595  // form the linear system, including eliminating essential BC's and applying
596  // static condensation. The system to solve is A2 X2 = B2
597  Array<int> thermal_ess_tdof_list;
598  HDivFESpace.GetEssentialTrueDofs(thermal_ess_bdr, thermal_ess_tdof_list);
599  a2->FormLinearSystem(thermal_ess_tdof_list,F_gf,*v2,*A2,*X2,*B2);
600 
601  // We only need to create the solver and preconditioner once
602  if ( ads_a2 == NULL )
603  {
604  ParFiniteElementSpace *prec_fespace =
606  ads_a2 = new HypreADS(*A2, prec_fespace);
607  }
608  if ( pcg_a2 == NULL )
609  {
610  pcg_a2 = new HyprePCG(*A2);
611  pcg_a2->SetTol(SOLVER_TOL);
612  pcg_a2->SetMaxIter(SOLVER_MAX_IT);
613  pcg_a2->SetPrintLevel(SOLVER_PRINT_LEVEL);
615  }
616  // solve for dF from a2 dF = v2
617  // dF = (A2)^-1 [S2*F + rhs]
618  pcg_a2->Mult(*B2, *X2);
619 
620  // this is required because of static condensation
621  a2->RecoverFEMSolution(*X2,*v2,F);
622 
623  // c dT = [W - div F]
624  //
625  // <u,u> dT = <1/c W,u> - <1/c div v,u>
626  //
627  // where W is Joule heating and F is the flux that we just computed
628  //
629  // note: if div is a BilinearForm, then W should be converted to a LoadVector
630  // compute load vector <1/c W, u> where W is the Joule heating GF
631 
632  // create the Coefficient 1/c W
633  //ScaledGFCoefficient Wcoeff(&W, *InvTcap);
634  GridFunctionCoefficient Wcoeff(&W);
635 
636  // compute <W,u>
637  ParLinearForm temp_lf(&L2FESpace);
638  temp_lf.AddDomainIntegrator(new DomainLFIntegrator(Wcoeff));
639  temp_lf.Assemble();
640 
641  // lf = lf - div F
642  weakDiv->AddMult(F, temp_lf, -1.0);
643 
644  // need to perform mass matrix solve to get temperature T
645  // <c u, u> Tdot = -<div v, u> F + <1/c W, u>
646  // NOTE: supposedly we can just invert any L2 matrix, could do that here
647  // instead of a solve
648 
649  if (dsp_m3 == NULL) { dsp_m3 = new HypreDiagScale(*M3); }
650  if (pcg_m3 == NULL)
651  {
652  pcg_m3 = new HyprePCG(*M3);
653  pcg_m3->SetTol(SOLVER_TOL);
654  pcg_m3->SetMaxIter(SOLVER_MAX_IT);
655  pcg_m3->SetPrintLevel(SOLVER_PRINT_LEVEL);
657  }
658 
659  // solve for dT from M3 dT = lf
660  // no boundary conditions on this solve
661  pcg_m3->Mult(temp_lf, dT);
662 }
663 
665 {
666  if ( a0 != NULL ) { delete a0; }
667 
668  // First create and assemble the bilinear form. For now we assume the mesh
669  // isn't moving, the materials are time independent, and dt is constant. So
670  // we only need to do this once.
671 
672  // ConstantCoefficient Sigma(sigma);
675  if (STATIC_COND == 1) { a0->EnableStaticCondensation(); }
676  a0->Assemble();
677 
678  // Don't finalize or parallel assemble this is done in FormLinearSystem.
679 }
680 
683  double dt)
684 {
685  if ( a1 != NULL ) { delete a1; }
686 
687  // First create and assemble the bilinear form. For now we assume the mesh
688  // isn't moving, the materials are time independent, and dt is constant. So
689  // we only need to do this once.
690 
691  ConstantCoefficient dtMuInv(dt*muInv);
695  if (STATIC_COND == 1) { a1->EnableStaticCondensation(); }
696  a1->Assemble();
697 
698  // Don't finalize or parallel assemble this is done in FormLinearSystem.
699 
700  dt_A1 = dt;
701 }
702 
704  MeshDependentCoefficient &InvTcap,
705  double dt)
706 {
707  if ( a2 != NULL ) { delete a2; }
708 
709  InvTcap.SetScaleFactor(dt);
712  a2->AddDomainIntegrator(new DivDivIntegrator(InvTcap));
713  if (STATIC_COND == 1) { a2->EnableStaticCondensation(); }
714  a2->Assemble();
715 
716  // Don't finalize or parallel assemble this is done in FormLinearSystem.
717 
718  dt_A2 = dt;
719 }
720 
722 {
723  if ( m1 != NULL ) { delete m1; }
724 
727  m1->Assemble();
728 
729  // Don't finalize or parallel assemble this is done in FormLinearSystem.
730 }
731 
733 {
734  if ( m2 != NULL ) { delete m2; }
735 
736  // ConstantCoefficient MuInv(muInv);
739  m2->Assemble();
740 
741  // Don't finalize or parallel assemble this is done in FormLinearSystem.
742 }
743 
745 {
746  if ( m3 != NULL ) { delete m3; }
747 
748  // ConstantCoefficient Sigma(sigma);
749  m3 = new ParBilinearForm(&L2FESpace);
750  m3->AddDomainIntegrator(new MassIntegrator(Tcapacity));
751  m3->Assemble();
752  m3->Finalize();
753  M3 = m3->ParallelAssemble();
754 }
755 
757 {
758  if ( s1 != NULL ) { delete s1; }
759 
760  ConstantCoefficient MuInv(muInv);
763  s1->Assemble();
764 }
765 
767 {
768  if ( s2 != NULL ) { delete s2; }
769 
770  // ConstantCoefficient param(a);
772  s2->AddDomainIntegrator(new DivDivIntegrator(InvTcap));
773  s2->Assemble();
774 }
775 
777 {
778  if ( curl != NULL ) { delete curl; }
779  if ( weakCurl != NULL ) { delete weakCurl; }
780 
783  curl->Assemble();
784 
785  ConstantCoefficient MuInv(muInv);
788  weakCurl->Assemble();
789 
790  // no ParallelAssemble since this will be applied to GridFunctions
791 }
792 
794 {
795  if ( weakDiv != NULL ) { delete weakDiv; }
796  if ( weakDivC != NULL ) { delete weakDivC; }
797 
800  weakDivC->Assemble();
801 
804  weakDiv->Assemble();
805 
806  // no ParallelAssemble since this will be applied to GridFunctions
807 }
808 
810 {
811  if ( grad != NULL ) { delete grad; }
812 
815  grad->Assemble();
816 
817  // no ParallelAssemble since this will be applied to GridFunctions
818 }
819 
821 {
822  double el = m1->InnerProduct(E_gf,E_gf);
823 
824  double global_el;
825  MPI_Allreduce(&el, &global_el, 1, MPI_DOUBLE, MPI_SUM,
826  m2->ParFESpace()->GetComm());
827 
828  return el;
829 }
830 
831 // E is the input GF, w is the output GF which is assumed to be an L2 scalar
832 // representing the Joule heating
834  ParGridFunction &w_gf) const
835 {
836  // The w_coeff object stashes a reference to sigma and E, and it has
837  // an Eval method that will be used by ProjectCoefficient.
838  JouleHeatingCoefficient w_coeff(*sigma, E_gf);
839 
840  // This applies the definition of the finite element degrees-of-freedom
841  // to convert the function to a set of discrete values
842  w_gf.ProjectCoefficient(w_coeff);
843 }
844 
846 { t = _t; }
847 
849 {
850  if ( ams_a1 != NULL ) { delete ams_a1; }
851  if ( pcg_a1 != NULL ) { delete pcg_a1; }
852 
853  if ( dsp_m1 != NULL ) { delete dsp_m1; }
854  if ( pcg_m1 != NULL ) { delete pcg_m1; }
855 
856  if ( dsp_m2 != NULL ) { delete dsp_m2; }
857  if ( pcg_m2 != NULL ) { delete pcg_m2; }
858 
859  if ( curl != NULL ) { delete curl; }
860  if ( weakDiv != NULL ) { delete weakDiv; }
861  if ( weakDivC != NULL ) { delete weakDivC; }
862  if ( weakCurl != NULL ) { delete weakCurl; }
863  if ( grad != NULL ) { delete grad; }
864 
865  if ( a0 != NULL ) { delete a0; }
866  if ( a1 != NULL ) { delete a1; }
867  if ( a2 != NULL ) { delete a2; }
868  if ( m1 != NULL ) { delete m1; }
869  if ( m2 != NULL ) { delete m2; }
870  if ( s1 != NULL ) { delete s1; }
871  if ( s2 != NULL ) { delete s2; }
872 
873  if ( A0 != NULL ) { delete A0; }
874  if ( X0 != NULL ) { delete X0; }
875  if ( B0 != NULL ) { delete B0; }
876 
877  if ( A1 != NULL ) { delete A1; }
878  if ( X1 != NULL ) { delete X1; }
879  if ( B1 != NULL ) { delete B1; }
880 
881  if ( A2 != NULL ) { delete A2; }
882  if ( X2 != NULL ) { delete X2; }
883  if ( B2 != NULL ) { delete B2; }
884 
885  if ( v1 != NULL ) { delete v1; }
886  if ( v2 != NULL ) { delete v2; }
887 
888  if (sigma != NULL) { delete sigma; }
889  if (Tcapacity != NULL) { delete Tcapacity; }
890  if (InvTcap != NULL) { delete InvTcap; }
891  if (InvTcond != NULL) { delete InvTcond; }
892 }
893 
894 void MagneticDiffusionEOperator::Debug(const char *base, double)
895 {
896  {
897  hypre_ParCSRMatrixPrint(*A1,"A1_");
898  HypreParVector tempB1(A1->GetComm(),A1->N(),B1->GetData(),A1->ColPart());
899  tempB1.Print("B1_");
900  HypreParVector tempX1(A1->GetComm(),A1->N(),X1->GetData(),A1->ColPart());
901  tempX1.Print("X1_");
902  }
903 
904  {
905  hypre_ParCSRMatrixPrint(*A2,"A2_");
906  HypreParVector tempB2(A2->GetComm(),A2->N(),B2->GetData(),A2->ColPart());
907  tempB2.Print("B2_");
908  HypreParVector tempX2(A2->GetComm(),A2->N(),X2->GetData(),A2->ColPart());
909  tempX2.Print("X2_");
910  }
911 }
912 
914  const IntegrationPoint &ip)
915 {
916  Vector E;
917  double thisSigma;
918  E_gf.GetVectorValue(T.ElementNo, ip, E);
919  thisSigma = sigma.Eval(T, ip);
920  return thisSigma*(E*E);
921 }
922 
924  const std::map<int, double> &inputMap, double scale)
925  : Coefficient()
926 {
927  // make a copy of the magic attribute-value map for later use
928  materialMap = new std::map<int, double>(inputMap);
929  scaleFactor = scale;
930 }
931 
933  const MeshDependentCoefficient &cloneMe)
934  : Coefficient()
935 {
936  // make a copy of the magic attribute-value map for later use
937  materialMap = new std::map<int, double>(*(cloneMe.materialMap));
938  scaleFactor = cloneMe.scaleFactor;
939 }
940 
942  const IntegrationPoint &ip)
943 {
944  // given the attribute, extract the coefficient value from the map
945  std::map<int, double>::iterator it;
946  int thisAtt = T.Attribute;
947  double value;
948  it = materialMap->find(thisAtt);
949  if (it != materialMap->end())
950  {
951  value = it->second;
952  }
953  else
954  {
955  value = 0.0; // avoid compile warning
956  std::cerr << "MeshDependentCoefficient attribute " << thisAtt
957  << " not found" << std::endl;
958  mfem_error();
959  }
960 
961  return value*scaleFactor;
962 }
963 
965  MeshDependentCoefficient &input_mdc)
966  : GridFunctionCoefficient(gf), mdc(input_mdc) {}
967 
969  const IntegrationPoint &ip)
970 {
971  return mdc.Eval(T,ip) * GridFunctionCoefficient::Eval(T,ip);
972 }
973 
974 } // namespace electromagnetics
975 
976 } // namespace mfem
977 
978 #endif // MFEM_USE_MPI
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
void buildDiv(MeshDependentCoefficient &InvTcap)
void SetTol(double tol)
Definition: hypre.cpp:2088
virtual double GetTime() const
Read the currently set time.
Definition: operator.hpp:169
int Size() const
Logical size of the array.
Definition: array.hpp:110
Class for domain integration L(v) := (f, v)
Definition: lininteg.hpp:93
int GetVSize() const
Definition: fespace.hpp:163
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:529
MPI_Comm GetComm() const
MPI communicator.
Definition: hypre.hpp:298
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:846
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:27
The Auxiliary-space Divergence Solver in hypre.
Definition: hypre.hpp:878
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_dofs, int component=-1) const
Determine the boundary degrees of freedom.
Definition: pfespace.cpp:515
Subclass constant coefficient.
Definition: coefficient.hpp:57
HYPRE_Int N() const
Returns the global number of columns.
Definition: hypre.hpp:342
void FormLinearSystem(const Array< int > &ess_tdof_list, Vector &x, Vector &b, OperatorHandle &A, Vector &X, Vector &B, int copy_interior=0)
void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:370
HYPRE_Int * ColPart()
Returns the column partitioning.
Definition: hypre.hpp:334
virtual void AddMult(const Vector &x, Vector &y, const double a=1.0) const
Integrator for (curl u, curl v) for Nedelec elements.
Base abstract class for time dependent operators.
Definition: operator.hpp:142
void Assemble()
Assembles the linear form i.e. sums over all domain/bdr integrators.
Definition: linearform.cpp:51
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:31
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: pgridfunc.cpp:275
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
double * GetData() const
Definition: vector.hpp:121
(Q div u, div v) for RT elements
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:2103
void AddDomainInterpolator(DiscreteInterpolator *di)
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Debug(const char *basefilename, double time)
double t
Current time.
Definition: operator.hpp:153
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:169
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Definition: pgridfunc.cpp:86
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:765
The BoomerAMG solver in hypre.
Definition: hypre.hpp:783
Class for parallel linear form.
Definition: plinearform.hpp:26
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Definition: gridfunc.cpp:1549
MPI_Comm GetComm() const
Definition: pfespace.hpp:166
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:744
ParFiniteElementSpace * ParFESpace() const
Return the parallel FE space associated with the ParBilinearForm.
void buildM1(MeshDependentCoefficient &sigma)
void Assemble(int skip_zeros=1)
void Assemble(int skip_zeros=1)
Assemble the local matrix.
void SetTime(double t)
Definition: coefficient.hpp:39
void SetMaxIter(int max_iter)
Definition: hypre.cpp:2093
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:73
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator.
Definition: linearform.cpp:19
double ElectricLosses(ParGridFunction &E_gf) const
bool StaticCondensationIsEnabled() const
MeshDependentCoefficient(const std::map< int, double > &inputMap, double scale=1.0)
PCG solver in hypre.
Definition: hypre.hpp:643
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
void mfem_error(const char *msg)
Definition: error.cpp:107
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:503
double InnerProduct(const Vector &x, const Vector &y) const
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)
Definition: coefficient.cpp:49
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.
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.
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 ...
void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Definition: gridfunc.cpp:1620
Class for parallel bilinear form.
void edot_bc(const Vector &x, Vector &E)
Definition: joule.cpp:740
class for C-function coefficient
Vector data type.
Definition: vector.hpp:41
Class for parallel grid function.
Definition: pgridfunc.hpp:32
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
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
void buildA0(MeshDependentCoefficient &sigma)
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
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...
Integrator for (Q u, v) for VectorFiniteElements.
void EnableStaticCondensation()
void SetTime(const double _t)
Set the current time.