MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
adjoint_advection_diffusion.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// ----------------------------------------------------------
13// Advection-Diffusion Miniapp: Parallel MFEM CVODES Example
14// ----------------------------------------------------------
15//
16// Compile with: make adjoint_advection_diffusion
17//
18// Sample runs: adjoint_advection_diffusion -dt 0.01 -tf 2.5
19// adjoint_advection_diffusion -dt 0.005
20//
21// Description: This example is a port of cvodes/parallel/cvsAdvDiff_ASAp_non_p
22// example that is part of SUNDIALS. The goal is to demonstrate
23// how to use the adjoint SUNDIALS CVODES interface with MFEM.
24// Below is an excerpt description from the aforementioned file.
25//
26// Example problem:
27//
28// The following is a simple example problem, with the program for its solution
29// by CVODE. The problem is the semi-discrete form of the advection-diffusion
30// equation in 1-D:
31//
32// du/dt = p1 * d^2u / dx^2 + p2 * du / dx
33//
34// on the interval 0 <= x <= 2, and the time interval 0 <= t <= 5. Homogeneous
35// Dirichlet boundary conditions are posed, and the initial condition is:
36//
37// u(x,t=0) = x(2-x)exp(2x).
38//
39// The nominal values of the two parameters are: p1=1.0, p2=0.5. The PDE is
40// discretized on a uniform grid of size MX+2 with central differencing, and
41// with boundary values eliminated, leaving an ODE system of size NEQ = MX.
42//
43// The program solves the problem with the option for nonstiff systems: ADAMS
44// method and functional iteration. It uses scalar relative and absolute
45// tolerances. In addition to the solution, sensitivities with respect to p1 and
46// p2 as well as with respect to initial conditions are computed for the
47// quantity:
48//
49// g(t, u, p) = int_x u(x,t) at t = 5
50//
51// These sensitivities are obtained by solving the adjoint system:
52//
53// dv/dt = -p1 * d^2 v / dx^2 + p2 * dv / dx
54//
55// with homogeneous Dirichlet boundary conditions and the final condition:
56//
57// v(x,t=5) = 1.0
58//
59// Then, v(x, t=0) represents the sensitivity of g(5) with respect to u(x, t=0)
60// and the gradient of g(5) with respect to p1, p2 is
61//
62// (dg/dp)^T = [ int_t int_x (v * d^2u / dx^2) dx dt ]
63// [ int_t int_x (v * du / dx) dx dt ]
64//
65// This version uses MPI for user routines.
66// Execute with number of processors = N, with 1 <= N <= MX.
67
68#include "mfem.hpp"
69#include <fstream>
70#include <iostream>
71#include <algorithm>
72
73#ifndef MFEM_USE_SUNDIALS
74#error This example requires that MFEM is built with MFEM_USE_SUNDIALS=YES
75#endif
76
77using namespace std;
78using namespace mfem;
79
80// Implement the adjoint rate equations in AdvDiffSUNDIALS
81class AdvDiffSUNDIALS : public TimeDependentAdjointOperator
82{
83public:
84 AdvDiffSUNDIALS(int ydot_dim, int ybdot_dim, Vector p,
85 ParFiniteElementSpace *fes, Array<int> & ess_tdof) :
86 TimeDependentAdjointOperator(ydot_dim, ybdot_dim),
87 p_(p),
88 ess_tdof_list(ess_tdof),
89 pfes(fes),
90 Mf(NULL),
91 M_solver(fes->GetComm())
92 {
93 int skip_zeros = 0;
94
95 cout << "Essential tdofs: " << endl;
96 ess_tdof_list.Print();
97
98 m = new ParBilinearForm(pfes);
100 m->Assemble(skip_zeros);
101 m->Finalize(skip_zeros);
102
103 // Define coefficients
104 mp0 = new ConstantCoefficient(-p_[0]);
105 p0 = new ConstantCoefficient(p_[0]);
106 Vector p2vec(fes->GetParMesh()->SpaceDimension());
107 p2vec = p_[1];
108 p2 = new VectorConstantCoefficient(p2vec);
109
110 k = new ParBilinearForm(pfes);
113 k->Assemble(skip_zeros);
114 k->Finalize(skip_zeros);
115
116 k1 = new ParBilinearForm(pfes);
119 k1->Assemble(skip_zeros);
120 k1->Finalize(skip_zeros);
121
122 M = m->ParallelAssemble();
123 HypreParMatrix *temp = M->EliminateRowsCols(ess_tdof_list);
124 delete temp;
125
126 K = k->ParallelAssemble();
127 temp = K->EliminateRowsCols(ess_tdof_list);
128 delete temp;
129
130 K_adj = k1->ParallelAssemble();
131 temp = K_adj->EliminateRowsCols(ess_tdof_list);
132 delete temp;
133
135 M_solver.SetPreconditioner(M_prec);
136 M_solver.SetOperator(*M);
137
138 M_solver.SetRelTol(1e-14);
139 M_solver.SetAbsTol(0.0);
140 M_solver.SetMaxIter(1000);
141 M_solver.SetPrintLevel(0);
142
143 }
144
145 virtual void Mult(const Vector &x, Vector &y) const;
146
147 virtual void AdjointRateMult(const Vector &y, Vector &yB,
148 Vector &yBdot) const;
149
150 virtual int SUNImplicitSetup(const Vector &y,
151 const Vector &fy, int jok, int *jcur,
152 real_t gamma);
153
154 virtual int SUNImplicitSolve(const Vector &b, Vector &x, real_t tol);
155
156 virtual void QuadratureSensitivityMult(const Vector &y, const Vector &yB,
157 Vector &qbdot) const;
158
159 ~AdvDiffSUNDIALS()
160 {
161 delete m;
162 delete k;
163 delete k1;
164 delete M;
165 delete K;
166 delete K_adj;
167 delete Mf;
168 delete p0;
169 delete mp0;
170 delete p2;
171 }
172
173protected:
174 Vector p_;
175 Array<int> ess_tdof_list;
177
178 // Internal matrices
181 ParBilinearForm *k1;
182
185 HypreParMatrix *K_adj;
186
187 HypreParMatrix *Mf;
188
189 CGSolver M_solver;
190 HypreSmoother M_prec;
194};
195
196// Initial conditions for the problem
198{
199 return x[0]*(2. - x[0])*exp(2.*x[0]);
200}
201
202int main(int argc, char *argv[])
203{
204 // Initialize MPI and HYPRE.
205 Mpi::Init(argc, argv);
206 int myid = Mpi::WorldRank();
207 Hypre::Init();
208
209 // Parse command-line options.
210 int ser_ref_levels = 0;
211 int par_ref_levels = 0;
212 real_t t_final = 2.5;
213 real_t dt = 0.01;
214 int mx = 20;
215 bool step_mode = true;
216
217 int precision = 8;
218 cout.precision(precision);
219
220 OptionsParser args(argc, argv);
221
222 args.AddOption(&mx, "-m", "--mx", "The number of mesh elements in the x-dir");
223 args.AddOption(&ser_ref_levels, "-r", "--refine",
224 "Number of times to refine the mesh uniformly.");
225 args.AddOption(&step_mode, "-a", "--adams", "-no-a","--no-adams",
226 "A switch to toggle between CV_ADAMS, and CV_BDF stepping modes");
227 args.AddOption(&t_final, "-tf", "--t-final",
228 "Final time; start time is 0.");
229 args.AddOption(&dt, "-dt", "--time-step",
230 "Time step.");
231
232 args.Parse();
233 if (!args.Good())
234 {
235 if (myid == 0)
236 {
237 args.PrintUsage(cout);
238 }
239 return 1;
240 }
241 if (myid == 0)
242 {
243 args.PrintOptions(cout);
244 }
245
246 // Create a small 1D mesh with a length of 2. This mesh corresponds with the
247 // cvsAdvDiff_ASA_p_non_p example.
248 Mesh mesh = Mesh::MakeCartesian1D(mx+1, 2.);
249
250 // Refine the mesh to increase the resolution. In this example we do
251 // 'ref_levels' of uniform refinement, where 'ref_levels' is a
252 // command-line parameter. If the mesh is of NURBS type, we convert it to
253 // a (piecewise-polynomial) high-order mesh.
254 for (int lev = 0; lev < ser_ref_levels; lev++)
255 {
256 mesh.UniformRefinement();
257 }
258
259 ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, mesh);
260 mesh.Clear();
261 for (int lev = 0; lev < par_ref_levels; lev++)
262 {
263 pmesh->UniformRefinement();
264 }
265
266 // Finite Element Spaces
267 H1_FECollection fec(1, pmesh->SpaceDimension());
268 ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
269
270 HYPRE_BigInt global_vSize = fes->GlobalTrueVSize();
271 if (myid == 0)
272 {
273 cout << "Number of unknowns: " << global_vSize << endl;
274 }
275
276 // Set up material properties, and primal and adjoint variables
277 // p are the fixed material properties
278 Vector p(2);
279 p[0] = 1.0;
280 p[1] = 0.5;
281
282 // U is the size of the solution/primal vector
283 // Set U with the initial conditions
284 ParGridFunction u(fes);
286 u.ProjectCoefficient(u0);
287
288 cout << "Init u: " << endl;
289 u.Print();
290
291 // TimeDependentOperators need to be TrueDOF Size
292 HypreParVector *U = u.GetTrueDofs();
293
294 // Get boundary conditions
295 Array<int> ess_tdof_list;
296 Array<int> essential_attr(pmesh->bdr_attributes.Size());
297 essential_attr[0] = 1;
298 essential_attr[1] = 1;
299 fes->GetEssentialTrueDofs(essential_attr, ess_tdof_list);
300
301 // Setup the TimeDependentAdjointOperator and the CVODESSolver
302 AdvDiffSUNDIALS adv(U->Size(), U->Size(), p, fes, ess_tdof_list);
303
304 // Set the initial time to the TimeDependentAdjointOperator
305 real_t t = 0.0;
306 adv.SetTime(t);
307
308 // Create the CVODES solver corresponding to the selected step method
309 CVODESSolver *cvodes = new CVODESSolver(fes->GetComm(),
310 step_mode ? CV_ADAMS : CV_BDF);
311 cvodes->Init(adv);
312 cvodes->UseSundialsLinearSolver();
313 cvodes->SetMaxNSteps(5000);
314
315 // Relative and absolute tolerances for CVODES
316 real_t reltol = 1e-8, abstol = 1e-6;
317 cvodes->SetSStolerances(reltol, abstol);
318
319 // Initialize adjoint problem settings
320 int checkpoint_steps = 50; // steps between checkpoints
321 cvodes->InitAdjointSolve(checkpoint_steps, CV_HERMITE);
322
323 // Perform time-integration for the problem (looping over the time
324 // iterations, ti, with a time-step dt).
325 bool done = false;
326 for (int ti = 0; !done; )
327 {
328 real_t dt_real = max(dt, t_final - t);
329 cvodes->Step(*U, t, dt_real);
330 ti++;
331
332 done = (t >= t_final - 1e-8*dt);
333
334 if (done && myid == 0)
335 {
336 cvodes->PrintInfo();
337 }
338 }
339
340 u = *U;
341 if (myid == 0)
342 {
343 cout << "Final Solution: " << t << endl;
344 }
345
346 cout << "u (" << myid << "):" << endl;
347 u.Print();
348 cout << flush;
349 MPI_Barrier(MPI_COMM_WORLD);
350
351 // Calculate the quadrature int_x u dx at t = 5
352 // Since it's only a spatial quadrature we evaluate it at t=5
353 ParLinearForm obj(fes);
354 ConstantCoefficient one(1.0);
356 obj.Assemble();
357
358 real_t g = obj(u);
359 if (myid == 0)
360 {
361 cout << "g: " << g << endl;
362 }
363
364 // Solve the adjoint problem. v is the adjoint solution
365 ParGridFunction v(fes);
366 v = 1.;
367 v.SetSubVector(ess_tdof_list, 0.0);
369
370 // Initialize quadrature sensitivity values to zero
371 Vector qBdot(p.Size());
372 qBdot = 0.;
373
374 t = t_final;
375 cvodes->InitB(adv);
376 cvodes->InitQuadIntegrationB(qBdot, 1.e-6, 1.e-6);
377 cvodes->UseSundialsLinearSolverB();
378 cvodes->SetSStolerancesB(reltol, abstol);
379
380 // Results at time TBout1
381 real_t dt_real = max(dt, t);
382 cvodes->StepB(*V, t, dt_real);
383 if (myid == 0)
384 {
385 cout << "t: " << t << endl;
386 }
387
388 cout << "v (" << myid << "):" << endl;
389 V->Vector::Print();
390 cout << flush;
391 MPI_Barrier(MPI_COMM_WORLD);
392
393 // Evaluate the sensitivity
394 cvodes->EvalQuadIntegrationB(t, qBdot);
395
396 MPI_Barrier(MPI_COMM_WORLD);
397 if (myid == 0)
398 {
399 cout << "sensitivity:" << endl;
400 qBdot.Print();
401 }
402
403 // Free the used memory.
404 delete fes;
405 delete pmesh;
406 delete U;
407 delete V;
408 delete cvodes;
409
410 return 0;
411}
412
413// AdvDiff rate equation
414void AdvDiffSUNDIALS::Mult(const Vector &x, Vector &y) const
415{
416 Vector z(x.Size());
417 Vector x1(x);
418
419 // Set boundary conditions to zero
420 x1.SetSubVector(ess_tdof_list, 0.0);
421
422 K->Mult(x1, z);
423
424 y = 0.;
425 M_solver.Mult(z, y);
426}
427
428// AdvDiff Rate equation setup
429int AdvDiffSUNDIALS::SUNImplicitSetup(const Vector &y,
430 const Vector &fy,
431 int jok, int *jcur, real_t gamma)
432{
433 // Mf = M(I - gamma J) = M - gamma * M * J
434 // J = df/dy => K
435 *jcur = 1; // We've updated the jacobian
436
437 delete Mf;
438 Mf = Add(1., *M, -gamma, *K);
439 HypreParMatrix *temp = Mf->EliminateRowsCols(ess_tdof_list);
440 delete temp;
441 return 0;
442}
443
444// AdvDiff Rate equation solve
445int AdvDiffSUNDIALS::SUNImplicitSolve(const Vector &b, Vector &x, real_t tol)
446{
447 Vector z(b.Size());
448 M->Mult(b,z);
449
450 CGSolver solver(pfes->GetComm());
451 HypreSmoother prec;
453 solver.SetPreconditioner(prec);
454 solver.SetOperator(*Mf);
455 solver.SetRelTol(1E-14);
456 solver.SetMaxIter(1000);
457
458 solver.Mult(z, x);
459
460 return (0);
461}
462
463// AdvDiff adjoint rate equation
464void AdvDiffSUNDIALS::AdjointRateMult(const Vector &y, Vector & yB,
465 Vector &yBdot) const
466{
467 Vector z(yB.Size());
468
469 // Set boundary conditions to zero
470 yB.SetSubVector(ess_tdof_list, 0.0);
471 K_adj->Mult(yB, z);
472 M_solver.Mult(z, yBdot);
473}
474
475// AdvDiff quadrature sensitivity rate equation
476void AdvDiffSUNDIALS::QuadratureSensitivityMult(const Vector &y,
477 const Vector &yB,
478 Vector &qBdot) const
479{
480 // Now we have both the adjoint, yB, and y, at the same point in time
481 // We calculate
482 /*
483 * to u(x, t=0) and the gradient of g(5) with respect to p1, p2 is
484 * (dg/dp)^T = [ int_t int_x (v * d^2u / dx^2) dx dt ]
485 * [ int_t int_x (v * du / dx) dx dt ]
486 */
487
488 ParBilinearForm dp1(pfes);
489 ConstantCoefficient mone(-1.);
490 dp1.AddDomainIntegrator(new DiffusionIntegrator(mone));
491 dp1.Assemble();
492 dp1.Finalize();
493
494 HypreParMatrix * dP1 = dp1.ParallelAssemble();
495 HypreParMatrix *temp = dP1->EliminateRowsCols(ess_tdof_list);
496 delete temp;
497
498 Vector b1(y.Size());
499 dP1->Mult(y, b1);
500 delete dP1;
501
502 ParBilinearForm dp2(pfes);
503 Vector p2vec(pfes->GetParMesh()->SpaceDimension()); p2vec = 1.;
504 VectorConstantCoefficient dp2_coef(p2vec);
505 dp2.AddDomainIntegrator(new ConvectionIntegrator(dp2_coef));
506 dp2.Assemble();
507 dp2.Finalize();
508
509 HypreParMatrix * dP2 = dp2.ParallelAssemble();
510 temp = dP2->EliminateRowsCols(ess_tdof_list);
511 delete temp;
512
513 Vector b2(y.Size());
514 dP2->Mult(y, b2);
515 delete dP2;
516
517 real_t dp1_result = InnerProduct(pfes->GetComm(), yB, b1);
518 real_t dp2_result = InnerProduct(pfes->GetComm(), yB, b2);
519
520 qBdot[0] = -dp1_result;
521 qBdot[1] = -dp2_result;
522}
real_t u_init(const Vector &x)
int Size() const
Return the logical size of the array.
Definition array.hpp:144
void Print(std::ostream &out=mfem::out, int width=4) const
Prints array to stream with width elements per row.
Definition array.cpp:23
virtual void Finalize(int skip_zeros=1)
Finalizes the matrix initialization if the AssemblyLevel is AssemblyLevel::LEGACY....
void AddDomainIntegrator(BilinearFormIntegrator *bfi)
Adds new Domain Integrator. Assumes ownership of bfi.
Conjugate gradient method.
Definition solvers.hpp:513
virtual void SetOperator(const Operator &op)
Also calls SetOperator for the preconditioner if there is one.
Definition solvers.hpp:526
virtual void Mult(const Vector &b, Vector &x) const
Iterative solution of the linear system using the Conjugate Gradient method.
Definition solvers.cpp:718
void EvalQuadIntegrationB(double t, Vector &dG_dp)
Evaluate Quadrature solution.
void InitB(TimeDependentAdjointOperator &f_)
Initialize the adjoint problem.
virtual void Step(Vector &x, double &t, double &dt)
void UseSundialsLinearSolverB()
Use built in SUNDIALS Newton solver.
void Init(TimeDependentAdjointOperator &f_)
virtual void StepB(Vector &w, double &t, double &dt)
Solve one adjoint time step.
void InitQuadIntegrationB(mfem::Vector &qB0, double reltolQB=1e-3, double abstolQB=1e-8)
Initialize Quadrature Integration (Adjoint)
void InitAdjointSolve(int steps, int interpolation)
Initialize Adjoint.
void SetSStolerancesB(double reltol, double abstol)
Tolerance specification functions for the adjoint problem.
void SetSStolerances(double reltol, double abstol)
Set the scalar relative and scalar absolute tolerances.
Definition sundials.cpp:875
void SetMaxNSteps(int steps)
Set the maximum number of time steps.
Definition sundials.cpp:899
void PrintInfo() const
Print various CVODE statistics.
Definition sundials.cpp:919
void UseSundialsLinearSolver()
Attach SUNDIALS GMRES linear solver to CVODE.
Definition sundials.cpp:855
A coefficient that is constant across space and time.
Class for domain integration .
Definition lininteg.hpp:109
A general function coefficient.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
Wrapper for hypre's ParCSR matrix class.
Definition hypre.hpp:388
void EliminateRowsCols(const Array< int > &rows_cols, const HypreParVector &X, HypreParVector &B)
Definition hypre.cpp:2356
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, real_t alpha=1.0, real_t beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition hypre.cpp:1815
Wrapper for hypre's parallel vector class.
Definition hypre.hpp:206
void Print(const char *fname) const
Prints the locally owned rows in parallel.
Definition hypre.cpp:413
Parallel smoothers in hypre.
Definition hypre.hpp:1020
void SetType(HypreSmoother::Type type, int relax_times=1)
Set the relaxation type and number of sweeps.
Definition hypre.cpp:3550
static void Init()
Initialize hypre by calling HYPRE_Init() and set default options. After calling Hypre::Init(),...
Definition hypre.hpp:74
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPreconditioner(Solver &pr)
This should be called before SetOperator.
Definition solvers.cpp:173
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
void AddDomainIntegrator(LinearFormIntegrator *lfi)
Adds new Domain Integrator. Assumes ownership of lfi.
Mesh data type.
Definition mesh.hpp:56
Array< int > bdr_attributes
A list of all unique boundary attributes used by the Mesh.
Definition mesh.hpp:282
static Mesh MakeCartesian1D(int n, real_t sx=1.0)
Creates 1D mesh, divided into n equal intervals.
Definition mesh.cpp:4235
void Clear()
Clear the contents of the Mesh.
Definition mesh.hpp:730
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
void UniformRefinement(int i, const DSTable &, int *, int *, int *)
Definition mesh.cpp:10970
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
static void Init(int &argc, char **&argv, int required=default_thread_required, int *provided=nullptr)
Singleton creation with Mpi::Init(argc, argv).
void Parse()
Parse the command-line options. Note that this function expects all the options provided through the ...
void PrintUsage(std::ostream &out) const
Print the usage message.
void PrintOptions(std::ostream &out) const
Print the options.
void AddOption(bool *var, const char *enable_short_name, const char *enable_long_name, const char *disable_short_name, const char *disable_long_name, const char *description, bool required=false)
Add a boolean option and set 'var' to receive the value. Enable/disable tags are used to set the bool...
Definition optparser.hpp:82
bool Good() const
Return true if the command line options were parsed successfully.
Class for parallel bilinear form.
HypreParMatrix * ParallelAssemble()
Returns the matrix assembled on the true dofs, i.e. P^t A P.
void Assemble(int skip_zeros=1)
Assemble the local matrix.
Abstract parallel finite element space.
Definition pfespace.hpp:29
MPI_Comm GetComm() const
Definition pfespace.hpp:273
void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1) const override
HYPRE_BigInt GlobalTrueVSize() const
Definition pfespace.hpp:285
ParMesh * GetParMesh() const
Definition pfespace.hpp:277
Class for parallel grid function.
Definition pgridfunc.hpp:33
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Class for parallel linear form.
void Assemble()
Assembles the ParLinearForm i.e. sums over all domain/bdr integrators.
Class for parallel meshes.
Definition pmesh.hpp:34
virtual void SetTime(const real_t t_)
Set the current time.
Definition operator.hpp:360
Vector coefficient that is constant in space and time.
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
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:604
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
int main()
HYPRE_Int HYPRE_BigInt
real_t b
Definition lissajous.cpp:42
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
real_t InnerProduct(HypreParVector *x, HypreParVector *y)
Definition hypre.cpp:439
float real_t
Definition config.hpp:43
void Add(const DenseMatrix &A, const DenseMatrix &B, real_t alpha, DenseMatrix &C)
C = A + alpha*B.
real_t p(const Vector &x, real_t t)
RefCoord t[3]