0
$\begingroup$

I am trying to calculate and plot the x-axis of a volume centroid for the region between zx-plane and a function squared $u(t,x,z)$ from time $t=15$ to $25$. The function is a numerical solution to a particular wave equation, codes of which I will attach at the end.

The code I used to obtain the position of this volume centroid is as follows:

centroid[t_] := 
  RegionPlot3D[
     0 <= z <= 75 && 0 <= y <= uwave[t, x, z]^2, {z, 0, 75}, {y, 
      0,.2}, {x, -1, 1}] // DiscretizeGraphics // RegionCentroid;

Take one particular time-instance, for example $t=18$. The plot of our function squared in 3D looks like this: enter image description here

and its code is

Plot3D[uwave[18, x, z]^2, {x, -1, 1}, {z, 0, 75}, PlotPoints -> 200, 
 PlotRange -> All, AxesLabel -> {"x", "z", "u"}]

Now clearly we would expect the centroid to have a x-axis slightly before $x=20$ as can be seen from the plot given. However, when I calculate it using the centroid function defined above, I find that

First[centroid[18]]

gives a value of 37.4976 which is clearly wrong. Where did I do wrong in the code?

For completeness, the numerical solution $uwave(t,x,z)$ is obtained by the following code which you may copy paste run, and get a result within 15 minutes.

Lx = 85;
Lz = 75;
Domain = Region[Rectangle[{-Lx, -Lz}, {Lx, Lz}]];

\[Rho]in = 5;
\[Rho]ex = 1;
w = 1;
p = 8;
\[Rho][x_] := (\[Rho]in - \[Rho]ex) (Sech[(x/w)^p])^2 + \[Rho]ex;
wx = 1;
wz = Sqrt[2];
u0saus[x_, z_] := x/wx*Exp[-(x/wx)^2 - (z/wz)^2];
VA[x_] := Sqrt[B0^2/\[Rho][x]]/Sqrt[B0^2/\[Rho][0]];
icX = u[0, x, z] == u0saus[x, z];
dicX = Derivative[1, 0, 0][u][0, x, z] == 0;
DX = DirichletCondition[u[t, x, z] == 0, True];
EqX = D[u[t, x, z], {t, 
      2}] - (VA[x])^2*(D[u[t, x, z], {x, 2}] + 
       D[u[t, x, z], {z, 2}]) == 0;
uwave = NDSolveValue[{EqX, icX, dicX, DX}, 
  u, {t, 0, 25}, {x, z} \[Element] Domain, Method -> {
    "PDEDiscretization" -> {"MethodOfLines",
      "SpatialDiscretization" -> {"FiniteElement", 
        "MeshOptions" -> {"MaxCellMeasure" -> 0.3}, 
        "InterpolationOrder" -> {u -> 2}}}}]
$\endgroup$
8
  • 1
    $\begingroup$ Please add definitions of B0 and u0saus (and any other missing values) to your code. $\endgroup$
    – bill s
    Apr 28, 2022 at 19:13
  • $\begingroup$ I have edited the original question to include a definition of u0saus. B0 is undefined because it will cancel out, so you may define whatever you like for B0. I have only included this thing B0 for physical completeness, a term which represents the background magnetic field. $\endgroup$
    – Rescy_
    Apr 28, 2022 at 19:16
  • $\begingroup$ You define u0sausX but use u0saus $\endgroup$ Apr 29, 2022 at 7:42
  • $\begingroup$ Further, there is something wrong with your function "centriod". The centroid should be a 3 dim point, not a scalar. $\endgroup$ Apr 29, 2022 at 8:22
  • $\begingroup$ Indeed, it is a point with 3 components as you said. I used first to extract its x-component $\endgroup$
    – Rescy_
    Apr 29, 2022 at 8:26

0

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.