L2 projection after non-uniform refinement

[jchan192]

Hi all, The question here is similar to #1072 . Given a Gridfunction u_gf , I want to project u_gf to new finite element space, where the mesh is refined. Unlike #1072, the refinement is not performed to the whole mesh using UniformRefinemen, but only at certain regions using the Thresholdrefiner.

I wonder if L2ProjectionGridTransfer can still project u_gf to space that is not refined uniformly.

Thank you

[psocratis]

Hi @jchan192,

after a mesh refinement (uniform or not) you need to call FiniteElementSpace::Update(). This updates the FiniteElementSpace but also creates internally the transfer operator (the embedding) from the old FiniteElementSpace to the new one. If you then call u_gf.Update(), this transfer operator is used to map the u_gf to the new FiniteElementSpace.

Hope this helps.

Best,

Socratis

[jchan192]

Thanks for the reply! @psocratis Yes. My current code is using u_gf.Update() after refinement. However, the new problem that I am solving now requires using projection operator, L2- projection more specifically.

Please let me extend my question here. What kind of interpolation scheme does u_gf.Update() adopt? Is it nodal interpolation or projection?

[psocratis]

Hi @jchan192,

The Update() uses nodal interpolation. I still fail to see the need for the projection. If you are doing a refinement then the coarse space is a subspace of the fine space, hence the solution already belongs to the fine space. A Projection would give you back the same thing, i.e., the embedding of u_gf in the fine space.

Unless you are doing LOR, where the refined space is of lower order, then I don't see the need for a projection. is that the case?

[tzanio]

Just to add to @psocratis's answer, the general intent of the GridTransfer classes is to be used for transfer between say a high-order space on a high-order mesh and a low-order space on a low-order refined mesh, see https://arxiv.org/abs/2103.05283.

That being said, one can consider using them with AMR, e.g. for derefinement. Do you need to do something like this?

If the answer is no, you are much better off with the built-in AMR interpolation and restriction algorithms.

[jchan192]

Thanks @psocratis @tzanio for the reply. I am not doing a transfer from HO space and mesh to LO space and refined mesh. So looks like I actually don't need GridTransfer. Thanks for clarifying.

I am implementing a modified Crank-Nicolson solver for Schorindger Eq #2811. In pointwise form, it looks like this where U^n is the approximated solution at timestep n and π is the projection/ interpolation from space V^{n-1} to V^{n}. Literature applied L2 projection for π, which is why I brought up the topic here.

So looks like Update() is all I need, since Im not dealing with LOR. The problem is the RHS term with π * Laplacian U^{n-1}. In finite element, does it mean I will need to perform an Update() after applying the DiffusionIntergrator?

This is a short demonstration of my plan. Do you think the implementation look correct?

  // U is actually complex valued, but lets assume real here for simplicity
   K = new ParBilinearForm(&fespace);                                                                                                                                                                              
   K->AddDomainIntegrator(new DiffusionIntegrator(KCst));                                                                                                                                                          
   K->Assemble();                                                                                                                                                                                                  
   K->Finalize();                                                                                                                                                                                                  
   Kmat = K->ParallelAssemble(); 

   Kmat->Mult(U,Utmp); // Laplacian U^{n-1}
   RHS_gf.SetFromTrueDofs(Utmp); // vector cant be interpolated, but gridfunction can

   // After refinement....

   fespace.Update();
   RHS_gf.Update(); // interpolated Laplacian U^{n-1}?
   fespace.UpdatesFinished();
   Utmp.Update(fespace.GetTrueVSize());
   RHSU_gf.GetTrueDofs(Utmp); // Utmp now is π * Laplacian U^{n-1} ?

Thanks all for your time. Since I have already implemented the standard Crank-Nicolson code, I believe that the projection/interpolation on the Laplacian is the main problem here.

[jchan192]

I will close the issue here, because the question is answered. I will open a new issue regarding the question above. Thank you all

[jchan192]

@psocratis sorry for coming back. I have consulted with an FE analysis expert regarding my problem. I now have a better understanding but still couldn't implement it correctly in MFEM....

We suppose to compute the following L2 projection of ΔU from a mesh to a refined/derefined mesh (same order of FE space). You suggest that interpolation from Update will give the same effect, because the refinement is performed hierarchically, I think. But How should I interpolate ΔU correctly? I have an example above, where I get the interpolation of ΔU from 1) multiplying U with the stiffness matrix 2) then put it into the gridfunction 3) do an update.

Does this approach sound good?

One more last question. I am actaully also derefining. Do we expect Update in the case to give same result as L2 projection?

Thank you for your time. I really appreciate it if there is any comment.

[psocratis]

Hi @jchan192,

I think what you are doing would not work. Multiplying U with stiffness gives you a dual vector on the coarse space. The update method is intended to work for a primal vector, i.e., a GridFunction.

So I see the following choices.

