MFEM  v4.5.2
Finite element discretization library
tmop_pa_h2s.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, 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 "../tmop.hpp"
13 #include "tmop_pa.hpp"
14 #include "../../general/forall.hpp"
15 #include "../../linalg/kernels.hpp"
16 #include "../../linalg/dinvariants.hpp"
17 
18 namespace mfem
19 {
20 
22 
23 // weight * ddI1
24 static MFEM_HOST_DEVICE inline
25 void EvalH_001(const int e, const int qx, const int qy,
26  const double weight, const double *Jpt,
28 {
29  constexpr int DIM = 2;
30  double ddI1[4];
31  kernels::InvariantsEvaluator2D ie(Args().J(Jpt).ddI1(ddI1));
32  for (int i = 0; i < DIM; i++)
33  {
34  for (int j = 0; j < DIM; j++)
35  {
36  ConstDeviceMatrix ddi1(ie.Get_ddI1(i,j),DIM,DIM);
37  for (int r = 0; r < DIM; r++)
38  {
39  for (int c = 0; c < DIM; c++)
40  {
41  const double h = ddi1(r,c);
42  H(r,c,i,j,qx,qy,e) = weight * h;
43  }
44  }
45  }
46  }
47 }
48 
49 // 0.5 * weight * dI1b
50 static MFEM_HOST_DEVICE inline
51 void EvalH_002(const int e, const int qx, const int qy,
52  const double weight, const double *Jpt,
53  DeviceTensor<7,double> H)
54 {
55  constexpr int DIM = 2;
56  double ddI1[4], ddI1b[4], dI2b[4];
57  kernels::InvariantsEvaluator2D ie(Args()
58  .J(Jpt)
59  .ddI1(ddI1)
60  .ddI1b(ddI1b)
61  .dI2b(dI2b));
62  const double w = 0.5 * weight;
63  for (int i = 0; i < DIM; i++)
64  {
65  for (int j = 0; j < DIM; j++)
66  {
67  ConstDeviceMatrix ddi1b(ie.Get_ddI1b(i,j),DIM,DIM);
68  for (int r = 0; r < DIM; r++)
69  {
70  for (int c = 0; c < DIM; c++)
71  {
72  const double h = ddi1b(r,c);
73  H(r,c,i,j,qx,qy,e) = w * h;
74  }
75  }
76  }
77  }
78 }
79 
80 static MFEM_HOST_DEVICE inline
81 void EvalH_007(const int e, const int qx, const int qy,
82  const double weight, const double *Jpt,
83  DeviceTensor<7,double> H)
84 {
85  constexpr int DIM = 2;
86  double ddI1[4], ddI2[4], dI1[4], dI2[4], dI2b[4];
87  kernels::InvariantsEvaluator2D ie(Args()
88  .J(Jpt)
89  .ddI1(ddI1)
90  .ddI2(ddI2)
91  .dI1(dI1)
92  .dI2(dI2)
93  .dI2b(dI2b));
94  const double c1 = 1./ie.Get_I2();
95  const double c2 = weight*c1*c1;
96  const double c3 = ie.Get_I1()*c2;
97  ConstDeviceMatrix di1(ie.Get_dI1(),DIM,DIM);
98  ConstDeviceMatrix di2(ie.Get_dI2(),DIM,DIM);
99 
100  for (int i = 0; i < DIM; i++)
101  {
102  for (int j = 0; j < DIM; j++)
103  {
104  ConstDeviceMatrix ddi1(ie.Get_ddI1(i,j),DIM,DIM);
105  ConstDeviceMatrix ddi2(ie.Get_ddI2(i,j),DIM,DIM);
106  for (int r = 0; r < DIM; r++)
107  {
108  for (int c = 0; c < DIM; c++)
109  {
110  H(r,c,i,j,qx,qy,e) =
111  weight * (1.0 + c1) * ddi1(r,c)
112  - c3 * ddi2(r,c)
113  - c2 * ( di1(i,j) * di2(r,c) + di2(i,j) * di1(r,c) )
114  + 2.0 * c1 * c3 * di2(r,c) * di2(i,j);
115  }
116  }
117  }
118  }
119 }
120 
121 static MFEM_HOST_DEVICE inline
122 void EvalH_077(const int e, const int qx, const int qy,
123  const double weight, const double *Jpt,
124  DeviceTensor<7,double> H)
125 {
126  constexpr int DIM = 2;
127  double dI2[4], dI2b[4], ddI2[4];
128  kernels::InvariantsEvaluator2D ie(Args()
129  .J(Jpt)
130  .dI2(dI2)
131  .dI2b(dI2b)
132  .ddI2(ddI2));
133  const double I2 = ie.Get_I2(), I2inv_sq = 1.0 / (I2 * I2);
134  ConstDeviceMatrix di2(ie.Get_dI2(),DIM,DIM);
135  for (int i = 0; i < DIM; i++)
136  {
137  for (int j = 0; j < DIM; j++)
138  {
139  ConstDeviceMatrix ddi2(ie.Get_ddI2(i,j),DIM,DIM);
140  for (int r = 0; r < DIM; r++)
141  {
142  for (int c = 0; c < DIM; c++)
143  {
144  H(r,c,i,j,qx,qy,e) =
145  weight * 0.5 * (1.0 - I2inv_sq) * ddi2(r,c)
146  + weight * (I2inv_sq / I2) * di2(r,c) * di2(i,j);
147  }
148  }
149  }
150  }
151 }
152 
153 static MFEM_HOST_DEVICE inline
154 void EvalH_080(const int e, const int qx, const int qy,
155  const double weight, const double gamma, const double *Jpt,
156  DeviceTensor<7,double> H)
157 {
158  // h_80 = (1-gamma) h_2 + gamma h_77.
159 
160  constexpr int DIM = 2;
161  double ddI1[4], ddI1b[4], dI2[4], dI2b[4], ddI2[4];
162  kernels::InvariantsEvaluator2D ie(Args()
163  .J(Jpt)
164  .dI2(dI2)
165  .ddI1(ddI1)
166  .ddI1b(ddI1b)
167  .dI2b(dI2b)
168  .ddI2(ddI2));
169 
170  const double I2 = ie.Get_I2(), I2inv_sq = 1.0 / (I2 * I2);
171  ConstDeviceMatrix di2(ie.Get_dI2(),DIM,DIM);
172  for (int i = 0; i < DIM; i++)
173  {
174  for (int j = 0; j < DIM; j++)
175  {
176  ConstDeviceMatrix ddi1b(ie.Get_ddI1b(i,j),DIM,DIM);
177  ConstDeviceMatrix ddi2(ie.Get_ddI2(i,j),DIM,DIM);
178  for (int r = 0; r < DIM; r++)
179  {
180  for (int c = 0; c < DIM; c++)
181  {
182  H(r,c,i,j,qx,qy,e) =
183  (1.0 - gamma) * 0.5 * weight * ddi1b(r,c) +
184  gamma * ( weight * 0.5 * (1.0 - I2inv_sq) * ddi2(r,c) +
185  weight * (I2inv_sq / I2) * di2(r,c) * di2(i,j) );
186  }
187  }
188  }
189  }
190 }
191 
192 MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_2D,
193  const Vector &x_,
194  const double metric_normal,
195  const double metric_param,
196  const int mid,
197  const int NE,
198  const Array<double> &w_,
199  const Array<double> &b_,
200  const Array<double> &g_,
201  const DenseTensor &j_,
202  Vector &h_,
203  const int d1d,
204  const int q1d)
205 {
206  MFEM_VERIFY(mid == 1 || mid == 2 || mid == 7 || mid == 77 || mid == 80,
207  "Metric not yet implemented!");
208 
209  constexpr int DIM = 2;
210  constexpr int NBZ = 1;
211  const int D1D = T_D1D ? T_D1D : d1d;
212  const int Q1D = T_Q1D ? T_Q1D : q1d;
213 
214  const auto W = Reshape(w_.Read(), Q1D, Q1D);
215  const auto b = Reshape(b_.Read(), Q1D, D1D);
216  const auto g = Reshape(g_.Read(), Q1D, D1D);
217  const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, NE);
218  const auto X = Reshape(x_.Read(), D1D, D1D, DIM, NE);
219  auto H = Reshape(h_.Write(), DIM, DIM, DIM, DIM, Q1D, Q1D, NE);
220 
221  MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
222  {
223  const int D1D = T_D1D ? T_D1D : d1d;
224  const int Q1D = T_Q1D ? T_Q1D : q1d;
225  constexpr int NBZ = 1;
226  constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
227  constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
228 
229  MFEM_SHARED double s_BG[2][MQ1*MD1];
230  MFEM_SHARED double s_X[2][NBZ][MD1*MD1];
231  MFEM_SHARED double s_DQ[4][NBZ][MD1*MQ1];
232  MFEM_SHARED double s_QQ[4][NBZ][MQ1*MQ1];
233 
234  kernels::internal::LoadX<MD1,NBZ>(e,D1D,X,s_X);
235  kernels::internal::LoadBG<MD1,MQ1>(D1D, Q1D, b, g, s_BG);
236 
237  kernels::internal::GradX<MD1,MQ1,NBZ>(D1D, Q1D, s_BG, s_X, s_DQ);
238  kernels::internal::GradY<MD1,MQ1,NBZ>(D1D, Q1D, s_BG, s_DQ, s_QQ);
239 
240  MFEM_FOREACH_THREAD(qy,y,Q1D)
241  {
242  MFEM_FOREACH_THREAD(qx,x,Q1D)
243  {
244  const double *Jtr = &J(0,0,qx,qy,e);
245  const double detJtr = kernels::Det<2>(Jtr);
246  const double weight = metric_normal * W(qx,qy) * detJtr;
247 
248  // Jrt = Jtr^{-1}
249  double Jrt[4];
250  kernels::CalcInverse<2>(Jtr, Jrt);
251 
252  // Jpr = X^t.DSh
253  double Jpr[4];
254  kernels::internal::PullGrad<MQ1,NBZ>(Q1D,qx,qy,s_QQ,Jpr);
255 
256  // Jpt = Jpr.Jrt
257  double Jpt[4];
258  kernels::Mult(2,2,2, Jpr, Jrt, Jpt);
259 
260  // metric->AssembleH
261  if (mid == 1) { EvalH_001(e,qx,qy,weight,Jpt,H); }
262  if (mid == 2) { EvalH_002(e,qx,qy,weight,Jpt,H); }
263  if (mid == 7) { EvalH_007(e,qx,qy,weight,Jpt,H); }
264  if (mid == 77) { EvalH_077(e,qx,qy,weight,Jpt,H); }
265  if (mid == 80) { EvalH_080(e,qx,qy,weight,metric_param,Jpt,H); }
266  } // qx
267  } // qy
268  });
269 }
270 
272 {
273  const int N = PA.ne;
274  const int M = metric->Id();
275  const int D1D = PA.maps->ndof;
276  const int Q1D = PA.maps->nqpt;
277  const int id = (D1D << 4 ) | Q1D;
278  const double mn = metric_normal;
279  const DenseTensor &J = PA.Jtr;
280  const Array<double> &W = PA.ir->GetWeights();
281  const Array<double> &B = PA.maps->B;
282  const Array<double> &G = PA.maps->G;
283  Vector &H = PA.H;
284 
285  double mp = 0.0;
286  if (auto m = dynamic_cast<TMOP_Metric_080 *>(metric)) { mp = m->GetGamma(); }
287 
288  MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_2D,id,X,mn,mp,M,N,W,B,G,J,H);
289 }
290 
291 } // namespace mfem
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
DeviceTensor< 2, const double > ConstDeviceMatrix
Definition: dtensor.hpp:144
struct mfem::TMOP_Integrator::@23 PA
void AssembleGradPA_2D(const Vector &) const
TMOP_QualityMetric * metric
Definition: tmop.hpp:1645
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:448
constexpr int DIM
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read( GetMemory(), TotalSize(), on_dev).
Definition: densemat.hpp:1112
virtual int Id() const
Return the metric ID.
Definition: tmop.hpp:78
MFEM_HOST_DEVICE double * Get_ddI1(int i, int j)
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:456
double b
Definition: lissajous.cpp:42
MFEM_REGISTER_TMOP_KERNELS(void, DatcSize, const int NE, const int ncomp, const int sizeidx, const DenseMatrix &w_, const Array< double > &b_, const Vector &x_, DenseTensor &j_, const int d1d, const int q1d)
Definition: tmop_pa_da3.cpp:20
A basic generic Tensor class, appropriate for use on the GPU.
Definition: dtensor.hpp:81
MFEM_HOST_DEVICE void Mult(const int height, const int width, const TA *data, const TX *x, TY *y)
Matrix vector multiplication: y = A x, where the matrix A is of size height x width with given data...
Definition: kernels.hpp:163
kernels::InvariantsEvaluator2D::Buffers Args
Definition: tmop_pa_h2s.cpp:21
Vector data type.
Definition: vector.hpp:60
MFEM_HOST_DEVICE DeviceTensor< sizeof...(Dims), T > Reshape(T *ptr, Dims... dims)
Wrap a pointer as a DeviceTensor with automatically deduced template parameters.
Definition: dtensor.hpp:131
Rank 3 tensor (array of matrices)
Definition: densemat.hpp:978