Extract profile spherical velocity and data load by vivi235711 · Pull Request #386 · gamer-project/gamer

vivi235711 · 2024-12-19T08:02:35Z

Enhance GAMER_ExtractProfile with reading hybrid scheme output, improved accuracy, shell-average velocity, and spherical velocity output.

…rdson extrapolation into second order

…f velocity

…ocity output

… = 0 region

…tput volume and mass in gamer_extract_profile

…xtract_profile_sph_vel_and_data_load

…tract_profile_sph_vel_and_data_load

hsinhaoHHuang

@vivi235711
Thank you for providing this useful tool.
I have finished the PR review.
Please see my comments when you have time and let me know if they are not clear.
Thanks!

hsinhaoHHuang · 2025-01-08T09:32:06Z

+   if ( OutputSphere )
+   {
+                        sprintf( FileName[Var++], "%s", "AveVr"    );
+                        sprintf( FileName[Var++], "%s", "AveVtheta");
+                        sprintf( FileName[Var++], "%s", "AveVphi"  );
+                        sprintf( FileName[Var++], "%s", "AveWr"    );
+                        sprintf( FileName[Var++], "%s", "AveWtheta");
+                        sprintf( FileName[Var++], "%s", "AveWphi"  );
+   }


It seems OutputSphere only supports ELBDM currently. Would we plan to also support HYDRO? The coordinate transformation should also work for the hydro velocity.

Now can support hydro.
May I ask what hydro test problem you use to test this code?

I think we can use Hydro/Plummer test problem.

Rename variables, remove useless line, modify some if statements, and add usage

vivi235711 · 2025-02-12T07:11:56Z

@hsinhaoHHuang Thanks for your detailed reviewed! I've updated the code based on your comments. Please take a look when you have time.

@koarakawaii

Based on @koarakawaii 's comment

hyschive · 2026-02-10T03:04:03Z

@vivi235711 Please update this PR to synchronize it with the latest main.

@hsinhaoHHuang Please take another look at this PR (if needed) and approve it if you think it's ready to be merged.

hsinhaoHHuang

Apologies for my late response. I have reviewed it again and left some comments. Please let me know your thoughts about them when you have time. Thank you!

hsinhaoHHuang · 2026-03-24T01:33:49Z

         char Key[MaxString];
         sprintf( Key, "FieldLabel%02d", v );

+#        if ( MODEL == ELBDM)


Suggested change

# if ( MODEL == ELBDM)

# if ( MODEL == ELBDM )

hsinhaoHHuang · 2026-03-24T01:39:01Z

