MFEM  v4.5.2
Finite element discretization library
tmop_pa_p3_c0.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 "../linearform.hpp"
15 #include "../../general/forall.hpp"
16 #include "../../linalg/kernels.hpp"
17 
18 namespace mfem
19 {
20 
21 MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
22  const double lim_normal,
23  const Vector &lim_dist,
24  const Vector &c0_,
25  const int NE,
26  const DenseTensor &j_,
27  const Array<double> &w_,
28  const Array<double> &b_,
29  const Array<double> &bld_,
30  const Vector &x0_,
31  const Vector &x1_,
32  Vector &y_,
33  const bool exp_lim,
34  const int d1d,
35  const int q1d)
36 {
37  const bool const_c0 = c0_.Size() == 1;
38 
39  constexpr int DIM = 3;
40  const int D1D = T_D1D ? T_D1D : d1d;
41  const int Q1D = T_Q1D ? T_Q1D : q1d;
42 
43  const auto C0 = const_c0 ?
44  Reshape(c0_.Read(), 1, 1, 1, 1) :
45  Reshape(c0_.Read(), Q1D, Q1D, Q1D, NE);
46  const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, NE);
47  const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
48  const auto b = Reshape(b_.Read(), Q1D, D1D);
49  const auto bld = Reshape(bld_.Read(), Q1D, D1D);
50  const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
51  const auto X0 = Reshape(x0_.Read(), D1D, D1D, D1D, DIM, NE);
52  const auto X1 = Reshape(x1_.Read(), D1D, D1D, D1D, DIM, NE);
53 
54  auto Y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, DIM, NE);
55 
56  MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
57  {
58  const int D1D = T_D1D ? T_D1D : d1d;
59  const int Q1D = T_Q1D ? T_Q1D : q1d;
60  constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
61  constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
62  constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
63 
64  MFEM_SHARED double B[MQ1*MD1];
65  MFEM_SHARED double sBLD[MQ1*MD1];
66  kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
67  ConstDeviceMatrix BLD(sBLD, D1D, Q1D);
68 
69  MFEM_SHARED double sm0[MDQ*MDQ*MDQ];
70  MFEM_SHARED double sm1[MDQ*MDQ*MDQ];
71  DeviceCube DDD(sm0, MD1,MD1,MD1);
72  DeviceCube DDQ(sm1, MD1,MD1,MQ1);
73  DeviceCube DQQ(sm0, MD1,MQ1,MQ1);
74  DeviceCube QQQ(sm1, MQ1,MQ1,MQ1);
75 
76  MFEM_SHARED double DDD0[3][MD1*MD1*MD1];
77  MFEM_SHARED double DDQ0[3][MD1*MD1*MQ1];
78  MFEM_SHARED double DQQ0[3][MD1*MQ1*MQ1];
79  MFEM_SHARED double QQQ0[3][MQ1*MQ1*MQ1];
80 
81  MFEM_SHARED double DDD1[3][MD1*MD1*MD1];
82  MFEM_SHARED double DDQ1[3][MD1*MD1*MQ1];
83  MFEM_SHARED double DQQ1[3][MD1*MQ1*MQ1];
84  MFEM_SHARED double QQQ1[3][MQ1*MQ1*MQ1];
85 
86  kernels::internal::LoadX(e,D1D,LD,DDD);
87  kernels::internal::LoadX<MD1>(e,D1D,X0,DDD0);
88  kernels::internal::LoadX<MD1>(e,D1D,X1,DDD1);
89 
90  kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
91 
92  kernels::internal::EvalX(D1D,Q1D,BLD,DDD,DDQ);
93  kernels::internal::EvalY(D1D,Q1D,BLD,DDQ,DQQ);
94  kernels::internal::EvalZ(D1D,Q1D,BLD,DQQ,QQQ);
95 
96  kernels::internal::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD0,DDQ0);
97  kernels::internal::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ0,DQQ0);
98  kernels::internal::EvalZ<MD1,MQ1>(D1D,Q1D,B,DQQ0,QQQ0);
99 
100  kernels::internal::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD1,DDQ1);
101  kernels::internal::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ1,DQQ1);
102  kernels::internal::EvalZ<MD1,MQ1>(D1D,Q1D,B,DQQ1,QQQ1);
103 
104  MFEM_FOREACH_THREAD(qz,z,Q1D)
105  {
106  MFEM_FOREACH_THREAD(qy,y,Q1D)
107  {
108  MFEM_FOREACH_THREAD(qx,x,Q1D)
109  {
110  const double *Jtr = &J(0,0,qx,qy,qz,e);
111  const double detJtr = kernels::Det<3>(Jtr);
112  const double weight = W(qx,qy,qz) * detJtr;
113 
114  double D, p0[3], p1[3];
115  const double coeff0 = const_c0 ? C0(0,0,0,0) : C0(qx,qy,qz,e);
116  kernels::internal::PullEval(qx,qy,qz,QQQ,D);
117  kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ0,p0);
118  kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ1,p1);
119 
120  double d1[3];
121  // Eval_d1 (Quadratic Limiter)
122  // subtract(1.0 / (dist * dist), x, x0, d1);
123  // z = a * (x - y)
124  // grad = a * (x - x0)
125 
126  // Eval_d1 (Exponential Limiter)
127  // double dist_squared = dist*dist;
128  // subtract(20.0*exp(10.0*((x.DistanceSquaredTo(x0) / dist_squared) - 1.0)) /
129  // dist_squared, x, x0, d1);
130  // z = a * (x - y)
131  // grad = a * (x - x0)
132  const double dist = D; // GetValues, default comp set to 0
133  double a = 0.0;
134  const double w = weight * lim_normal * coeff0;
135  const double dist_squared = dist * dist;
136 
137  if (!exp_lim)
138  {
139  a = 1.0 / dist_squared;
140  }
141  else
142  {
143  double dsq = kernels::DistanceSquared<3>(p1,p0) / dist_squared;
144  a = 20.0*exp(10.0*(dsq - 1.0))/dist_squared;
145  }
146 
147  kernels::Subtract<3>(w*a, p1, p0, d1);
148  kernels::internal::PushEval<MQ1>(Q1D,qx,qy,qz,d1,QQQ0);
149  }
150  }
151  }
152  MFEM_SYNC_THREAD;
153  kernels::internal::LoadBt<MD1,MQ1>(D1D,Q1D,b,B);
154  kernels::internal::EvalXt<MD1,MQ1>(D1D,Q1D,B,QQQ0,DQQ0);
155  kernels::internal::EvalYt<MD1,MQ1>(D1D,Q1D,B,DQQ0,DDQ0);
156  kernels::internal::EvalZt<MD1,MQ1>(D1D,Q1D,B,DDQ0,Y,e);
157  });
158 }
159 
161 {
162  const int N = PA.ne;
163  const int D1D = PA.maps->ndof;
164  const int Q1D = PA.maps->nqpt;
165  const int id = (D1D << 4 ) | Q1D;
166  const double ln = lim_normal;
167  const Vector &LD = PA.LD;
168  const DenseTensor &J = PA.Jtr;
169  const Array<double> &W = PA.ir->GetWeights();
170  const Array<double> &B = PA.maps->B;
171  const Array<double> &BLD = PA.maps_lim->B;
172  MFEM_VERIFY(PA.maps_lim->ndof == D1D, "");
173  MFEM_VERIFY(PA.maps_lim->nqpt == Q1D, "");
174  const Vector &X0 = PA.X0;
175  const Vector &C0 = PA.C0;
176  auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
177  const bool exp_lim = (el) ? true : false;
178 
179  MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y,
180  exp_lim);
181 }
182 
183 } // namespace mfem
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
struct mfem::TMOP_Integrator::@23 PA
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
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
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
void AddMultPA_C0_3D(const Vector &, Vector &) const
double a
Definition: lissajous.cpp:41
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:464
Exponential limiter function in TMOP_Integrator.
Definition: tmop.hpp:1124
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
TMOP_LimiterFunction * lim_func
Definition: tmop.hpp:1665