A redundant inverse sign in the comment here?

[Shangxin-Liu]

When I read through this part, the inverse sign in this comment seems redundant because the left side is angular momentum and the right side is moment of inertial times angular velocity (\omega)?

[Shangxin-Liu]

Another thing I'm a little confused. From the line 428 to 430, I see the definition on the I_xy, I_xz and I_yz of the non-diagonal components of the moment of inertia. Should we have a negative sign before I_xy, I_xz and I_yz here because the definition of the moment of inertial tensor is supposed to be
I_xx -I_xy -I_xz
-I_xy I_yy -I_yz
-I_zx -I_zy I_zz ?

I cannot find the negative sign in the later code. The negative sign at the line 460 seems means corresponding "subtracting" at the line 465. Maybe I made a mistake or the moment of inertial has a different definition, but can you let me know :) ? @gassmoeller @tjhei

[gassmoeller]

I am far from an expert in this part of the code, but I think you might be right. If you are right the 2D version of the angular momentum and net rotation removal worked, but the 3D version was wrong and would create a velocity correction that would not cancel the angular momentum (or net rotation)?
Could you test your idea by creating a 'rotation_statistics' postprocessor that computes net rotation and angular momentum? It should be relatively simple by copying the 'velocity statistics' and replace lines 64-67 by something that computes net rotation and angular momentum (as done in nullspace.cc 413 and 422). If you are correct the current geoid.prm test should have a varying angular momentum, and should have a constant angular momentum after fixing the tensor entries. Let me know if you do not have the time to finish this until the weekend, then I will take a stab at it. I would like to check this before the release.
Nice job finding this! And another reason why we need as many tests as possible.

[ian-r-rose]

You are right that the inverse sign in the comment is a typo (the correct treatment is in the manual here).

Also, looks like you are right that the products of inertia have a sign error. I think the reason we didn't notice earlier was that the products of inertia are sooooo much smaller than the diagonal terms (by a factor of maybe 10^-5) that the sign error didn't prevent the algorithm from still mostly working.

[Shangxin-Liu]

I have some other stuff to do so feel free to take a stab at this @gassmoeller . This indeed seems influencing the results pretty little but anyway, let us fix this:)

The reason why I looked at this patch is that I'd like to think of the possibility to add one more removing net rotation option to remove net rotation at each depth separately, instead of overall for the whole domain. Sometimes when the shallower part of the domain (say lithosphere) has a lateral viscosity variation, even if the net rotation is removed overall for the whole domain, the top part with LVV may still have a net rotation over the underlying deeper part. I guess Ian wrote this part to remove net rotation, so any insight/suggestions on how to write one more option to remove net rotation at each depth separately? @ian-r-rose From the line 463 to 465, the subtraction of the angular velocity is applied to the velocity of the whole domain. But any idea on how to apply this kind of subtraction to the velocity at each depth separately by "interpolate_onto_velocity_system"? AMR feature of ASPECT may be a challenge here. Or maybe the best way is to do this by my own after running ASPECT? @ian-r-rose

[gassmoeller]

I will take a look at it tomorrow. Lets test and merge this first (since it is only documentation).

/run-tests

[bangerth]

On 05/02/2018 07:11 PM, Shangxin Liu wrote: The reason why I looked at this patch is that I'd like to think of the possibility to add one more removing net rotation option to remove net rotation at each depth separately, instead overall for the whole domain. Sometimes when the shallower part of the domain (say lithosphere) has a lateral viscosity variation, even if the net rotation is removed overall for the whole domain, the top part with LVV may still have a net rotation over the underlying deeper part. I guess Ian wrote this part to remove net rotation, so any insight/suggestions on how to write one more option to remove net rotation at each depth separately? @ian-r-rose <https://github.com/ian-r-rose> From the line 463 to 465, the subtraction of the angular velocity is applied to the velocity of the whole domain. But any idea on how to apply this kind of subtraction to the velocity at each depth separately by "interpolate_onto_velocity_system"? AMR feature of ASPECT may be a challenge here. Or maybe the best way is to do this by my own after running ASPECT? @ian-r-rose <https://github.com/ian-r-rose>

Shangxin, I think before you think about how to implement this, you'll have to define what exactly you want to do. The net rotation correction adds a velocity field to the current solution that at every point corresponds to a solid body rotation. If you want to do this layer by layer, you will end up with adding a field that is not continuous across layer interfaces. That's presumably not what you want. Can you explain how you want to transition from the added velocity field in one layer to that of the next layer?

[Shangxin-Liu]

Hi Wolfgang, @bangerth
I'm still thinking about this. Directly removing the net rotation depth by depth is perhaps not a clever way. Let me further describe what I want to do in detail. Removing net rotation from the 3D whole domain is indeed physically reasonable and is a commonly used procedure. However, if lateral viscosity variation exists in the top shallower part (say lithosphere, which is commonly used in global model), even if the net rotation of the whole domain has been removed, the lithosphere my still has a net rotation with respect to the underlying mantle. So I'd like to perform a "second" net rotation removal step that further removing the net rotation within the lithosphere sublayer. After this, we can conveniently compare the output surface velocity field from ASEPCT with the plate motion observation in no-net-rotation frame. I'm thinking of adding this second net rotation removal in the shallower sublayer domain (say top 200 km). From the current null_space.cc code, I see that I can compute the rotation vector (\omega) from the lithosphere by defining the integral within the top shallower sublayer. However, I have no idea how to subtract this rotation vector only from the velocity field in the shallower part. From line 463 to 465, the subtraction of the angular velocity is applied to the velocity of the whole domain. But how can I extract the velocity solution field above an input depth to construct rotation like line 463 does, and perform the removal of the net rotation within this sublayer following what line 464-465 does?

[bangerth]

You would define a class like the Rotation class used there that only returns a rotational velocity if $r&gt;R0$ and a zero velocity otherwise.

But you didn't address my question. You are focused on how to do things that you forget what you want to do. If you subtract a rotation only for the lithosphere, then you get an infinite strain (a discontinuous velocity) at the interface between lithosphere and the region below. This is not physical.

[Shangxin-Liu]

Yes, Wolfgang. @bangerth I forgot to mention that this "sublayer net rotation removal" indeed can only be applied at the final velocity field for output. The velocity field resulted from this "sublayer net rotation removal" cannot go into the solver for next time step due to the "discontinuous velocity" problem you mention. Because the net rotation of lithosphere with respect to the underlying mantle is a natural result when there is lateral viscosity variation within the lithosphere, keeping it in the flow field is physically correct. I'm thinking of removing this only when I want to compare the final output surface flow with the plate motion data in no-net-rotation frame. Feel free to point out if you have more thoughts on this.

Can you give me a hint on how to define a class like the Rotation class used there that only returns a rotational velocity if $r&gt;R0$ and a zero rotation velocity otherwise? What line 463 to 465 does is outside of the above loop over each cell. Do I need to loop over each cell again to grab the shallower cells with $r&gt;R0$ to achieve this kind of special sublayer net rotation subtraction?

[bangerth]

The Rotation class looks like this:

    template <int dim>
    class Rotation : public TensorFunction<1,dim>
    {
      private:
        const Tensor<1,dim> axis;

      public:
        // Constructor for TensorFunction that takes cartesian direction (1,2, or 3)
        // and creates a solid body rotation around that axis.
        Rotation(const unsigned int a)
          :
          axis(Tensor<1,dim>(Point<dim>::unit_vector(a)))
        {}

        // Constructor for TensorFunction that takes an axis
        // and creates a solid body rotation around that axis.
        Rotation(const Tensor<1,dim> &rotation_axis)
          :
          axis(rotation_axis)
        {}

        virtual Tensor<1,dim> value (const Point<dim> &p) const
        {
          if ( dim == 2)
            return cross_product_2d(p);
          else
            return cross_product_3d(axis, p);
        }
    };

You only need to modify the last function as follows:

        virtual Tensor<1,dim> value (const Point<dim> &p) const
        {
          if (p.norm() > R0)
          {
            if ( dim == 2)
              return cross_product_2d(p);
            else
              return cross_product_3d(axis, p);
          }
          else
            return Tensor<1,dim>();   // zero vector
        }

[Shangxin-Liu]

Thanks, Wolfgang. Because I'd like to apply this "two-step" removal of the rotation (the first global removal like ASPECT currently does and the second top sublayer further removal), It seems that if I want both being applied in my simulation, I need to create another special "Rotation" class in null_space.cc file like you post.

One thing I'm still a little concerned and this links to the question I post on Rene's pull request to perform the rotation statistics @gassmoeller : For the velocity field used in the post processor (such as in dynamic topography and geoid to compute strain rate), is it the velocity solution already with the pure rotation removed or the original one without the removal of the pure rotation? If it is the one already with the rotation removed as null_space.cc does, adding this second further rotation removal on the top sublayer here is perhaps not what I want due to the discontinuous velocity issue that leads to the infinite strain rate. I may add another postprocessor to only calculate the rotation vector \omega from the top sublayer and apply this rotation subtraction by myself to the output shallower velocity field from ASPECT. Any suggestions? @bangerth

[bangerth]

See my answer here: #2212 (comment)