nonisotropic kernel for diffusion

[Xiaojingtang1234]

non converged in both x and y direction the reference of kernel , paper "ADAPTIVE SMOOTHED PARTICLE HYDRODYNAMICS : METHODOLOGY. II. J. MICHAEL OWEN1 LLNL, " the tensor formation of 2d: for B-cubline kernel

for Wenland kernel: (my conclusion)

refer of SPH dissipative dynamics, https://sci-hub.hkvisa.net/10.1103/PhysRevE.67.026705 @Xiangyu-Hu

[Xiangyu-Hu]

What is your conclusion?

[Xiaojingtang1234]

What is your conclusion?

I do not make or try it out ....confusing..

nonisotropic transform tensor : Vec2d kernel_vector = (1.0, 1.0 / ratio) , ratio can be 2.0, 3.0, 4.0...... transform_angle(1.0) the displacemetnt in formal parameter is posi-posj

Mat2d AnisotropicKernel::getCoordinateTransformationTensorG(Vec2d kernel_vector, Vec2d transform_angle)

         {
    Mat2d scaling_tensor = Mat2d(1.0 / (this->h_* kernel_vector[0]), 0.0, 0.0, 1.0 / (this->h_* kernel_vector[1]));
    Mat2d rotation_tensor = getTransformationMatrixFromAngle(transform_angle);
    Mat2d G_kernel_coordinate = rotation_tensor * scaling_tensor * (~rotation_tensor);
    return G_kernel_coordinate;
        }

this is how i make eij in nonisotropic kernel

Vec2d AnisotropicKernel::geteij(const Real& distance, const Vec2d& displacement) const

{
    Vec2d transformed_displacement_ = this->h_ * transformed_tensor_2d_ * displacement;
    return  transformed_tensor_2d_ * transformed_displacement_ / (( transformed_displacement_).norm() + TinyReal);
}

this is how we get dw of nonisotropic ,

Real AnisotropicKernel::dW(const Real& r_ij, const Vec2d& displacement) const

{
    Vec2d transformed_displacement = transformed_tensor_2d_ * displacement;
    Real q = transformed_displacement.norm();
    return this->factor_dW_2D_ * this->dW_2D(q);
}

factor_dW2D= this->h * this->h det(transformed_tensor2d) this->FactorW2D(); FactorW2D() : factor_W2D = invh invh 7.0 / (4.0 * Pi);

correspond to this formulation

[Xiaojingtang1234]

Google test ; this is when Laplacian equals to 2, saturation = xx, while when saturation = yy, Laplacian around 0.01

TEST(test_anisotropic_kernel, test_rij)
{
 int y_num = 8;
Real PH = 1.0;
Real ratio_ = 4.0;
Real PL = ratio_* PH;

Vec2d scaling_vector = Vec2d(1.0, 1.0 / ratio_);
Real resolution_ref = PH / Real(y_num);
Real resolution_ref_large = ratio_ * resolution_ref;
int x_num = PL / resolution_ref_large;

AnisotropicKernel<KernelWendlandC2>  wendland(1.15 * resolution_ref_large, scaling_vector,  Vec2d(0.0));
Real second_order_rate = 0.0;
Real sum = 0.0;
Vec2d first_order_rate =Vec2d(0.0) ;
Vec2d center = Vec2d(resolution_ref_large * 4.0, resolution_ref *4.0);

Real V_ = resolution_ref * resolution_ref_large;

   for (int i = 0; i < (x_num + 1); i++)
   {
    for (int j = 0; j < (y_num + 1); j++)
    {
        Real x = i * resolution_ref_large;
        Real y = j * resolution_ref;

        Vec2d displacement =  center - Vec2d(x, y) ;
        Real distance_ = displacement.norm();

        Real  sarutration_ = 1.0 * x * x;
        Real  sarutration_center = 1.0 * center[0] * center[0];

        if (wendland.checkIfWithinCutOffRadius(displacement))
        {
            Vec2d  eij_dwij_V = wendland.geteij(distance_, displacement)* wendland.dW(distance_, displacement) * V_;

            sum += wendland.W(distance_, displacement)* V_;
            first_order_rate -= (sarutration_center - sarutration_)* eij_dwij_V;

            second_order_rate +=  2.0 * (sarutration_center - sarutration_)  * dot(displacement, wendland.geteij(distance_, displacement) )
                        /  (wendland.getrij(distance_, displacement) + TinyReal)
                    * wendland.dW(distance_, displacement) * V_ * 1.0 / (wendland.getrij(distance_, displacement)+ TinyReal);
        }       
    }

}

EXPECT_FLOAT_EQ(1.0, sum); 
EXPECT_EQ(2.0, second_order_rate);

}