+      // Transfer dens and phase to real and imag
+#     if ( MODEL == ELBDM )
+      if ( ELBDM_Scheme == 2 )
+      {  


Suggested change

{

{

hsinhaoHHuang · 2026-03-24T02:08:41Z

+                  if ( ELBDM_RichExtrap )
+                  {
+//                   Richardson extrapolation 2 order
+                     GradD[0] = (    4*_2dh*( Field[p][DENS][k  ][j  ][ip ] - Field[p][DENS][k  ][j  ][im ] ) 


Suggested change

GradD[0] = ( 4*_2dh*( Field[p][DENS][k ][j ][ip ] - Field[p][DENS][k ][j ][im ] )

GradD[0] = ( 4*_2dh*( Field[p][DENS][k ][j ][ip ] - Field[p][DENS][k ][j ][im ] )

hsinhaoHHuang · 2026-03-24T02:30:51Z

+   if ( average_dens  != NULL )  delete [] average_dens;
+   cout << " Bulk CoM velocity = ("<< CMV[0]<< ", "<< CMV[1] << ", "<<CMV[2]<< ")" << endl;
+
+   cout << "Remove the motion of center-of-mass by phase shift" << endl;


I think this message cout << "Remove the motion of center-of-mass by phase shift" << endl; can be put inside # elif ( MODEL == ELBDM ).

As for MODEL == HYDRO, we should be able to do something like this to remove the center-of-mass velocity:

Dens = amr.patch[lv][PID]->fluid[DENS][k][j][i]; VelX = amr.patch[lv][PID]->fluid[MOMX][k][j][i] / Dens; VelY = amr.patch[lv][PID]->fluid[MOMY][k][j][i] / Dens; VelZ = amr.patch[lv][PID]->fluid[MOMZ][k][j][i] / Dens; amr.patch[lv][PID]->fluid[MOMX][k][j][i] = (VelX - CMV[0]) * Dens; amr.patch[lv][PID]->fluid[MOMX][k][j][i] = (VelY - CMV[1]) * Dens; amr.patch[lv][PID]->fluid[MOMX][k][j][i] = (VelZ - CMV[2]) * Dens;

hsinhaoHHuang · 2026-03-24T02:31:52Z

+#   else
+#   error : ERROR : unsupported MODEL !!
+#   endif // MODEL
+} // FUNCTION : Remove_CMVel 


Suggested change

} // FUNCTION : Remove_CMVel

} // FUNCTION : Remove_CMVel

hsinhaoHHuang · 2026-03-24T04:25:04Z

 double     *ELBDM_RhoUr2      = NULL;     // Rho*(vr^2 + wr^2) in the virial surface terms
 double     *ELBDM_dRho_dr     = NULL;     // dRho/dr in the virial surface terms
 double     *ELBDM_LapRho      = NULL;     // Lap(Rho) in the virial surface terms
+bool        ELBDM_RichExtrap  = true;     // Richardson extrapolation for the velocity


I feel this name may be confusing. I think we didn't really do Richardson extrapolation. What this option means in this code is actually to use a 5-point stencil central difference method to compute the gradient and Laplacian. Although the adopted formula can be derived by Richardson extrapolation, it seems to me that "Richardson extrapolation" is a more general term referring to estimating the value using two different resolutions (not necessarily dh and 2*dh).
Do I understand it correctly?

I would think a name like ELBDM_5PointDerivatives (just an example) would at least be better to convey that the derivatives will be computed by a 5-point stencil instead of a 3-point stencil.

hsinhaoHHuang · 2026-03-24T04:26:08Z

 bool        GetAvePot         = false;    // calculate the spherically symmetric gravitational potential
 double      NewtonG           = WRONG;    // gravitational constant (must be set if it cannot be determined from the input data)
+bool        RemoveCMV         = false;    // remove the center-of-mass velocity 
+bool        OutputSphVel      = false;    // output the sphere coordinate (r, theta, phi) velocity


Suggested change

bool OutputSphVel = false; // output the sphere coordinate (r, theta, phi) velocity

bool OutputSphVel = false; // output the spherical coordinate (r, theta, phi) velocity

hsinhaoHHuang · 2026-03-24T05:33:37Z

+
+#                 if   ( MODEL == HYDRO )
+//                evaluate the values on the shell
+#                 warning : WAIT HYDRO !!!


I think the calculation is straightforward in HYDRO:

Dens = Field[p][DENS][k][j][i]; MomX = Field[p][MOMX][k][j][i]; MomY = Field[p][MOMY][k][j][i]; MomZ = Field[p][MOMZ][k][j][i]; CMV[0] += (double)(dv*MomX); CMV[1] += (double)(dv*MomY); CMV[2] += (double)(dv*MomZ); massAll += (double)(dv*Dens);

hsinhaoHHuang · 2026-03-24T05:57:06Z

+                        GradD[0] = (    4*_2dh*( Field[p][DENS][k  ][j  ][ip ] - Field[p][DENS][k  ][j  ][im ] ) 
+                                    - 0.5*_2dh*( Field[p][DENS][k  ][j  ][ipp] - Field[p][DENS][k  ][j  ][imm] ))/3;
+                        GradD[1] = (    4*_2dh*( Field[p][DENS][k  ][jp ][i  ] - Field[p][DENS][k  ][jm ][i  ] )
+                                    - 0.5*_2dh*( Field[p][DENS][k  ][jpp][i  ] - Field[p][DENS][k  ][jmm][i  ] ))/3;
+                        GradD[2] = (    4*_2dh*( Field[p][DENS][kp ][j  ][i  ] - Field[p][DENS][km ][j  ][i  ] )
+                                    - 0.5*_2dh*( Field[p][DENS][kpp][j  ][i  ] - Field[p][DENS][kmm][j  ][i  ] ))/3;


Do you think it would be cleaner and easier to maintain if we wrap the calculation as a function?

void grad( real* field_out, const real* field_in, const int i, const int j, const int k, const double dh) { if ( ELBDM_RichExtrap ) { field_out[0] = ( 4*_2dh*( field_in[k ][j ][i+1] - field_in[k ][j ][i-1] ) -0.5*_2dh*( field_in[k ][j ][i+2] - field_in[k ][j ][i-2] ))/3; field_out[1] = ( 4*_2dh*( field_in[k ][j+1][i ] - field_in[k ][j-1][i ] ) -0.5*_2dh*( field_in[k ][j+2][i ] - field_in[k ][j-2][i ] ))/3; field_out[2] = ( 4*_2dh*( field_in[k+1][j ][i ] - field_in[k-1][j ][i ] ) -0.5*_2dh*( field_in[k+2][j ][i ] - field_in[k-2][j ][i ] ))/3; } else { field_out[0] = _2dh*( field_in[k ][j ][i+1] - field_out[k ][j ][i-1] ); field_out[1] = _2dh*( field_in[k ][j+1][i ] - field_out[k ][j-1][i ] ); field_out[2] = _2dh*( field_in[k+1][j ][i ] - field_out[k-1][j ][i ] ); } } real laplacian( const real* field_in, const int i, const int j, const int k, const double dh) { // ... return laplacian; }

Then, we can invoke it simply by

grad( GradD, Field[p][DENS], i, j, k, dh); grad( GradR, Field[p][REAL], i, j, k, dh); grad( GradI, Field[p][IMAG], i, j, k, dh); LapR = laplacian( Field[p][REAL], i, j, k, dh); LapI = laplacian( Field[p][IMAG], i, j, k, dh);

hsinhaoHHuang · 2026-03-24T06:00:24Z

+   if ( OutputSphere )
+   {
+                        sprintf( FileName[Var++], "%s", "AveVr"    );
+                        sprintf( FileName[Var++], "%s", "AveVtheta");
+                        sprintf( FileName[Var++], "%s", "AveVphi"  );
+                        sprintf( FileName[Var++], "%s", "AveWr"    );
+                        sprintf( FileName[Var++], "%s", "AveWtheta");
+                        sprintf( FileName[Var++], "%s", "AveWphi"  );
+   }


I think we can use Hydro/Plummer test problem.

hyschive · 2026-04-14T04:27:41Z

@vivi235711

The compilation fails for MODEL=HYDRO
Please synchronize with the latest main branch.

koarakawaii · 2026-04-30T10:35:35Z

When running the code with MODEL = ELBDM with -V option on, the RMS for AveDens_* is different when -O option is on/off (probably related to this line, which should not touch the density data when -O option is off, but seems to still touch it because NOut seems to be 12 in this case).

vivi235711 added 13 commits September 12, 2023 10:10

Enhance 'GAMER_ExtractProfile' to read hybrid scheme data

a7e5179

Improve the differential accuracy in GAMER_ExtractProfile using Richa…

70bd31d

…rdson extrapolation into second order

Remove the center of mass velocity for evaluating the shell-average o…

b8e2b7c

…f velocity

Add option 'O' to 'GAMER_ExtractProfile' for generating spherical vel…

ff02372

…ocity output

Enable differential accuracy as an option in GAMER_ExtractProfile

73abe31

Correct CMV calculation, changed sin&cos to SIN&COS, and skip density…

cad6826

… = 0 region

Remove CMV using phase shift

b2b39d5

Change to wave function derivative to calculate internal velocity; Ou…

bcfdb80

…tput volume and mass in gamer_extract_profile

Correct velocity calculation

c475678

Remove center of mass velocity as an option 'C'

67aae0c

Fix error of unsupported model

dc0f00b

Merge branch 'psidm' of https://github.com/hyschive/gamer-fork into e…

7d26a85

…xtract_profile_sph_vel_and_data_load

Merge branch 'main' of https://github.com/gamer-project/gamer into ex…

9003366

…tract_profile_sph_vel_and_data_load

hyschive requested a review from hsinhaoHHuang December 20, 2024 13:16

hyschive assigned hsinhaoHHuang Dec 20, 2024

hyschive added enhancement tool Useful analysis and visualization tools labels Dec 20, 2024

hsinhaoHHuang suggested changes Jan 8, 2025

View reviewed changes

vivi235711 added 4 commits February 8, 2025 15:30

Style: Align code

372ff76

Fix: HDF5 loading with ELBDM_Scheme

bed46c1

Refactor code based on @hsinhaoHHuang's comment

b4d08b7

Rename variables, remove useless line, modify some if statements, and add usage

Support: OutputSphVel in Hydro

59ad25d

vivi235711 added 2 commits March 6, 2025 14:43

Change RMS velocity to mass-weighted

f9262bf

Turn on ELBDM_GetVir when -O is on in ELBDM

239d21c

Based on @koarakawaii 's comment

hsinhaoHHuang suggested changes Mar 24, 2026

View reviewed changes

	GradD[0] = ( 4_2dh( Field[p][DENS][k ][j ][ip ] - Field[p][DENS][k ][j ][im ] )
	GradD[0] = ( 4_2dh( Field[p][DENS][k ][j ][ip ] - Field[p][DENS][k ][j ][im ] )

	bool OutputSphVel = false; // output the sphere coordinate (r, theta, phi) velocity
	bool OutputSphVel = false; // output the spherical coordinate (r, theta, phi) velocity

Conversation

vivi235711 commented Dec 19, 2024

Uh oh!

hsinhaoHHuang left a comment

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

vivi235711 commented Feb 12, 2025

Uh oh!

hyschive commented Feb 10, 2026

Uh oh!

hsinhaoHHuang left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

hyschive commented Apr 14, 2026

Uh oh!

koarakawaii commented Apr 30, 2026

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants