Hybrid velocity boundary conditions

ASPECT
1

Dear all,

I’m developing a data assimilation method with hybrid velocity boundary conditions based on ASPECT. To study the internal stress in the continental plates, I’m trying to fix the boundary velocities on the top surface above oceanic plates and apply a free slip boundary condition for the top surface above continental lithosphere. To realize the purpose, I add such a function in Simulator class:


  template <int dim>
  void
  Simulator <dim>::update_hybrid_vbc_boundary_ids()
  {
    const std::vector<Point<dim>> points = finite_element.get_unit_support_points();
    const Quadrature<dim> quadrature(points);
    FEValues<dim> fe_values (finite_element, quadrature, update_quadrature_points);
    typename DoFHandler<dim>::active_cell_iterator cell;

    for (cell = dof_handler.begin_active();
         cell != dof_handler.end(); ++cell)
    {
      if (!cell->is_artificial())
      {
        for (const unsigned int face_number : cell->face_indices())
        {
          const typename DoFHandler<dim>::face_iterator face = cell->face(face_number);
          if (face->at_boundary())
          {
            if (face->boundary_id() != 0)
            {
              face->set_boundary_id(1);
              fe_values.reinit(cell);
              for (unsigned int q = 0; q < quadrature.size(); ++q)
              {
                const Point<dim> position = fe_values.quadrature_point(q);
                const double lithosphere_age = get_lithosphere_age (position);
                const double distance_from_trench = get_distance_from_trench (position);
                if (!(lithosphere_age < 400 && std::abs(distance_from_trench) > 0.4))
                {
                  face->set_boundary_id(2);
                  break;
                }
              }
            }
          }
        }
      }
    }
  }

The function get_lithosphere_age is to determine oceanic or continental plates and the function get_distance_from_trench is to determine the distance of a point from the trench. The top surface is divided into two boundaries, 1 and 2. One for the gplates boundary and one for the free slip boundary.

However, the results show that the velocities have been applied into the results, but not constrained in the linear systems (see the Figure).

I’d like to know how to apply the hybrid velocity boundary conditions correctly.

Thanks a lot!

figure1
figure1822×456 24.7 KB

figure2
figure2822×462 33.2 KB

figure3
figure3990×556 38.8 KB

2

Hi @lixinyv,

Thank you for posting the question and associated code to the forum.

A quick response:

  1. There is a cookbook (prescribed velocity) that shows how to constrain velocity (or other) degrees of freedom within the model solution using a plugin. Is your approach based on this plugin? If not, I suggest looking through the prescribed velocity and seeing if that can be modified to accomplish your intended goals.
  2. Can you expand on how the figures show the applied velocity is not being constrained in the linear system?

Cheers,

John

3

Hi John,

Thanks for your reply!

I’m trying using the cookbook called prescribed velocity. However, I met some difficulties when compiling the plugin.

/fs2/home/ljliu01/plugin/prescribed_velocity.cc:1:10: fatal error: deal.II/base/parameter_handler.h: No such file or directory
    1 | #include <deal.II/base/parameter_handler.h>
      |          ^~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
compilation terminated.
make[2]: *** [CMakeFiles/prescribed_velocity.dir/build.make:76: CMakeFiles/prescribed_velocity.dir/prescribed_velocity.cc.o] Error 1
make[1]: *** [CMakeFiles/Makefile2:111: CMakeFiles/prescribed_velocity.dir/all] Error 2

Also I add this plugin into the Simulator class in core.cc, but it doesn’t work, either.

For the second questions, you can see the first figure, velocity magnitude. You can find that the surface velocity on the top of oceanic plate (left part) is the same with the real velocity. However, in the internal part of oceanic plates, there is a large discrepancy with the surface velocity.

Thanks a lot!
Xinyu

4

@lixinyv How are you compiling the plugin? Are you just calling the compiler yourself, or are you using the CMake files provided with ASPECT to compile plugins?

Best

W.