[Xiaojingtang1234]

*I checked Paul W. Cleary∗ and Joseph J. Monaghan† paper, and use equation 16 with q = transformed_displacement = transformed_tensor2d displacement; it does not work right**. I did not find any other direct materials that have concluded this Laplacian in nonisotropic kernel

https://sci-hub.hkvisa.net/10.1006/jcph.1998.6118

[Xiangyu-Hu]

We may try to derive by find the relation factor between isotropic and anisotropic laplacian in its analytical form first. Then apply this relation to the numerical formulation.

[Xiangyu-Hu]

We may try to derive by find the relation factor between isotropic and anisotropic laplacian in its analytical form first. Then apply this relation to the numerical formulation.

[Xiangyu-Hu]

Could you make a branch with only the anisotropic kernel and the Google test?

[Xiaojingtang1234]

sure, i am making

[Xiaojingtang1234]

refer : https://iopscience.iop.org/article/10.1086/313100/pdf https://sci-hub.hkvisa.net/10.1007/s00466-004-0620-y

[Xiangyu-Hu]

Here is my derivation of the ASPH kernel approximations.

[Xiaojingtang1234]

ok i will see , and i already pushed the documents and google test on the branch xiaojing-anisotropic
in my library. (https://github.com/Xiaojingtang1234/SPHinXsys-xiaojing/tree/xiaojing-anisotropic/tests/user_examples)

I use another namespace to seperate the orginal kernel and my kernel. Some funtions are virtual and checkIfWithinCutOffRadius(Vec2d displacement) ; should be written into base_kernel.h

[Xiangyu-Hu]

Please have a check on the new formulation and we can discuss it.

[Xiaojingtang1234]

I checked, it seems not right.
I update this part code in my branch according to my understanding. if i understand correctly, 2|G|G:G is a scalar and the latter actually is the same with previous Laplacian discrization. it cannot modify the nonisotrpic nature in x and y direction.

https://github.com/Xiaojingtang1234/SPHinXsys-xiaojing/blob/0215d8ea286efc0fa622f9459241bf0c26b6c878/tests/user_examples/test_nonisotropic_kernel/test_nonisotropic_kernel.cpp#L68

[Xiangyu-Hu]

You can check if the derivation is wrong somewhere.

[Xiaojingtang1234]

could you update the new formulation? it seems not right in code. i did not find anything wrong. need to check the theory formulation again

[Xiangyu-Hu]

[Xiaojingtang1234]

.

[Xiangyu-Hu]

Yes. There is Something wrong.

[Xiaojingtang1234]

[Xiangyu-Hu]

There is the new formulation, it seems correct now.

[Xiaojingtang1234]

It seems something wrong . We still consider point (0,0) and (4, 1) for nonisotropic, Phi(r) = x^2+y^2, then we have phiij = Phi(0,0) - Phi_(4,1) = 0- (4x4 -1x1) =-17. but in isotropic case, it is phiij = Phi(0,0) - Phi_(1,1) = 0- (1x1 -1x1) =-2. Weired number 17 cannot be modified by G^2 : phi_ij Tensor;

[Xiangyu-Hu]

We do not need modify 17 to 2 for computing the laplacian.

[Xiangyu-Hu]

Update the newest formulation. We can have a try.

[Xiaojingtang1234]

I am more confused....i will try to have a try.. also i am checking how Laplacian operator works in polar coordinates, trying to find how different coordinates transform. Maybe we can get some useful information from polar coordinates

[Xiangyu-Hu]

Where is your confusing ?

[Xiaojingtang1234]

The result is bad. The profile is T =x^2

I checked my code implementation, when it is isotropic, it works right. i checked the tensor, it is right, dW is right. I put the code in the branch xiaojing_anisotropic kernel in SPHinxsys resposity, including the kernel. cpp and the google test.

[Xiaojingtang1234]

Note we refer this paper, this is the reason why there is a h in G . https://iopscience.iop.org/article/10.1086/313100/pdf

[Xiaojingtang1234]

ok, i checke the code, at line 42 in google test , we should use Vecd pos_i = Vecd(resolution_x * x_num / 2.0, resolution_y * y_num /2.0 ); // Particle i location to keep i is always at the center. Also i am wondering if we could first caculate the gradient T and then caculate the gradient gradient T. Hourglass may exist but it may work.

OK. i checked, result is ugly. the profile is not smooth at all. This does not work

[Xiangyu-Hu]

it is ok. We just use this formula and check the diffusion convergence.