MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
seq_test.cpp
Go to the documentation of this file.
1// Copyright (c) 2010-2024, 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
15template<typename TDataType, typename TParamVector, typename TStateVector
16 , int state_size, int param_size>
17class DiffusionFunctional
18{
19public:
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
32template<typename TDataType, typename TParamVector, typename TStateVector,
33 int residual_size, int state_size, int param_size>
34class DiffusionResidual
35{
36public:
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
53int 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 mfem::real_t kappa = param[0];
134 mfem::real_t 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
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}
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
virtual void Print(std::ostream &out=mfem::out, int width_=4) const
Prints matrix to stream out.
void Hessian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
Definition admfem.hpp:345
real_t Eval(const mfem::Vector &vparam, mfem::Vector &uu)
Definition admfem.hpp:248
void Grad(const mfem::Vector &vparam, mfem::Vector &uu, mfem::Vector &rr)
Definition admfem.hpp:261
void Jacobian(mfem::Vector &vparam, mfem::Vector &uu, mfem::DenseMatrix &jac)
Definition admfem.hpp:159
Templated vector data type.
Definition tadvector.hpp:41
void Jacobian(mfem::Vector &vparam, mfem::Vector &vstate, mfem::DenseMatrix &jac)
Definition admfem.hpp:70
Vector data type.
Definition vector.hpp:80
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition vector.cpp:755
real_t kappa
Definition ex24.cpp:54
int main()
float real_t
Definition config.hpp:43