MFEM  v4.6.0
Finite element discretization library
seq_test.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 "admfem.hpp"
13 #include "mfem.hpp"
14 
15 template<typename TDataType, typename TParamVector, typename TStateVector
16  , int state_size, int param_size>
17 class DiffusionFunctional
18 {
19 public:
20  TDataType operator() (TParamVector& vparam, TStateVector& uu)
21  {
22  MFEM_ASSERT(state_size==4,"ExampleFunctor state_size should be equal to 4!");
23  MFEM_ASSERT(param_size==2,"ExampleFunctor param_size should be equal to 2!");
24  auto kappa = vparam[0]; // diffusion coefficient
25  auto load = vparam[1]; // volumetric influx
26  TDataType rez = kappa*(uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2])/2.0 - load*uu[3];
27  return rez;
28  }
29 
30 };
31 
32 template<typename TDataType, typename TParamVector, typename TStateVector,
33  int residual_size, int state_size, int param_size>
34 class DiffusionResidual
35 {
36 public:
37  void operator ()(TParamVector& vparam, TStateVector& uu, TStateVector& rr)
38  {
39  MFEM_ASSERT(residual_size==4,
40  "DiffusionResidual residual_size should be equal to 4!");
41  MFEM_ASSERT(state_size==4,"ExampleFunctor state_size should be equal to 4!");
42  MFEM_ASSERT(param_size==2,"ExampleFunctor param_size should be equal to 2!");
43  auto kappa = vparam[0]; // diffusion coefficient
44  auto load = vparam[1]; // volumetric influx
45 
46  rr[0] = kappa * uu[0];
47  rr[1] = kappa * uu[1];
48  rr[2] = kappa * uu[2];
49  rr[3] = -load;
50  }
51 };
52 
53 int main(int argc, char *argv[])
54 {
55 
56 #ifdef MFEM_USE_ADFORWARD
57  std::cout<<"MFEM_USE_ADFORWARD == true"<<std::endl;
58 #else
59  std::cout<<"MFEM_USE_ADFORWARD == false"<<std::endl;
60 #endif
61 
62 #ifdef MFEM_USE_CALIPER
63  cali::ConfigManager mgr;
64 #endif
65  // Caliper instrumentation
66  MFEM_PERF_FUNCTION;
67 #ifdef MFEM_USE_CALIPER
68  const char* cali_config = "runtime-report";
69  mgr.add(cali_config);
70  mgr.start();
71 #endif
72  mfem::Vector param(2);
73  param[0]=3.0; // diffusion coefficient
74  param[1]=2.0; // volumetric influx
75 
76  mfem::Vector state(4);
77  state[0]=1.0; // grad_x
78  state[1]=2.0; // grad_y
79  state[2]=3.0; // grad_z
80  state[3]=4.0; // state value
81 
83 
84  mfem::Vector rr0(4);
85  mfem::DenseMatrix hh0(4,4);
86 
87  mfem::Vector rr1(4);
88  mfem::DenseMatrix hh1(4,4);
89  MFEM_PERF_BEGIN("Grad");
90  adf.Grad(param,state,rr0);
91  MFEM_PERF_END("Grad");
92  MFEM_PERF_BEGIN("Hessian");
93  adf.Hessian(param, state, hh0);
94  MFEM_PERF_END("Hessian");
95  // dump out the results
96  std::cout<<"FunctionAutoDiff"<<std::endl;
97  std::cout<< adf.Eval(param,state)<<std::endl;
98  rr0.Print(std::cout);
99  hh0.Print(std::cout);
100 
102  MFEM_PERF_BEGIN("Jacobian");
103  rdf.Jacobian(param, state, hh1);
104  MFEM_PERF_END("Jacobian");
105 
106  std::cout<<"ResidualAutoDiff"<<std::endl;
107  hh1.Print(std::cout);
108 
109  // using lambda expression
110  auto func = [](mfem::Vector& vparam,
113  {
114  // auto func = [](auto& vparam, auto& uu, auto& vres) { //c++14
115  auto kappa = vparam[0]; // diffusion coefficient
116  auto load = vparam[1]; // volumetric influx
117 
118  vres[0] = kappa * uu[0];
119  vres[1] = kappa * uu[1];
120  vres[2] = kappa * uu[2];
121  vres[3] = -load;
122  };
123 
125  MFEM_PERF_BEGIN("JacobianV");
126  fdr.Jacobian(param,state,
127  hh1); // computes the gradient of func and stores the result in hh1
128  MFEM_PERF_END("JacobianV");
129  std::cout<<"LambdaAutoDiff"<<std::endl;
130  hh1.Print(std::cout);
131 
132 
133  double kappa = param[0];
134  double load = param[1];
135  // using lambda expression
136  auto func01 = [&kappa,&load](mfem::Vector& vparam,
139  {
140  // auto func = [](auto& vparam, auto& uu, auto& vres) { //c++14
141 
142  vres[0] = kappa * uu[0];
143  vres[1] = kappa * uu[1];
144  vres[2] = kappa * uu[2];
145  vres[3] = -load;
146  };
147 
148  mfem::VectorFuncAutoDiff<4,4,2> fdr01(func01);
149  MFEM_PERF_BEGIN("Jacobian1");
150  fdr01.Jacobian(param,state,hh1);
151  MFEM_PERF_END("Jacobian1");
152  std::cout<<"LambdaAutoDiff 01"<<std::endl;
153  hh1.Print(std::cout);
154 
155 #ifdef MFEM_USE_CALIPER
156  mgr.flush();
157 #endif
158 }
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:756
void Grad(const mfem::Vector &vparam, mfem::Vector &uu, mfem::Vector &rr)
Definition: admfem.hpp:261
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
double kappa
Definition: ex24.cpp:54
virtual void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
Definition: densemat.cpp:2221
void Jacobian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
Definition: admfem.hpp:159
int main(int argc, char *argv[])
Definition: seq_test.cpp:53
Templated vector data type.
Definition: tadvector.hpp:40
Vector data type.
Definition: vector.hpp:58
void Hessian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
Definition: admfem.hpp:345
void Jacobian(mfem::Vector &vparam, mfem::Vector &vstate, mfem::DenseMatrix &jac)
Definition: admfem.hpp:70
double Eval(const mfem::Vector &vparam, mfem::Vector &uu)
Definition: admfem.hpp:248