1. Perform the Updateon U^{n-1} and then compute its Laplacian on the fine space.
2. Compute the Gradient of U^{n-1} on the coarse space using the GradientInterpolator. This will give you a GridFunction in H(curl) (ND) (or in Vector L2, your choice). You then perform the Update on the gradient and you get an ND GridFunction on the fine space. With this you can construct a VectorGridFunctionCoeffiecient and use it with DomainLFGradIntegrator to get your RHS.
3. Lastly you can produce a primal vector out of your dual vector in the coarse space by solving a projection problem. i.e., for L2-projection you multiply the dual vector with the inverse of the mass matrix. This would give you an H1 representation of the Laplacian and you can then use the Update to transfer to the fine space, and finally you multiply with the stiffness matrix on the fine space to get your RHS.

In case of a derefinement the Update will perform a projection.

Hope this helps.

Best,

Socratis

[jchan192]

@psocratis Thanks for the reply!

My problem requires specificlly applying the laplacian on the coarse space first then interpolate on fine, instead of the reverse. So 1) may not be what I want.

I think the other two solutions seem to be what im looking for!

2) I am computing a laplacian, but the GradientInterpolatoronly computes the 1st derivative. For this method, I don't understnad which step acutally computes the laplacian.

3) I think I have tried this method, partially. Let me clarify that U is a primal vector. you suggested me to 3.1- first solve the projection problem based on the dual vector, which is the stiffness matrix^{n-1} U^{n-1}. 3.2 - Update the projected primal vector to fine space 3.3 - multiply stifness matrix on fine space, to get RHS My quesiton is why would 3.3 be neccessary , when we already did stiffness matrix U at step 3.1?

Thanks again for the help.

[psocratis]

Hi @jchan192,

For (2), you compute the Gradient but then when you transfer to the fine space you compute the Laplacian in the weak sense, i.e., (\nabla U, \nabla \v) (\nabla U is what you transferred). If you really need the Laplacian to be computed on the coarse grid, then this would not do it.

For (3), yes you are right, I meant multiply by the mass matrix (not the stiffness). Alternatively, you can use it in a GridFunctionCoefficient with a DomainLFintegrator.

Best, Socratis

[jchan192]

Hi @psocratis Based on your suggestion, I implemented the following code. I didn't show the full code, because its too complicated, but do you think below is what you suggested?

I still don't understand why we will need to do step 3.3, multiplying the updated vector by mass matrix.

   // prepare Mass matrix for projection     
    Mproj = new ParBilinearForm(&fespace);                                                                                                                                                                          
   Mproj->AddDomainIntegrator(new MassIntegrator(one));                                                                                                                                                            
   Mproj->Assemble();                                                                                                                                                                                              
   Mproj->Finalize();                                                                                                                                                                                              
   Mpmat = Mproj->ParallelAssemble();                                                                                                                                                     

   // prepare stiffness matrix                                                                                                                                                  
   K = new ParBilinearForm(&fespace);                                                                                                                                                                              
   K->AddDomainIntegrator(new DiffusionIntegrator(KCst));                                                                                                                                                          
   K->Assemble();                                                                                                                                                                                                  
   K->Finalize();                                                                                                                                                                                                  
   Kmat = K->ParallelAssemble();   

   //  Step 3.1: To do on coarse grid    //
   Utmp.Update(block_toffsets);                                                                                                                                                                                                                                                 
   Kmat->Mult(U,Utmp); // laplacian U^{n-1} but dual vector                                                                                                                                                                                                      

   // l2 projection to obtain primal vector
   CGSolver *l2projection = new CGSolver(MPI_COMM_WORLD);                                                                                                                                                                                                                       
   Vector Uprim(Utmp); 
   l2projection->SetRelTol(1e-12);                                                                                                                                                                                                                                              
   l2projection->SetAbsTol(1e-12);                                                                                                                                                                                                                                              
   l2projection->SetMaxIter(150);                                                                                                                                                                                                                                               
   l2projection->SetPrintLevel(-1);                                                                                                                                                                                                                                             
   l2projection->SetOperator(*Mpmat);                                                                                                                                                                                                                                           
   l2projection->iterative_mode = false;                                                                                                                                                                                                                                        
   l2projection->Mult(Utmp,Uprim); // get primal vector Utmp2 

   RHS_gf.SetFromTrueDofs(Uprim);   // vector cant be updated, so put it into Gridfunction RHS_gf

   ////////////////   
   //    Step 3.2: After refining and derefining mesh, Update RHS.........  //
   ///////////////

   //     Step 3.3: To do on fine grid    //
   Vector Uproj(Uprim);
   RHS_gf.GetTrueDofs(Uprim); //  Get Projected laplacian U^{n-1}. Done?
   Mpmat->Mult(Uprim,Uproj); // Why do we need this step?

Thank you for your time.

[psocratis]

@jchan192,

after the refinement/derefinement you transfer the RHS_gf which is now a GridFunction using RHS_gf.Update(). Now you have a GridFuction (or a primal vector) in the refined space. If you want to produce the RHS i.e. (RHS_gf, v) you can either multiply it by the mass matrix on the fine space or you can use a LinearForm along with a DomainLFIntegrator and a GridFunctionCoefficient produced by RHS_gf.

Best,

Socratis