Closed tu1620808124 closed 1 year ago
jwpeterson
commented
1 year ago I'm not exactly sure what you meant by this:
substitute the coordinates into the exact solution， I get the interpolate values
I think you need to use the function:
System::point_value(unsigned int var, const Point & p);to evaluate the L2_LAGRANGE solution at the desired points. Can you try that and let us know if it works?
roystgnr
commented
1 year ago if my understanding of L2_LAGRANGE is wrong
One misunderstanding, though I'm not sure if it's you misunderstanding L2_LAGRANGE or if I've misunderstood you:
There are actually 9 different coordinates on the grid
Though there are 9 different coordinates on the grid, there are 16 independent values of an L2_LAGRANGE solution at those coordinates. It's a discontinuous FE type, so at nodes shared by multiple elements it doesn't take a unique value; instead it has a unique directional limit from each of the elements sharing that node.
And one misunderstanding, of project_solution, which probably explains your results:
project_solution is a complicated projection of the solution in general, not an interpolation of it. For elements with vertex degrees of freedom, we do interpolate at those vertices - but discontinuous elements don't have vertex degrees of freedom, so we don't do that there. For LAGRANGE elements, we also interpolate at other nodes - but L2_LAGRANGE != LAGRANGE, so we don't do that there. For any other elements, we do local Hilbert projections (H1 for C1 elements, L2 for others) to evaluate non-vertex DoFs: first along edges, holding vertex values constant, then along faces, holding edge+vertex values constant, then on element interiors, holding edge+vertex+face values constant. A discontinuous FE type like L2_LAGRANGE only has element interior DoFs, so the value on any element ought to be the FE value that minimizes L2 error vs the input (based on the default quadrature rule, not exact integration). So we do expect 16 distinct values for those 16 DoFs, at least in the case where the input isn't already in your FE space. The vertex values won't be an exact match at vertices, because the projection can add a little error there to obtain a greater reduction in error over the rest of the element.
roystgnr
commented
1 year ago Even if this is just a misunderstanding or two, I really appreciate the test case to make your issue easy to replicate. And any misunderstanding here is my fault; I really need to get this into the documentation. If you think the paragraph above is readable enough then I'll put something like it into the GenericProjector docs and I'll add a reference to it to project_solution, project_vector, etc.
tu1620808124
commented
1 year ago #include <iostream>
#include <algorithm>
#include <sstream>
#include <math.h>
#include <cstring>
// Various include files needed for the mesh & functionality.
#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/elem.h"
#include "libmesh/exodusII_io.h"
#include "libmesh/equation_systems.h"
#include "libmesh/quadrature_gauss.h"
#include "libmesh/getpot.h"
#include "libmesh/parameters.h"
#include "libmesh/numeric_vector.h"
#include "libmesh/dof_map.h"
// The nonlinear solver and system we will be using
#include "libmesh/nonlinear_solver.h"
#include "libmesh/nonlinear_implicit_system.h"
#include "libmesh/transient_system.h"
#include "libmesh/enum_solver_package.h"
// Necessary for programmatically setting petsc options
#ifdef LIBMESH_HAVE_PETSC
#include <petsc.h>
#include "libmesh/petsc_macro.h"
#include "libmesh/petsc_matrix.h"
#endif
#include "libmesh/getpot.h"
// error_estimate
#include "libmesh/exact_solution.h"
#include "libmesh/exact_error_estimator.h"
#include "libmesh/error_vector.h"
using namespace libMesh;
Real exact_solution(const Real x,
const Real y,
const Real t)
{
return 2*sin(x)+sin(y)+t;;
}
Number init_value(const Point &p, const Parameters ¶meters, const std::string &, const std::string &var_name)
{
Real time = parameters.get<Real>("time");
return exact_solution(p(0), p(1), time);
}
void init_cd(EquationSystems &es, const std::string &system_name)
{
TransientNonlinearImplicitSystem &system = es.get_system<TransientNonlinearImplicitSystem>("project");
system.project_solution(init_value, NULL, es.parameters);
}
int main(int argc, char **argv)
{
LibMeshInit init(argc, argv);
libmesh_example_requires(libMesh::default_solver_package() == PETSC_SOLVERS ||
libMesh::default_solver_package() == TRILINOS_SOLVERS,
"--enable-petsc or --enable-trilinos");
// define the mesh
Mesh mesh(init.comm());
const unsigned int nx = 2;
MeshTools::Generation::build_cube(mesh, nx, nx, 0, 0., 1., 0.,
1., 0., 0., QUAD4);
mesh.all_first_order();
mesh.print_info();
unsigned int order = 1;
EquationSystems equation_systems(mesh);
TransientNonlinearImplicitSystem &system =
equation_systems.add_system<TransientNonlinearImplicitSystem>("project");
system.time = 0.;
equation_systems.parameters.set<Real>("time") = system.time;
system.add_variable("u", static_cast<Order>(order), L2_LAGRANGE);
const unsigned int u_var = system.variable_number("u");
system.attach_init_function(init_cd);
equation_systems.init();
const Point p(0,0,0);
std::cout << system.point_value(u_var, p) << std::endl;
std::cout << exact_solution(0,0,0) << std::endl;
}0.0103302
0For example, when my exact_solution function expression is $f (x, y, t)= 2sin(x)+sin(y)+t$, and everyone knows $f(0,0,0) = 0$, but use the project_solution, the result can be $0.0103302$ it seems it was not the solution in position (0,0), maybe the value of the nearest Gauss–Legendre point.
jwpeterson
commented
1 year ago maybe the value of the nearest Gauss–Legendre point.
Like @roystgnr said, the projection for L2_LAGRANGE is based on minimizing the L2-error between the true solution and the projected solution, so there is nothing that guarantees it will interpolate the true solution value. I would suggest trying your code again while increasing nx and you should see that the result converges to 0.
tu1620808124
commented
1 year ago I seem to understand most of it, When I am doing the convergence test, will this error have a relatively large impact? I looked at the shape function of L2_LAGRANGE and it seems to be the same as LANGRANGE, Can I write an interpolation function myself as an initial value?
roystgnr
commented
1 year ago As an aside: how did I not notice when Github added MathJax support? That's awesome!
Anyway: IIRC this is equivalent to interpolation at Gaussian quadrature points for special cases (e.g. biquadratic onto bilinear), but for something like a trig function onto bilinear I don't know of a simpler interpretation than L2 projection.
Imagine you're in 1D, you're trying to project $f(x)=1-x^2$ to an L2_LAGRANGE space, and you just have one EDGE2 on $[-1,1]$. Obviously the interpolation is $f_h(x)=0$ ... but that's kind of a lousy projection, isn't it? $||f-f_h||_2 = \sqrt{16/15}$? Take $f_h(x)=2/3$, way more than 0.01 off in this case ... and you get $||f-f_h||_2 = \sqrt{8/45}$, much better over most of the element, at the cost of being much worse at the nodes. (But if you want to minimize error at nodes, why not use LAGRANGE? The whole point of the discontinuity is to add degrees of freedom that lets your approximation take on multiple values on inter-element boundaries.)
You can definitely interpolate yourself if you prefer. Loop over elements, and on each element with L2_LAGRANGE the shape function indexed by dof_indices[i] is the Lagrange function equal to $1$ at elem->point(i) and $0$ at the other nodes, so just evaluate at that point and set() your solution vector with the result, then when all done do a close() to sync it up with other processors in parallel and then an update() to get your current_local_solution in sync.
In the region $(0,1)\times (0, 1)$, The grid is $2\times2$, when I use L2_LAGRANGE to attach the init condition, the 16-point value of the project values are
There are actually 9 different coordinates on the grid， substitute the coordinates into the exact solution， I get the interpolate values:
the difference between interpolated and projected values seems too big, For example, the value at (0,0) is obviously different. when I use LAGRANGE the project values are:
It seems the project results on LANGRANGE are right. I don't know if this means the projection is not accurate on L2_LAGRANGE, or if my understanding of L2_LAGRANGE is wrong,